2.3.1. Effective Stress and Poroelasticity

The effective stress for incompressible rock and the concept of one-dimensional consolidation are the fundamentals of poroelasticity formulated by Terzaghi in 1923 [17]. Successively, using the basic principles of continuum mechanics and applying the concept of the coupling of stress and pore pressure in a porous medium, Biot developed a comprehensive three-dimensional theory of consolidation [17,20]. Biot's theory and the papers that he published are more aligned towards geomechanics than flow models, due to which they are rarely compatible with the coupling of geomechanics and flow models. By introducing the so-called Skempton pore pressure parameters (A and B), Skempton in 1954 procured a relationship between the total stress and the pore pressure under undrained initial loading [17,21]. Later, however, the relationships among pore pressure, stress, and volume and the concept of compressibility in a porous medium were better clarified by Geerstma in 1957 [17,22]. Later, Van der Knaap (1959) extended Geerstma's work to nonlinear elastic geomaterials only for dense and uncemented sands [17,23]. By applying Biot's theory, Geerstma in 1966 examined subsidence problems in oil fields and published prototype geomechanical modelling, which is probably the first-ever coupled analysis of fluid flow [17,20]. Nur and Byerlee (1971) demonstrated that the effective stress law proposed by Biot is far more general and precise than that proposed by Terzaghi [17,24]. Nevertheless, there are certain limitations (one-dimensional analysis, neglect of the compressibility of

fluids and rocks, etc.) that Terzaghi recognized in the assumptions that he made in the 1920s to solve problems of applied rock mechanics in clay consolidation [17]. Later in the 1970s, there were further developments on coupled flow stress issues; e.g., fluid compressibility was introduced into the classical soil mechanical consolidation theory of Ghaboussi and Wilson [17,25]. Rice and Cleary (1976) showed how poroelasticity problems could be solved using pore pressure and stress as primary variables, instead of the displacements used by Biot [17,26].

#### 2.3.2. Simulation Concept and Governing Equations

The numerical modelling of an underground gas reservoir can contribute to the understanding of the interaction mechanisms between the injected gas and the deformation of the reservoir. The injection of cold foreign gas into the reservoir leads to thermal and mechanical disequilibrium in the reservoir by altering the transport properties, including porosity and permeability. A three-dimensional THM coupling model of a reservoir is created by incorporating the mechanical equilibrium equation, the fluid flow or seepage equation (Darcy equation), the heat transfer equation (Fourier equation) of the formation's rock matrix, and the THM stress equation. These equations are based on the porosity, permeability, thermal diffusivity, and other physical and mechanical parameters of the sandstone formation [27].

The mechanical equilibrium equation can be expressed as follows [27]:

$$S\_{i,j,j} + f\_i = 0,\tag{3}$$

where, S*i*,*j*,*<sup>j</sup>* is the total stress tensor (N/m2), and *fi* is the body force (N/m3). The equation of the continuity of the fluid flow in the rock can be written as [27]:

$$\frac{\partial \rho\_l}{\partial t} + \frac{\partial (\rho\_l r v\_r)}{r \partial r} + \frac{\partial (\rho\_l v\_\theta)}{r \partial \theta} + \frac{\partial (\rho\_l w)}{\partial z} = 0,\tag{4}$$

In this equation *vr*, *vθ*, and *w* are the Darcy velocities (m/s) along the radial, hoop, and well-depth directions, respectively, in the porous reservoir. The relationship between stress and porosity/permeability changes in porous rock can be described with the following equations [27,28],

$$
\phi = \phi\_r + (\phi\_0 - \phi\_r) \exp(eS\_M)\_r \tag{5}
$$

$$k = k\_0 \exp(c.(\frac{\Phi}{\Phi\_0} - 1)),\tag{6}$$

In the above equations, *SM* denotes mean effective stress; *φ*<sup>0</sup> and *k*<sup>0</sup> are the porosity and permeability at zero stress, respectively; *φ<sup>r</sup>* represents the residual porosity at high stress; and the exponents *e* and *c* are determined experimentally.

The heat transfer process and the total energy conservation can be express by rewriting the Fourier equation [29]:

$$(\rho c)\_t \frac{\partial T}{\partial t} - \frac{1}{r} \frac{\partial}{\partial r} (k\_l r \frac{\partial T}{\partial t}) - \frac{1}{r^2} \frac{\partial}{\partial \theta} (k\_t \frac{\partial T}{\partial \theta}) - \frac{\partial}{\partial z} (k\_t \frac{\partial T}{\partial z}) - q\_l r = 0,\tag{7}$$

$$(\rho c)\_t = (1 - \phi)c\_s\rho\_s + \phi c\_f\rho\_{f'} \tag{8}$$

$$k\_{total} = \phi k\_f + (1 - \phi)k\_{s\prime} \tag{9}$$

In the above equations, the total heat capacity of the solid and fluid phases is denoted by (*pc*)*<sup>t</sup>* in units (J/(m3. ◦C)); *ktotal* is the total thermal conductivity (J/(m.s.◦C)); *qt* is the intensity of the internal heat source (J/(m3.s)); *cs* and *cf* are the specific heat capacities of the formation and the fluid, respectively (J/kg.K); *ρ<sup>s</sup>* and *ρ<sup>f</sup>* are the density of the formation and the fluid (Kg/m3), respectively; and *ks* and *k <sup>f</sup>* are the thermal conductivity of the formation and the fluid (W/(m·K)), respectively [29].

Finally, the governing equation involving all important thermal-hydro-mechanical parameters is as follows [30]:

$$2a\frac{(1-2\nu)}{1+\nu}\nabla^2 p + 6\beta \mathcal{K}\_B \frac{(1-2\nu)}{1+\nu} \nabla^2 T - \nabla \cdot f - 3\frac{(1-\nu)}{1+\nu} \nabla^2 \mathcal{S}\_{m(total)} = 0,\tag{10}$$

where *Sm*(*total*) is the mean total stress (MPa), *ν* is the Poisson's ratio of the rock mass, *α* is the Biot coefficient, *β* is the coefficient of linear thermal expansion (1/◦C), and *KB* is the bulk modulus of the rock (MPa). The term (2*α* (1−2*ν*) <sup>1</sup>+*<sup>ν</sup>* ∇<sup>2</sup> *<sup>p</sup>*) describes the effect of poroelastic stress; (6*βKB* (1−2*ν*) <sup>1</sup>+*<sup>ν</sup>* <sup>∇</sup>2*T*) represents thermalelastic stress; and the term (∇. *<sup>f</sup>* <sup>−</sup> <sup>3</sup> (1−*ν*) <sup>1</sup>+*<sup>ν</sup>* ∇<sup>2</sup>*Sm*(*total*)) shows the body force [30].
