Analyzing wellbore stability in chemically-active anisotropic formations under thermal, hydraulic, mechanical and chemical loadings

Accepted Manuscript Analyzing wellbore stability in chemically-active anisotropic formations under thermal, hydraulic, mechanical and chemical loadings Majed F. Kanfar, Z. Chen, S.S. Rahman PII:

S1875-5100(17)30051-3

DOI:

10.1016/j.jngse.2017.02.006

Reference:

JNGSE 2059

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 31 March 2016 Revised Date:

29 January 2017

Accepted Date: 1 February 2017

Please cite this article as: Kanfar, M.F., Chen, Z., Rahman, S.S., Analyzing wellbore stability in chemically-active anisotropic formations under thermal, hydraulic, mechanical and chemical loadings, Journal of Natural Gas Science & Engineering (2017), doi: 10.1016/j.jngse.2017.02.006. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Title: Analyzing Wellbore Stability in Chemically-Active Anisotropic Formations under

2

Thermal, Hydraulic, Mechanical and Chemical Loadings

3

Authors Order: Majed F. Kanfar *, Z. Chen, S.S. Rahman,

4

School of Petroleum Engineering, University of New South Wales, Australia

5

*Corresponding author: Email address: [email protected] (Majed Kanfar)

6

Keywords: anisotropy; porochemothermoelasticity; THMC; thermo-chemo-hydro-

7

mechanical; wellbore stability; shale

8

Abstract

RI PT

1

Several factors influence the stress state in subsurface rocks. In addition to pore pressure

10

and far-field in-situ stresses, thermal and chemical gradients have substantial bearings on the

11

in-situ stress state during and after drilling operations. Drilling inclined boreholes through

12

laminated formations such as shaley sands presents several challenges. Due to their laminated

13

nature, shales demonstrate high degree of elastic anisotropy, and are often chemically-

14

reactive to drilling fluid. The majority of models utilized in the industry assume a

15

homogeneous isotropic elastic static model that fails to give an accurate depiction of the

16

observed borehole failure.

17

porochemothermoelasticity are developed to simulate the drilling of an inclined borehole

18

problem in a transversely isotropic rock. Using the developed constitutive and transport

19

equations, a finite element method based numerical model is constructed to estimate the pore

20

pressure, temperature, solute concentration and stress distribution. Nonlinear conductive-

21

convective heat transfer and diffusive-advective solute transfer models are considered in this

22

work to address both high and low permeability formations. Finally, the model is used to

23

assess time-dependent wellbore stability during and after drilling. The novel pseudo-3D

24

analysis developed in this work is advantageous for real-time operations due to its

25

computational speed and stability.

M AN U

TE D

In this work, governing equations for anisotropic

EP

AC C

26

SC

9

ACCEPTED MANUSCRIPT 27

1. Introduction With skyrocketing demand for energy supply, the oil and gas industry has been

29

aggressively pursuing more ambitious ventures. Recent technological developments enabled

30

the exploitation of oil and gas prospects previously deemed unattainable, such as

31

unconventional shale reservoirs and deep water horizons. Drilling applications extend outside

32

of the fossil fuel industry, including for CO2 sequestration operations and harvesting

33

alternative sources of energy such as geothermal. These recent complex applications require

34

vigorous knowledge of the reservoir fluid and rock properties to accurately simulate coupled

35

dynamic problems. These rocks are chemically active and the physico-chemical nature should

36

be carefully analyzed. The borehole-reservoir system in most cases is a non-isothermal one.

37

Bespoke simulators are crucial to account for the solute concentration and temperature

38

gradient, and their effect on pore pressure and in-situ stresses. These coupled phenomena

39

produce complex unpredictable environments during drilling operations if they are not

40

properly modeled.

M AN U

SC

RI PT

28

Most simulators that couple fluid flow with geomechanics assume an isotropic rock

42

despite the fact that only 10% of subsurface formations exhibit true isotropic behavior (Ong,

43

1994). More than 30% of rocks classified as anisotropic have an anisotropy ratio of 1.5 for

44

Young’s modulus (Batugin and Nirenburg, 1972). The majority of shales are deemed

45

anisotropic. Shales make up approximately 75% of the drilled section and are the culprit for

46

90% of borehole instability problems (Mody and Hale, 1993). Several additional factors

47

affect shales behavior under stress such as permeability variation under different confining

48

pressures (Vishal et al., 2015, 2013), and its brittleness which is dependent on its mechanical

49

properties (Zhang et al., 2016).

EP

TE D

41

It is crucial to model the dynamic nature of a drilling problem during and after drilling.

51

Static elastic models (Aadnoy and Chenevert, 1987; McLean and Addis, 1990) fall short as

52

they are not time-dependent. Modelling the transient pore pressure behavior and its influence

53

on effective stresses were addressed by several authors using poroelasticity (Abousleiman et

54

al., 1995; Detournay and Cheng, 1988). The generalized anisotropic poroelasticity was

55

developed by Biot (1955). The majority of subsurface rocks exhibit a transversely isotropic

56

(TI) behavior caused by compaction and gravitational forces (Biot, 1955). Analytical solution

57

for pore pressure and stress distribution around a borehole in a poroelastic TI material was

58

developed by Abousleiman et al. (1995). They assume that the plane of isotropy of the TI

59

material is always perpendicular to the borehole axis.

AC C

50

ACCEPTED MANUSCRIPT Several publications have extended Biot's poroelasticity to a fully coupled

61

porothermoelastic model (Bear and Corapcioglu, 1981; Coussy, 1989; Kurashige, 1989;

62

Schiffman, 1971). A fully coupled isotropic conductive-convective porothermoelastic heat

63

and fluid flow model for boreholes with constant pressure and temperature was developed by

64

Wang and Dusseault (2003). A special formulation for isotropic porothermoelasticity for

65

swelling shale has been developed by Ghassemi and Diek (2002) to account for thermal

66

osmosis. Abousleiman and Ekbote (2005) developed an analytical solution for inclined

67

boreholes in transversely isotropic media. They assume that the plane of isotropy of the TI

68

material is always perpendicular to the borehole axis, and heat transfer is by conduction only.

69

Kanfar et al. (2016b) developed a fully coupled anisotropic porothermoelastic model using

70

the generalized plane strain (GPS) assumption that removes the abovementioned assumption

71

made by Abousleiman and Ekbote (2005).

SC

RI PT

60

Increased borehole collapse problems are observed during drilling in shale formations.

73

This is attributed to the chemically active nature of clays that are ubiquitous in shales and

74

shaley formations (Mody and Hale, 1993). Coupling this chemical phenomenon with Biot's

75

poroelastic model and its effect on chemically active rocks was first discussed by Sherwood

76

(1993). A chemo-mechanical model for wellbore stability that accounts for osmotic and

77

poroelastic effects was developed by Ghassemi and Wolfe (1999). Specialized physico-

78

chemical analytical solution for swelling shales were developed for pore pressure, stresses

79

and solute mass fraction around a vertical isotropic well (Ghassemi and Diek, 2003). Ekbote

80

and Abousleiman (2006) derived analytical solutions for porochemoelasticity in inclined

81

boreholes drilled in TI formations by assuming that the plane of isotropy is always

82

perpendicular to the borehole axis. The aforementioned authors assume solute transfer by

83

diffusion governed by Fick’s law only. A diffusion-advection solute transfer was developed

84

by Bader and Kooi (2005) that also accounts for membrane efficiency. The electrochemical

85

effect on the chemoelastic model was investigated in literature (Nguyen and Abousleiman,

86

2010; Tran and Abousleiman, 2013) and it is out of the scope of this work. Kanfar et al.

87

(2016a) investigated time-dependent wellbore stability by developing a porochemoelastic

88

model for anisotropic rocks that accounts for Fick’s diffusion and solute transfer by

89

advection. The aforementioned models do not consider the effect of the thermal gradient.

AC C

EP

TE D

M AN U

72

90

Due to the complexity of the problem, few authors have coupled the influence of

91

temperature gradient and chemical potential with Biot's poroelasticity (Ekbote and

92

Abousleiman, 2005; Ghassemi et al., 2009; Zhou and Ghassemi, 2009). Using the thermo-

93

chemo-hydro-mechanical (THMC) model, a numerical experiment that investigates CO2

ACCEPTED MANUSCRIPT sequestration by saturating it in water, and injecting the mixture into a carbonate aquifer was

95

conducted by Yin et al. (2011). Yasuhara et al. (2016) investigated the evolution of rock

96

permeability under THMC coupling. While these numerical and analytical models account

97

for the thermo-chemo-hydro-mechanical effects, none of which account for complete

98

material anisotropy to properly simulate case scenarios routinely encountered in field

99

applications. The effect of solute advection and thermal convection are also overlooked.

RI PT

94

In this work, a numerical model is developed to simulate fully coupled

101

porochemothermoelasticity for anisotropic formations. First, constitutive and transport

102

equations are developed. A numerical model using the finite element method (FEM) is built

103

using the governing equations. The GPS assumption is used for full three dimensional

104

analysis and stability of the numerical computation (Kanfar et al., 2015a). An inclined

105

borehole problem drilled in in shaley sand TI formation is mainly discussed in this work, but

106

the model can be applied to different facets of rock mechanics. The model is not exclusive to

107

shales or shaley sand, but can be generalized to other rock types as long as the material

108

properties are representative. It is assumed that the rock is homogeneous anisotropic and

109

fully-saturated with formation water. The model can be extended to other single phase liquid

110

fluids such as dead oil. Finally, time-dependent wellbore stability after drilling is analyzed

111

using the developed multi-physics anisotropic chemo-thermo-hydro-mechanical model.

112

2. Governing Equations for Porochemothermoelasticity

113

2.1. Constitutive Equations

EP

TE D

M AN U

SC

100

The constitutive equations for porochemothermoelasticity are governed by the stress-

115

strain relation, and the variation of the species content that comprise the solution of the

116

chemical system. The classical stress-strain relation is modified to account for thermal

117

expansion, and change in stresses due to the solute mass fraction perturbation. The variation

118

of fluid content per unit referential volume ( ) is also modified to incorporate indirect flows

119

of solute (

120

The constitutive equations are given by (Ekbote and Abousleiman, 2005; Ghassemi et al.,

121

2009):

AC C

114

) and diluent (

=

=

) by chemical osmosis, and fluid mixture thermal expansion.

+

− 1

+

+

∆

+ Ξ ∆

− Λ ∆

− Ω ∆

(1) (2)

ACCEPTED MANUSCRIPT 122

where

is the total stress,

is the total strain,

is the temperature,

is Biot’s modulus and

123

is the pore pressure. A repeated subscript refers to the Einstein summation. The elastic,

124

thermal and chemical coefficients presented in the constitutive equations (1) and (2) are

125

defined as follows:

̅

=$ −

=

=

# = +1 1 1 Ξ=

̅

(

.

(3)

3&' (

)*

0 0

/1 −

̅

̅

0

1 − #2

0-.

= #. + 1 − #2- / Λ = 3' + # /

4(

̅

̅

0

̅ !

TE D

4( Ω = 3 ! + #. + 1 − #2- / 3! =

3' = '

+1

!

'

−

'

!

(10) (11) (12) (13)

where

127

coefficient (chemoelastic coupling term),

128

universal gas constant,

129

temperature and 5 is the rock porosity. It is assumed that the swelling parameters for the

'

is the fluid density,

(

126

EP

is Biot’s effective stress coefficient vector,

(6)

(9)

0

25

(5)

(8)

0

̅

(4)

(7)

M AN U

=

∙#

̅ !

RI PT

−

SC

=

is the swelling

is the molar mass of the solute, R is the

is the solid thermal expansion coefficient vector,

!

is the fluid

131

solute (

132

value and the transpose of the matrix, respectively. The drained elastic tensor ( ) and the

133

compliance tensor ( ) in the directions of the orthotropic material are defined as:

AC C

130

thermal expansion coefficient, 4( is the fluid specific entropy measured at the average system ) and diluent (

) are equal. The overbar and the superscript

denote the average

ACCEPTED MANUSCRIPT

134 135

where ;, ? and B are the Young’s modulus, Poisson’s ratio and shear modulus, respectively,

and the subscripts refer to the coordinate systems described in Fig. 1. A detailed explanation

136

on the parameters for anisotropic poroelasticity is discussed by Cheng (1997). The strain can

137

be related to the solid displacement (F) as follows:

1 = HF , + F , J 2

0 0 0 0 1⁄B
0 E 0 D 0 D 0 D D 0 D 1⁄B<@= C

67

(14)

RI PT

0 0 0 1⁄B@A= 0 0

SC

−?A<= ⁄;A= −?A@= ⁄;A= 1⁄;A= 0 0 0

67

1⁄;<= : 9 −?@<= ⁄;@= 9−?A<= ⁄;A= =9 0 9 0 9 0 8

(15)

where a subscript followed by a comma denotes the spatial differentiation.

M AN U

138

−?@<= ⁄;@= 1⁄;@= −?A@= ⁄;A= 0 0 0

=

139

It is important to note that all material property tensors and stress tensor are rotated to the

140

borehole coordinate systems. This is achieved by defining three coordinate systems that are

141

referenced to the North-East-Vertical (NEV) global coordinates: 1. The far-field stresses

143 144

coordinates 1K' , L' , M' 2 that are controlled by angles N' and O' , 2. The intrinsic rock

1KP , LP , MP 2 that are controlled by angles NP and OP and 3. The

borehole coordinates 1KQ , LQ , MQ 2 that are controlled by angles NQ and OQ . Figure 1

properties coordinates

TE D

142

147

borehole 1R2 are referenced to the bottom of the hole (BOH), which always coincides with

148

the KQ axis. A schematic of the finite element method mesh is shown in Fig. 2, where the red lines indicate the cross-sectional plane of the FEM mesh. Property rotation is discussed

149

thoroughly in literature (Aadnoy and Chenevert, 1987; Amadei, 1983; Lekhnitskii, 1981).

EP

146

summarizes the coordinate systems and the associated rotation angles. Angles around the

AC C

145

150

The model formulation is simplified for a binary solution consisting of a single solute

151

(CS) and a diluent (CD) or solvent. It can be easily modified to include multiple solute species

152

as follows:

=S T7

153

= 1−

where k is the number of solute species.

(16)

ACCEPTED MANUSCRIPT 2.2. Transport Equations

155 156

A general conduction law or transport equation (Onsager’s relation) can be written as (Groot, 1951):

U = V K,

(17)

158

where U is the flux, V

159

law, chemical osmosis and thermal osmosis is given by:

161

162

gradient. The transport equation for fluid flux due to hydraulic potential governed by Darcy’s U = − ̅! WX !

,

− V X

− Y.

,

,

Z

(18)

where Y . is the thermal osmosis coefficient tensor, X is mobility tensor and V is the osmotic

pressure and are given by:

X =

V =

Y [

SC

160

is the phenomenological coefficient tensor and K, is the potential

) * !̅ ℜ ̅ ̅

M AN U

157

RI PT

154

(19) (20)

where Y is the rock permeability tensor, [ is the fluid viscosity measured at in-situ conditions

164

and ℜ is the reflection or membrane efficiency coefficient. The reflection coefficient can be

165

Kooi, 2005). Similarly, the solute transport equation due to solute diffusion governed by a

166

modified Fick’s law, thermal filtration and fluid advection is given by:

169

170

TE D

168

U ' = − ̅! 11 − ℜ2W



,

+

.



,

! ̅ −U .

Z

(21)

Fick’s law is modified by introducing the multiplier 11 − ℜ2 to account for the membrane

where

is the solute diffusion coefficient tensor and

EP

167

defined as the ratio of the hydraulic gradient and the osmotic pressure gradient (Bader and

is the thermal filtration coefficient.

efficiency (Bader and Kooi, 2005). Finally, the heat transfer equation is given by: U] = −^Q

,

+ ̅! ! U

where ^Q is the thermal conductivity and

AC C

163

!

!

(22)

is the fluid specific heat. The first term in the

171

heat transport equation accounts for bulk conduction by Fourier’s law, and the second for

172

heat transfer by fluid convection. The thermal conductivity (^Q ) and bulk heat capacity

173

(

Q Q)

are defined as:

and

^Q = 5^! + 11 − 52^'

Q Q

= 5 !̅

!

+ 11 − 52

' '

(23)

are the density and specific heat, respectively. The subscripts _, 4 and ` refer

174

where

175

to the bulk, solid matrix and fluid, respectively.

ACCEPTED MANUSCRIPT

177 178

179

2.3. Conservation Laws The conservation laws presented in this section ignore all body forces such as gravity and inertial forces. The momentum balance equation yields the equilibrium equation: ,

=0

(24)

A mass balance for the fluid and solute are given by: 1a ⁄ab2 + U !, = 0

5c1a

⁄ab2 + U ', = 0

RI PT

176

(25) (26)

181

where c is the solute retardation coefficient, which depends on the rock susceptibility to

182

assigned to this coefficient to neglect its effect on the analysis. Finally, the energy balance

183

equation is given by: Q Q

184

where

185

2.4. Field Equations

Q Q 1a

is the bulk heat capacity.

SC

absorb the chemical species and cation exchange (Bader and Kooi, 2005). A value of 1 can be

⁄ab2 + U], = 0

M AN U

180

(27)

The presented constitutive and transport equations are substituted into the conservation

187

laws to yield the field equations. The Navier-type equation for displacement in the x-y plane,

188

using the generalized plane strain assumption, is obtained by substituting the stress-strain

189

relation into the momentum balance equation as follows:

1 W1 2

77

−

7d 2F ,

TE D

186

+1

77

+

7d 2F ,

Z−

,

+

,

−Λ

,

=0

(28)

Incorporating the constitutive equation for the variation of fluid content, and the fluid

191

transport equation into the mass balance yields the local continuity equation:

1

+

AC C

192

1∇Ff 2 +

EP

190

5c f + 11 − ℜ2 1U '

f + Ξ f − Ω f − X∇d + V ∇d

+ Y . ∇d = 0

̅ ∇d = 0

(29)

Similarly, the field equation for the solute mass fraction is given by: 2, − 11 − ℜ2

− 11 − ℜ2

∇d

.

(30)

193

Substituting the heat flux equation into the energy balance relation yields the field equation

194

for the temperature:

195 196 197

Q Q

f − ^Q ∇d + ̅!

! HU

!

J, = 0

(31)

It is important to note the nonlinearity in equations (30) and (31) due to their dependence on

the fluid flux term 1U ! 2, which is responsible for the solute transfer by advection and heat

flux by convection, respectively. In this work, a homogenous anisotropic rock is assumed.

ACCEPTED MANUSCRIPT However, heterogeneity in rock properties such as permeability and porosity can be easily

199

incorporated into the model by varying them in each element of the FEM mesh. All rock

200

pores are assumed to be interconnected and fully saturated. In addition, no mass loss from the

201

rock matrix is allowed by the model by means such as sanding. It is also assumed that the

202

effluent from the borehole is fully miscible with the reservoir fluid.

203

3. Finite Element Formulae

RI PT

198

This section presents the anisotropic porochemothermoelastic formulae using the finite

205

element method. The generalized plane strain assumption will be used to allow for 3D

206

analysis and incorporation of material anisotropy (Amadei, 1983; Lekhnitskii, 1981). This

207

pseudo-3D analysis utilizes a 2D mesh, which is computationally far less intensive and

208

consumes significantly lower memory (Kaliakin, 2002). The classical plane strain assumption

209

only allows a 2D stress tensor, where the axial stress is always perpendicular to the mesh

210

plane and no anti-plane shear stresses can be applied. Moreover, the classical plane strain

211

assumption does not permit the use of material anisotropy. The GPS assumption and how

212

material anisotropy is incorporated was discussed thoroughly in a recent publication by the

213

authors (Kanfar et al., 2016b, 2015a). Quadratic eight node quadrilateral elements are used in

214

the examples presented in this work, and a simplified coarse prototype is shown in Fig. 3. An

215

inclined borehole problem drilled in a transversely isotropic shaley sand that is fully-saturated

216

with water will be mainly analyzed in this work. The unknown variables of the FEM system,

218

M AN U

TE D

F, ,

and , are given by:

F = gh Fi , = gj i,

= gk l m , = g. n

(32)

where gh , gj , gk l and g. are the shape functions of solid displacement, pore pressure,

EP

217

SC

204

solute mass fraction and temperature, respectively, and the tilde denotes the nodal values.

220

Using the Galerkin weighted residuals method for spatial discretization and the implicit finite

221

difference method for temporal discretization, the finite element formulation of the field

222

equations is given by:

AC C

219

oh Fif + oj if + ok m f + o. nf = `f

o.j Fif + ℙj if + ℙjq i + ℙk m f + ℙkq m + ℙ. nf + ℙ.q n = 0 ℂk m f + ℂkq m + ℂ.q n = 0 s . nf + s .q n = 0

(33)

ACCEPTED MANUSCRIPT

225

oh : . 9oj 9 0 9 8 0

oj

−ok

o.

in the appendix. The incremental form of the system of equations is given by: −Hℙj + ∆b ℙjq J −1ℙk − ∆b ℙkq 2 0 −1ℂk + ∆b ℂkq 2 0 0

E ∆Fi 1ℙ. − ∆b ℙ.q 2 D ∆ i D t∆ m u −∆b ℂ.q D −1s . + ∆b s .q 2C ∆ n

∆` : E m n 9∆b ℙjq i1vw 672 − ∆b ℙkq 1vw 672 + ∆b ℙ.q 1vw 672 D =9 D ∆b ℂkq m 1vw 672 + ∆b ℂ.q n1vw 672 9 D ∆b s .q n1v 672 8 C w

RI PT

224

The matrices oh , oj , ok , o. , ℙj , ℙjq , ℙk , ℙkq , ℙ. , ℙ.q , ℂk , ℂkq , ℂ.q , s . and s .q are defined

(34)

where ∆` is the change of the external traction forces, ∆b = bx − bx67 , bx and bx67 indicate

SC

223

the current and the previous time steps, respectively. First, the temperature equation is solved

227

separately. Then, the displacement, pore pressure and solute mass fraction are solved

228

simultaneously. However, the system of equation will then become nonsymmetrical. The

229

and s .q contain the fluid discharge term which drives the advection and convection

M AN U

226

230

symmetry is restored by employing a partitioned solution (Schrefler, 1985). The matrices ℂkq

231

contributions, respectively. An iterative scheme is developed to handle the nonlinearity in

232

equation (34). The solution converges when the following criteria are met: y∆ ivw − ∆ i1vw 672 y y∆ i1vw 672 y

< {,

y∆ mvw − ∆ m1vw 672 y y∆ nvw − ∆ n1vw 672 y < {, <{ y∆ n1v 672 y y∆ m1v 672 y

TE D

233

w

w

(35)

235

where { is the maximum error allowed by the iterative solver. In this work, the maximum

236

Along the borehole boundary, a constant pressure, solute concentration and temperature

237

are applied. Similarly, the in-situ far-field stresses are applied to the reservoir outer boundary

238

as traction forces. The FEM program boundary conditions (BC) are summarized in Table 1. It

239

is assumed throughout this paper that the borehole wall is permeable. An impermeable

240

borehole wall boundary condition can be easily implemented for cases where a “perfect”

241

filter cake is developed. In such cases, the pore pressure and solute concentration are not

242

influenced by the borehole boundary conditions (Fjaer, 2008). Moreover, the wellbore BC

243

will only affect the traction force acting on the borehole wall and heat transfer by conduction.

244

An impermeable BC was previously discussed in literature (Cui et al., 1998; Fjaer, 2008) and

245

it is out of the scope of the paper, though it can be easily implemented in the FEM code.

EP

error is set to 0.5%.

AC C

234

ACCEPTED MANUSCRIPT 246

4. Verification of Numerical Model In this section, the anisotropic porochemothermoelastic numerical model is validated

248

using published analytical solutions. The GPS assumption, anisotropic elasticity and

249

poroelasticity models were benchmarked in a previous publication by the authors (Kanfar et

250

al., 2015a). The modified Terzaghi or Biot’s effective stress concept (Nur and Byerlee, 1971)

251

is used throughout, and given by:

252

where

=

is the effective stress vector and

−

is the pore pressure.

RI PT

247

(36)

The porochemothermoelasticity is verified for a case where an inclined borehole is drilled

254

in a transversely isotropic formation. The analytical solution assumes that the plane of elastic

255

isotropy (x-y plane) is always perpendicular to the borehole axis (Ekbote and Abousleiman,

256

2005). This assumption clearly alienates most cases encountered in everyday drilling

257

operations, but will be enforced on the numerical model for the purpose of validating against

258

the aforementioned analytical solution. This type of anisotropy will be referred to as

259

“conditional anisotropy” from this point forward. For the purpose of this work, we define an

260

anisotropy ratio as

261

laboratory measurements are required to properly describe the anisotropic medium for real

262

field applications.

263

transformations parameters, far-field stresses, rock and fluid parameters are summarized in

264

Tables 2-4. The input parameters are representative of a low permeability shaley-sand in the

265

Middle East with structured clay that is fully-saturated with formation water, chemically-

266

active and has a membrane-like behavior, where osmotic flow is observed. The thermal

267

parameters

268

porochemothermoelastic model is solved for 0.001 day, 0.01 day and 0.1 day time steps. The

269

results of the pore pressure, solute mass fraction, temperature and total stresses are presented

270

in Figs. 4 and 5. The analytical and numerical solutions are represented by solid lines and

271

markers, respectively. The numerical FEM model demonstrates close agreement with the

272

analytical solution.

273

5. Sensitivity Analysis

|~=

=

•}€=

•}~=

=

‚ •} = ‚ •~ =

to estimate the properties in the MP -direction. Proper

EP

TE D

In this example, the anisotropy ratio is set to 2.0. The coordinate

slightly

modified

from

Ghassemi

and

Diek

(2002).

The

AC C

are

|}=

M AN U

SC

253

274

Several sensitivity analyses are conducted in this section on the porochemothermoelastic

275

model, which are used to emphasize the implications of different aspects of the model. First,

276

the combined effects of elastic and thermic anisotropy are discussed. The stresses are

ACCEPTED MANUSCRIPT 277

calculated for heating and cooling cases using different solute concentrations. The influence

278

of thermal convection and solute advection are accentuated with additional examples. Finally,

279

pore pressure and effective stresses are calculated using different constitutive models, and are

280

juxtaposed to illustrate the individual and overall combined effects of each model.

281

5.1. Effect of Rock Anisotropy The rock anisotropy only affects the elastic and thermic parameters. The chemical

283

parameters are mainly properties of the exchanged fluids, and the anisotropy only affects the

284

coupling chemoelastic coefficients. In the following example, the same input data are used,

285

and the anisotropy ratio is varied to 1.0 (isotropic), 1.5 and 2.0. The effective tangential and

286

axial stresses are calculated using the “conditional” and “complete” anisotropy cases, where

287

the plane of isotropy of the TI material is not necessarily perpendicular to the borehole axis.

288

The radial distributions of effective tangential and axial stresses are summarized in Fig. 6,

289

where the solid and dashed lines represent the conditional and complete anisotropy,

290

respectively. It is shown that the degree of anisotropy has a major influence on the effective

291

stresses around and away from the wellbore region. A noticeable divergence is also observed

292

between the results obtained using the conditional and complete anisotropy cases. Figure 7

293

closely examines the effective axial stress around the borehole wall. The isotropic results are

294

mainly tensile (negative stress). Moreover, the conditional anisotropy assumption

295

consistently overestimates the actual stresses when compared to the complete anisotropy

296

results. It is important to note that these observations should not be generalized, and should

297

be examined on a case-by-case basis. Consequently, assuming isotropy or conditional

298

anisotropy could have detrimental effects on wellbore stability by underestimating borehole

299

collapse pressure, or underestimating the fracture initiation pressure. The effect of anisotropy

300

on pore pressure is not presented in this example, as the effect was trivial.

301

5.2. Effect of Thermal Gradient and Solute Concentration

AC C

EP

TE D

M AN U

SC

RI PT

282

302

The effect of heating and cooling is explored in this subsection by varying the solute

303

mass fraction of the wellbore and the reservoir fluids. This will accentuate individual and

304

combined effects of the porochemothermoelastic model. The temperature differential

306

(∆ =

307

These values are reversed or equalized for the respective cases. The sensitivity analysis for

308

pore pressure, solute concentration, total radial and tangential stresses are presented in Figs.

305

ƒ

−

!)

is set to -50 and 50°C for the cooling and heating scenarios, respectively. In

addition, the solute mass fraction is set to 0.1 and 0.2 for the rock (

!

) and the borehole (

ƒ ).

ACCEPTED MANUSCRIPT 309

8-10. The solid lines represent the heating cases, and the dashed lines represent the cooling

310

ones. The sensitivity analysis for pore pressure, only, is done using two different

312

permeabilities, 1 mD and 1x10-5 mD. In Fig. 8a, the lower permeability makes the pore

313

pressure largely dominated by temperature. Heating the formation will result in thermal

314

expansion, which will subsequently increase the pore pressure. The inverse case, wherein

315

shrinkage occurs, will results in lowering the pore pressure and total stress. The lower

316

permeability allows the fluid exposed to the transient temperature region to be volumetrically

317

enhanced before it diffuses. On the other hand, Fig. 8b, where high permeability is used,

318

illustrates that pore pressure is heavily dominated by the chemical osmosis process. This is

319

evident when the effect of chemical osmotic flow is neglected (

320

becomes mainly a porothermoelastic one. The effect of temperature on pore pressure is

321

minute when high intrinsic permeability is present. A higher solute concentration in the

322

wellbore creates an osmotic backflow that lowers the pore pressure significantly in the area

323

around the borehole.

RI PT

311

=

ƒ)

and the problem

M AN U

SC

!

The total radial and tangential stresses summarized in Figs. 9 and 10. Heating shows a

325

significant increasing effect on the total stresses. On the other hand, a higher solute

326

concentration in the borehole has a lowering effect on the total stresses in the region around

327

the borehole. Injecting a cooler fluid with a higher solute concentration creates a considerable

328

stabilizing effect on the near wellbore region by lowering the pore pressure and stresses. This

329

effect can be exploited in wellbore stability or hydraulic fracturing problems.

330

5.3. Effect of Heat Transfer by Convection and Solute Transfer by Advection

EP

TE D

324

Most conventional subsurface rocks such as sandstone and carbonate reservoirs, exhibit a

332

high intrinsic permeability in the range of 1 mD to tens of Darcys. In which case, the heat

333

transfer by convection and solute transport by advection carried by the fluid plays a

334

significant role and should not be overlooked. In the following example, the effect of heat

335

convection and solute transport by advection are investigated. The process of estimating the

336

convection and advection contributions is a non-linear one as discussed in the Finite

337

Element Formulae section. The results of the conductive-advective model, which will be

338

referred to as “porochemothermoelastic C.A.,” are contrasted against the linear

339

porochemothermoelastic model in Figs. 11 and 12.

AC C

331

340

The pore pressure, solute mass fraction and temperature radial distributions shown in Fig.

341

11 illustrate that the linear porochemothermoelastic results significantly overestimates the

ACCEPTED MANUSCRIPT reservoir pore pressure, solute mass fraction and temperature. At longer time steps, the heat

343

transfer by convection and solute transport by advection eclipse the contributions by

344

conduction and diffusion as shown by the large separation between the two models. The

345

added solute contributed by advection increases the counteracting chemical osmotic

346

backflow. Subsequently, the incremental change in pore pressure translates into the effective

347

stresses.

RI PT

342

The total radial, tangential and axial stresses are summarized in Fig. 12. The linear

349

porochemothermoelastic model, shown in solid lines, consistently overestimates the stresses.

350

This is attributed to the underestimation of the cooling effect carried by fluid convection,

351

which consequently decreases the thermal stresses. In addition, the overestimation of solute

352

mass increases the effect of swelling, and consequently increases the total stresses. The effect

353

of advection and convection can be neglected if the borehole wall is impermeable, where the

354

borehole BC only affects the traction force and heat transfer by conduction. Advection and

355

convection problems are computationally intensive due to the non-linearity of the solution.

356

Kanfar et al. (2016a, 2016b) concluded that advection and convection contributions can be

357

neglected when the rock permeability is lower than 1x10−4 mD and 1x10−1 mD, respectively.

358

5.4. Effect of Different Constitutive Models

M AN U

SC

348

In this section, a comparison between the different constitutive models is investigated.

360

The poroelastic, porothermoelastic, porochemoelastic, porochemothermoelastic and the

361

convective-advective (C.A.) porochemothermoelastic models are used to calculate the pore

362

pressure, effective radial, tangential and axial stresses presented in Figs. 13-16. The same set

363

of inputs is used with the exception of rock permeability, which is changed to 1x10-4 mD to

364

magnify the effect of each model and avoid permeability bias as previously discussed.

EP

TE D

359

It was previously demonstrated that the thermal effect on pore pressure diminishes in

366

high permeability rocks. This can be attributed to the short time of exposure of the pore fluid

367

for it to thermally expand or shrink before it mobilizes outside of the radius of the thermally

368

enhanced area. In contrast, Fig. 13 presents the results for pore pressure away from the

369

borehole in a low permeability rock. Caused by cooling, the pore pressure estimated by the

370

porothermoelastic model shows significant divergence from the isothermal poroelastic model.

371

Similarly, the pore pressure obtained by the porochemothermoelastic model diverges from

372

the porochemoelastic model caused by the temperature enhancement which lowers the pore

373

pressure and the solute concentration.

AC C

365

ACCEPTED MANUSCRIPT The enhancement of the different models on the effective radial stresses in Fig. 14 is less

375

pronounced when compared to the tangential and axial stresses in Figs. 15 and 16,

376

respectively. It is clear that there is virtually no added effect by the conductive-advective

377

model when compared to the linear porochemothermoelastic model. The ultra-low

378

permeability makes solute concentration diffusion dominated, and the heat transfer mainly a

379

conduction problem. Consequently, the convection and advection contribution can be

380

neglected at lower permeabilities.

381

6. Wellbore Stability

RI PT

374

In this section, wellbore stability is analyzed using different constitutive models. The

383

chemical and thermal effects are mainly discussed. The effect of elastic anisotropy and rock

384

strength anisotropy on time-dependent wellbore stability was discussed previously by the

385

authors (Kanfar et al., 2015a). An improved mud window estimation for inclined boreholes

386

drilled in highly-stressed anisotropic formations to prevent shear and tensile failures using a

387

semi-analytical model was developed by the authors (Kanfar et al., 2015b).

M AN U

SC

382

Figure 17 summarizes the pore pressure, effective radial and tangential stress

389

distributions around the borehole using different constitutive models after 0.1 day from

390

drilling. The input data used in this section are summarized in Tables 1-4. Borehole shear

391

failure is assessed using the Mohr-Coulomb criterion (Jaeger et al., 2009). A linearly elastic-

392

brittle model is assumed, where shear failure occurs when the compressive strength of the

393

rock is exceeded. It is important to note that the failure analysis does not use a remeshing

394

technique or post-failure analysis between time steps to account for the failed elements. The

395

analysis provides a practical tool to prognosticate whether the borehole will stabilize or

396

worsen as time progresses. The strength parameters used for the Mohr-Coulomb criterion are

397

4.14 MPa and 20 degrees for the cohesive strength and angle of internal friction, respectively.

398

The corresponding borehole failure is presented in Fig. 18, where the failure zone is shown in

399

blue. If we consider the results of the poroelastic model (Fig. 18a) as a base case, the results

400

of the porothermoelastic model (Fig. 18b) and the porochemoelastic model (Fig. 18c) show a

401

more stable wellbore. This stabilizing effect is caused by temperature cooling and the higher

402

solute mass injected, respectively. The combined effects of the porochemothermoelastic

403

model (Fig. 18d) show a significant added improvement to the stability of the wellbore when

404

juxtaposed against the poroelastic model. As shown in Fig. 17d, the effective tangential stress

405

is considerably lower for the porochemothermoelastic model, which creates a stable

AC C

EP

TE D

388

ACCEPTED MANUSCRIPT 406

environment in the near wellbore region. Similar concepts can be used during hydraulic

407

fracture operations to create a more favorable environment to induce a tensile failure.

408

7. Conclusions In this work, constitutive and transport equations for fully coupled chemo-thermo-hydro-

410

mechanical processes are developed for homogenous anisotropic formations. A numerical

411

model using the FEM is developed to simulate the drilling of an inclined borehole in

412

transversely isotopic shaley sand. It was assumed that the formation pores are fully-saturated

413

with water, and the rock has a membrane-like characteristics, where chemical osmosis is

414

applicable. A novel pseudo 3D analysis using the generalized plane strain concept is used for

415

more efficient and stable computing without compromising the accuracy of the results.

417

•

418 419

SC

Several conclusions are drawn from the work presented in this paper as follows:

M AN U

416

RI PT

409

Anisotropy should be properly characterized by acquiring core data in the formation of interest to assess the degree of elastic anisotropy.

•

The presented examples demonstrated that assuming isotropy or conditional

420

anisotropy could potentially lead to large errors in the estimated stresses if high

421

degree of elastic anisotropy is present, especially in the region around borehole. •

It is assumed that the rock is fully saturated with formation brine, but the model

TE D

422 423

can be generalized for live oil or gas-filled formations for qualitative estimates of

424

optimum borehole pressure selection to avoid borehole collapse. •

426

428 429 430 431

model inputs are representative of the rock in question. •

Coupling heat flow is especially relevant in deep wells, where there is a significant temperature differential between the drilling fluid and the original

AC C

427

The model is inclusive of all rock types (clastics or carbonates) as long as the

EP

425

geothermal gradient of the formation.

•

The presented examples demonstrated that heating the formation increases the

pore pressure and the total stresses, and the opposite is observed when the

432

formation is cooled. Assuming an isothermal model will result in significant

433

errors in the estimated stresses.

434

•

Coupling solute transport, chemical potential and formation swelling should be

435

carefully investigated as the model is not always applicable. The model should be

436

used when the formations exhibits a membrane like behavior where the above

437

mentioned elements are relevant.

ACCEPTED MANUSCRIPT 438

•

It is observed from the presented examples, a higher solute concentration in the borehole fluid creates a higher chemical potential in the reservoir fluid. This

440

results in a reverse osmosis and subsequently a drop in pore pressure.

441

Concurrently, the effective radial stress increased, and the effective tangential

442

stress decreased. This has a stabilizing effect on the borehole wall, which is in

443

line with the drilling industry common practices.

444

•

RI PT

439

Advection and convection are computationally intensive and should only be used when the permeability thresholds are reached (1x10−4 mD and 1x10−1 mD,

446

respectively).Injecting a cooler fluid with a higher salinity than the intrinsic

447

formation fluid adds a significant stabilizing effect on the near wellbore region by

448

lowering the effective tangential stress.

SC

445

The combined effects of the thermo-chemo-hydro-mechanical model produce different

450

results from the individual thermoelastic or chemoelastic models. The different elements of

451

the anisotropic porochemothermoelastic model are interwoven with one another, and each

452

component should not be considered in isolation. In this work, the model was applied for an

453

inclined borehole problem, but can be generalized to a wide range of applications.

454

References

455

Aadnoy, B., Chenevert, M., 1987. Stability of highly inclined boreholes. SPE Drill. Eng.

456 457

Abousleiman, Y., Ekbote, S., 2005. Solutions for the Inclined Borehole in a Porothermoelastic Transversely Isotropic Medium. J. Appl. Mech. 72, 102–114.

458 459 460

Abousleiman, Y., Roegiers, J., Cui, L., Cheng, A., 1995. Poroelastic solution of an inclined borehole in a transversely isotropic medium, in: Rock Mechanics: Proceedings of the 35th US Symposium on Rock Mechanics. Taylor & Francis Group.

461 462

Amadei, B., 1983. Rock anisotropy and the theory of stress measurements. Springer-Verlag, Berlin; New York.

463 464

Bader, S., Kooi, H., 2005. Modelling of solute and water transport in semi-permeable clay membranes: comparison with experiments. Adv. Water Resour. 28, 203–214.

465 466

Batugin, S.A., Nirenburg, R.K., 1972. Approximate relation between the elastic constants of anisotropic rocks and the anisotropy parameters. Sov. Min. Sci. Sov. Min. Sci. 8, 5–9.

467 468

Bear, J., Corapcioglu, M.Y., 1981. A mathematical model for consolidation in a thermoelastic aquifer due to hot water injection or pumping. Water Resour. Res. 17, 723–736.

469 470

Biot, M., 1955. Theory of Elasticity and Consolidation for a Porous Anisotropic Solid. J. Appl. Phys. 26, 182–185.

471 472

Cheng, A., 1997. Material coefficients of anisotropic poroelasticity. Int. J. Rock Mech. Min. Sci. 34, 199–205.

AC C

EP

TE D

M AN U

449

ACCEPTED MANUSCRIPT Coussy, O., 1989. A general theory of thermoporoelastoplasticity for saturated porous materials. Transp. Porous Media 4, 281–293.

475 476 477

Cui, L., Ekbote, S., Abousleiman, Y., Zaman, M., Roegiers, J., 1998. Borehole stability analyses in fluid saturated formations with impermeable walls. Int. J. Rock Mech. Min. Sci. 35, 582–583.

478 479

Detournay, E., Cheng, A., 1988. Poroelastic response of a borehole in a non-hydrostatic stress field. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 25, 171–182.

480 481

Ekbote, S., Abousleiman, Y., 2005. Porochemothermoelastic Solution for an Inclined Borehole in a Transversely Isotropic Formation. J. Eng. Mech. 131, 522–533.

482 483

Ekbote, S., Abousleiman, Y., 2006. Porochemoelastic Solution for an Inclined Borehole in a Transversely Isotropic Formation. J. Eng. Mech. 132, 754–763.

484

Fjaer, E., 2008. Petroleum related rock mechanics. Elsevier, Amsterdam; Boston.

485 486

Ghassemi, A., Diek, A., 2002. Porothermoelasticity for swelling shales. J. Pet. Sci. Eng. 34, 123–135.

487 488

Ghassemi, A., Diek, A., 2003. Linear chemo-poroelasticity for swelling shales: theory and application. J. Pet. Sci. Eng. 38, 199–212.

489 490 491

Ghassemi, A., Tao, Q., Diek, A., 2009. Influence of coupled chemo-poro-thermoelastic processes on pore pressure and stress distributions around a wellbore in swelling shale. J. Pet. Sci. Eng. 67, 57–64.

492 493

Ghassemi, A., Wolfe, G., 1999. A chemo-mechanical model for borehole stability analyses, in: USRMS. American Rock Mechanics Association, pp. 239–246.

494 495

Groot, S.R. de, 1951. Thermodynamics of Irreversible Processes. North-holland Pub. Co.; [sole distributors for U.S.A.: Interscience Publishers, New York], Amsterdam.

496 497

Jaeger, J., Cook, N., Zimmerman, R., 2009. Fundamentals of rock mechanics. Chapman and Hall ; Wiley, London; New York.

498 499

Kaliakin, V.N., 2002. Introduction to approximate solution techniques, numerical modeling, and finite element methods. Marcel Dekker, New York.

500 501

Kanfar, M.F., Chen, Z., Rahman, S., 2016a. Anisotropic Diffusive-Advective Porochemoelasticity Modeling for Inclined Boreholes. Int. J. Geomech. 6016025.

502 503

Kanfar, M.F., Chen, Z., Rahman, S.S., 2015a. Effect of material anisotropy on timedependent wellbore stability. Int. J. Rock Mech. Min. Sci. 78, 36–45.

504 505

Kanfar, M.F., Chen, Z., Rahman, S.S., 2015b. Risk-controlled wellbore stability analysis in anisotropic Formations. J. Pet. Sci. Eng. 134, 214–222.

506 507 508

Kanfar, M.F., Chen, Z., Rahman, S.S., 2016b. Fully coupled 3D anisotropic conductiveconvective porothermoelasticity modeling for inclined boreholes. Geothermics 61, 135– 148.

509 510

Kurashige, M., 1989. A thermoelastic theory of fluid-filled porous materials. Int. J. Solids Struct. 25, 1039–1052.

511

Lekhnitskii, S.G., 1981. Theory of elasticity of an anisotropic body. Mir Publishers, Moscow.

512

McLean, M.R., Addis, M.A., 1990. Wellbore Stability: The Effect of Strength Criteria on

AC C

EP

TE D

M AN U

SC

RI PT

473 474

ACCEPTED MANUSCRIPT Mud Weight Recommendations, in: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.

515 516

Mody, F.K., Hale, A., 1993. Borehole-stability model to couple the mechanics and chemistry of drilling-fluid/shale interactions. J. Pet. Technol. 45.

517 518

Nguyen, V.X., Abousleiman, Y.N., 2010. Incorporating electrokinetic effects in the porochemoelastic inclined wellbore formulation and solution. An. Acad. Bras. Cienc.

519 520

Nur, A., Byerlee, J.D., 1971. An exact effective stress law for elastic deformation of rock with fluids. J. Geophys. Res. 76, 6414–6419.

521

Ong, S.H., 1994. Borehole Stability, Ph.D. Thesis. The University of Oklahoma.

522 523

Schiffman, R., 1971. A thermoelastic theory of consolidation. Environ. Geophys. Heat Transf. 78–84.

524 525

Schrefler, B., 1985. A partitioned solution procedure for geothermal reservoir analysis. Commun. Appl. Numer. methods 1, 53–56.

526 527

Sherwood, J.D., 1993. Biot Poroelasticity of a Chemically Active Shale. procmathphysscie Proc. Math. Phys. Sci. 440, 365–377.

528 529

Tran, M.H., Abousleiman, Y., 2013. Anisotropic Porochemoelectroelastic Solution for an Inclined Wellbore Drilled in Shale. J. Appl. Mech. 80, 20912.

530 531 532

Vishal, V., Jain, N., Singh, T.N., 2015. Three dimensional modelling of propagation of hydraulic fractures in shale at different injection pressures. Sustain. Environ. Res. 25, 217–225.

533 534

Vishal, V., Ranjith, P.G., Singh, T.N., 2013. CO2 permeability of Indian bituminous coals: Implications for carbon sequestration. Int. J. Coal Geol. 105, 36–47.

535 536

Wang, Y., Dusseault, M.B., 2003. A coupled conductive convective thermo poroelastic solution and implications for wellbore stability. J. Pet. Sci. Eng. 38, 187–198.

537 538 539

Yasuhara, H., Kinoshita, N., Ogata, S., Cheon, D.-S., Kishida, K., 2016. Coupled thermohydro-mechanical-chemical modeling by incorporating pressure solution for estimating the evolution of rock permeability. Int. J. Rock Mech. Min. Sci. 86, 104–114.

540 541

Yin, S., Dusseault, M.B., Rothenburg, L., 2011. Coupled THMC modeling of CO2 injection by finite element methods. J. Pet. Sci. Eng. 80, 53–60.

542 543 544

Zhang, D., Ranjith, P.G., Perera, M.S.A., 2016. The brittleness indices used in rock mechanics and their application in shale hydraulic fracturing: A review. J. Pet. Sci. Eng. 143, 158–170.

545 546

Zhou, X., Ghassemi, A., 2009. Finite element analysis of coupled effects around a wellbore in swelling shale. Int. J. Rock Mech. Min. Sci. 46, 769–778.

547 548

AC C

EP

TE D

M AN U

SC

RI PT

513 514

ACCEPTED MANUSCRIPT 549

Appendix

550

The finite element matrices are defined as follows: oh = „ … . … †Ω ‡

oj = „ … . ‡

(A.1)

gj †Ω

(A.2)

RI PT

ok = „ … . gk l †Ω

(A.3)

‡

o. = „ … . Λ g. †Ω ℙj = „ gj. ‡

1

(A.4)

+

gj †Ω

SC

‡

‡

.

M AN U

ℙjq = „ H∇gj J X H∇gj J †Ω ℙk = „ gj. Ξ gk l †Ω ‡

ℙkq = „ H∇gj J V H∇gk l J †Ω ℙ. = „ gj. Ω g. †Ω ‡

TE D

‡

.

‡

‡

ℂkq = „ H∇gk l J 11 − ℜ2

AC C

.

‡

ℂ.q = „ H∇gk l J 11 − ℜ2 ‡

s . = „ g.. ‡

551

.

Q Q g.

†Ω

1∇g. 2 †Ω

s .q = „ 1∇g. 2. ^Q 1∇g. 2 †Ω − 1∇g. 2. !̅ ‡

(A.8)

(A.11)

H∇gk l J †Ω − „ H∇gk l J 11 − ℜ2 U ! Hgk l J †Ω .

(A.7)

(A.10)

EP

ℂk = „ gk. l 5c gk l †Ω

(A.6)

(A.9)

ℙ.q = „ H∇gj J Y . 1∇g. 2 †Ω .

(A.5)

.

‡

(A.13) (A.13) (A.14)

! U

!

1g. 2 †Ω

(A.15)

ACCEPTED MANUSCRIPT Fig. 1 Coordinate transformation systems (Kanfar et al., 2016). Fig. 2 The green shapes represent the well trajectory, and the red shapes represent the FEM mesh. The formation beddings are shown in blue (Kanfar et al., 2016). Fig. 3 Eight-node quadrilateral FEM mesh (Kanfar et al., 2016).

RI PT

Fig. 4 Comparison between numerical and analytical results for pore pressure, solute mass fraction and temperature at different time steps. Fig. 5 Comparison between numerical and analytical results for total radial, tangential and axial stresses at different time steps.

SC

Fig. 6 Effect of elastic and thermal anisotropy on effective tangential and axial stresses with radial distance (t = 0.01 day).

M AN U

Fig. 7 Effect of elastic and thermal anisotropy on effective axial stress around the borehole wall (t = 0.01 day). Fig. 8 Effect of reservoir heating and cooling on pore pressure using different solute gradients (t = 0.01 day). Fig. 9 Effect of reservoir heating and cooling on total radial stresses using different solute gradients (t = 0.01 day).

TE D

Fig. 10 Effect of reservoir heating and cooling on total tangential stresses using different solute gradients (t = 0.01 day). Fig. 11 Comparison between conductive and conductive-convective heat transfer models on pore pressure and temperature.

EP

Fig. 12 Comparison between conductive and conductive-convective models on effective radial, tangential and axial stresses. Fig. 13 Pore pressure calculation using different models (t = 0.01 day).

AC C

Fig. 14 Effective radial stress calculation using different models (t = 0.01 day). Fig. 15 Effective tangential calculation using different models (t = 0.01 day). Fig. 16 Effective axial calculation using different models (t = 0.01 day).

Fig. 17 Pore pressure, effective radial and tangential stress (in MPa units) distribution using a) poroelasticity, b) porothermoelasticity, c) porochemoelasticity and d) porochemothermoelasticity (t = 0. 1 day). Fig. 18 Wellbore shear failure after 0.1 day using a) poroelasticity, b) porothermoelasticity, c) porochemoelasticity and d) porochemothermoelasticity (t = 0. 1 day). Failure zone is shown in blue.

ACCEPTED MANUSCRIPT

10 15 25 20 30

MPa MPa MPa MPa MPa

AC C

EP

TE D

Table 4 Rock and fluid properties. Porosity, Permeability, Fluid viscosity, Matrix bulk modulus, Fluid bulk modulus, Initial reservoir temperature, Wellbore temperature, Matrix thermal conductivity, Fluid thermal conductivity, Matrix density, Fluid density, ̅ Solid matrix specific heat, Fluid specific heat, Matrix thermal expansion coefficient, Fluid thermal expansion coefficient, Young's modulus, , Poisson's ratio, Specific entropy, Reflection coefficient, ℜ Retardation coefficient, " Thermal filtration coefficient, # $ Formation solute mass fraction,

Degree Degree Degree Degree Degree Degree

SC

Table 3 Inputs for far-field stresses. Pore pressure, Well pressure, Maximum horizontal stress, Minimum horizontal stress, Vertical (overburden) stress,

30 60 50 30 70 7

M AN U

Table 2 Inputs for stress and property rotations. Well azimuth, Well inclination, Rock dip azimuth, Rock dipping, Maximum horizontal stress azimuth, Overburden stress inclination,

RI PT

Table 1 FEM program boundary conditions. Type Well Nodes Outer Boundary Nodes Boundary Type Traction force Well pressure Far-field stresses Neumann Pore pressure Well pressure Initial pore pressure Dirichlet Temperature Well temperature Initial temperature Dirichlet = Solute mass fraction Dirichlet =

,

0.10 1 3 x 10-4 70 2.30 394 344 1.30 0.580 2640 1113 768 4181 1.8 x 10-5 3 x 10-4 24.14 0.30 3686 0.07 1 6 x 10-12 0.1

mD Pa ∙ s GPa GPa K K W/m/K W/m/K Kg/m3 Kg/m3 J/Kg/K J/Kg/K 1/K 1/K GPa J/Kg/K m2/s/K -

ACCEPTED MANUSCRIPT 0.2 5 x 10-9 0.0585 1.50

m2/s Kg/mol MPa

AC C

EP

TE D

M AN U

SC

RI PT

Wellbore solute mass fraction, Solute diffusion coefficient, # Solute molar mass, % Osmotic (swelling) coefficient, &

ACCEPTED MANUSCRIPT

Zs

Zb

Zr

V γs

Ys

V γr

Yr

γb

Yb Xb

Xr

TE D

M AN U

SC

Xs

RI PT

V

AC C

EP

E ψs N

Far-field stresses

N

E

E

ψr

Rock Properties

ψb N

Borehole

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

𝝈𝒁𝒔

Zb

AC C

𝝈𝑿𝒔

EP

TE D

Yb

Θ

BOH

Xb

𝝈𝒀𝒔

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

10

5

RI PT

Pore Pressure [MPa]

15

150

SC Analytical solution Numerical solution t = 0.001 Day t = 0.01 Day t = 0.1 Day

100

50 0.1

1

EP

0.05 0.1

TE D

0.15

200

Temperature [C]

M AN U

0.2

0.1

1

AC C

Solute Mass Fraction [-]

0 0.1

Radial Distance [m]

(at angle = 90 degrees)

1

ACCEPTED MANUSCRIPT

22 20

RI PT

18 16

M AN U

40 35

30 25

Analytical solution Numerical solution t = 0.001 Day t = 0.01 Day t = 0.1 Day

20 15 10 0.1

1

EP

20 0.1

TE D

30 25

1

SC

14 0.1

AC C

Total Radial Stress [MPa] Total Tangential Stress [MPa] Total Axial Stress [MPa]

24

Radial Distance [m]

(at angle = 90 degrees)

1

ACCEPTED MANUSCRIPT

RI PT

25

M AN U

SC

20

5

0

EP

10

TE D

15

AC C

Effective Tangential and Axial Stresses [MPa]

30

Conditional Anisotropy Complete Anisotropy Isotropic

EXr / EZr EXr / EZr

vXYr / vXZr vXYr / vXZr

-5 0.1

1.5 2.0 1

Radial Distance [m]

(at angle = 90 degrees)

ACCEPTED MANUSCRIPT

4

RI PT SC

0

M AN U

-2

TE D

-4

EP

-6 -8 -10 -12

Conditional Anisotropy Complete Anisotropy Isotropic

AC C

Effective Axial Stress [MPa]

2

0

60

120

EXr / EZr EXr / EZr 180

240

vXYr / vXZr vXYr / vXZr 300

Angle Around the Borehole Wall [Degree]

1.5 2.0 360

18

18

ACCEPTED MANUSCRIPT

14

16

14

12

RI PT

12

10

8

M AN U

6

10

8

4 0.1

Radial Distance [m]

(at angle = 90 degrees)

EP

12

4

16

14

AC C

14

1

18

TE D

18

16

6

Permeability = 1x10-5 mD

(a)

Pore Pressure [MPa]

50 C ‐50 C

SC

Pore Pressure [MPa]

16

Heating (∆T Cooling (∆T CSw > CSf CSw = CSf CSw < CSf

12

10

10

8

8

6

6

Permeability = 1 mD

(b) 4 0.1

1

Radial Distance [m]

10

(at angle = 90 degrees)

4 100

ACCEPTED MANUSCRIPT

25

RI PT

25

M AN U

SC

23

21

TE D

21

19

EP

19

17

15 0.1

Heating (∆T Cooling (∆T CSw > CSf CSw = CSf CSw < CSf

AC C

Total Radial Stress [MPa]

23

1

Radial Distance [m]

10

(at angle = 90 degrees)

50 C ‐50 C

17

15 100

ACCEPTED MANUSCRIPT

50 C ‐50 C 54

SC

RI PT

54

47

TE D

M AN U

47

40

EP

40

61

Heating (∆T Cooling (∆T CSw > CSf CSw = CSf CSw < CSf

33

26 0.1

AC C

Total Tangential Stress [MPa]

61

Radial Distance [m]

33

(at angle = 90 degrees)

1

26

ACCEPTED MANUSCRIPT

16 14 12 10

SC

0.2

150

1

10

EP

0.05 0.1

TE D

0.15 0.1

10

M AN U

0.25

200

Temperature [C]

1

AC C

Solute Mass Fraction [-]

8 0.1

RI PT

Pore Pressure [MPa]

18

Porochemothermoelastic Porochemothermoelastic (C.A.) t = 0.001 Day t = 0.01 Day t = 0.1 Day

100

50 0.1

1

Radial Distance [m]

(at angle = 90 degrees)

10

ACCEPTED MANUSCRIPT

20

10 0.1

M AN U

50 40

40 30

EP

10 0.1

TE D

30 20

10

SC

1

RI PT

15

1

10

AC C

Total Radial Stress [MPa] Total Tangential Stress [MPa] Total Axial Stress [MPa]

25

Porochemothermoelastic Porochemothermoelastic (C.A.) t = 0.001 Day t = 0.01 Day t = 0.1 Day

20 10 0 0.1

1

Radial Distance [m]

(at angle = 90 degrees)

10

ACCEPTED MANUSCRIPT

15

M AN U

13

10

9 0.1

EP

TE D

12

11

RI PT SC

+ x

AC C

Pore Pressure [MPa]

14

Poroelastic Porothermoelastic Porochemoelastic Porochemothermoelastic Porochemothermoelastic (C.A.)

Radial Distance [m]

(at angle = 90 degrees)

1

ACCEPTED MANUSCRIPT

16

RI PT M AN U

SC

12

6

4

2 0.1

EP

8

TE D

10

AC C

Effective Radial Stress [MPa]

14

+ x Radial Distance [m]

Poroelastic Porothermoelastic Porochemoelastic Porochemothermoelastic Porochemothermoelastic (C.A.)

(at angle = 90 degrees)

1

ACCEPTED MANUSCRIPT

38

SC M AN U

32

24

EP

26

TE D

30 28

RI PT

+ x

34

Poroelastic Porothermoelastic Porochemoelastic Porochemothermoelastic Porochemothermoelastic (C.A.)

AC C

Effective Tangential Stress [MPa]

36

22 20 0.1

Radial Distance [m]

(at angle = 90 degrees)

1

ACCEPTED MANUSCRIPT

18

RI PT

17

M AN U

SC

15 14

11 10

EP

12

TE D

13

AC C

Effective Axial Stress [MPa]

16

+ x

9 8 0.1

Radial Distance [m]

Poroelastic Porothermoelastic Porochemoelastic Porochemothermoelastic Porochemothermoelastic (C.A.)

(at angle = 90 degrees)

1

Effective Tangential Stress

SC

RI PT

(a)

EffectiveMANUSCRIPT Radial Stress ACCEPTED

0.40 m

Pore Pressure

TE D EP AC C

(c)

M AN U

(b)

(d)

RI PT

(a)

EP

TE D

M AN U

SC

(b)

AC C

0.30 m

ACCEPTED MANUSCRIPT

(c)

(d)

ACCEPTED MANUSCRIPT •

New numerical model for simulating anisotropic formations under THMC loadings.

•

The model addresses chemically-active, anisotropic non-isothermal formations.

•

The model can be used to assess time dependent wellbore stability.

•

The model addresses inclined boreholes drilled in arbitrarily oriented rock

EP

TE D

M AN U

SC

A novel pseudo-3D analysis is used to models 3D problems using 2D FEM mesh.

AC C

•

RI PT

laminations.

ACCEPTED MANUSCRIPT

Λ

Unit Pa Pa 1/Pa Pa Pa Pa/K K Pa 1/K mD Pa ∙ s Pa Pa W/m/K W/m/K Kg/m3 Kg/m3 J/Kg/K J/Kg/K 1/K 1/K Pa Pa J/Kg/K m2/s/K m2/s Kg/mol Pa Degree Degree Degree Degree Degree Degree Degree Degree Degree Pa Pa Pa Pa Pa

! "# $# "% $% " $ "# $# "% &' &( )* +* ,*

TE D

EP

ℜ

AC C

̅

M AN U

SC

Ξ Ω

Definition Stress tensor Drained stress-strain stiffness tensor Drained stress-strain Compliance tensor Strain tensor Pore pressure Chemo-mechanical coupling term Solute mass fraction Thermo-mechanical coupling term Temperature Biot’s effective stress coefficient Variation of fluid content per unit referential volume Biot’s modulus Hydro-chemical coupling term Hydro-thermal coupling term Porosity Permeability Fluid viscosity Matrix bulk modulus Fluid bulk modulus Matrix thermal conductivity Fluid thermal conductivity Matrix density Fluid density Solid matrix specific heat Fluid specific heat Matrix thermal expansion coefficient Fluid thermal expansion coefficient Young's modulus Poisson's ratio Shear modulus Specific entropy Reflection coefficient Retardation coefficient Thermal filtration coefficient Solute diffusion coefficient Solute molar mass Osmotic (swelling) coefficient Well azimuth Well inclination Rock dip azimuth Rock dipping Maximum horizontal stress azimuth Overburden stress inclination Well azimuth Well inclination Rock dip azimuth Pore pressure Well pressure Maximum horizontal stress Minimum horizontal stress Vertical (overburden) stress

RI PT

Symbol

Non-Business Use

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT