An uncoupled radial thermoporoelastic model with local thermal non-equilibrium

A new analytic model is introduced to describe rock deformation produced by fluid injection / extraction in geothermal reservoirs. This model is fully coupled under isothermal conditions and thermally uncoupled if local thermal non-equilibrium (LTNE) is considered. The model describes both fluid flow with conductive-convective heat transfer and rock deformation in cylindrical coordinates. There are 13 unknowns in the model: fluid pressure, fluid content variation, solid skeleton displacement, radial, tangential and volumetric strains, two stresses, porosity, solid deformation velocity, fluid velocity, rock and fluid temperatures. Considering LTNE, there is an effective heat transfer between the solid skeleton and the liquid. The porosity is estimated as a function of fluid pressure and temperature. The solid radial deformation is an irrotational vector field, thus the fluid content variation, is proportional to pore pressure, which is calculated using the classical Theis solution. Therefore, the solid radial displacement can be obtained in analytical form. Once the fluid velocity is computed, the fluid temperature is deduced from a new analytical solution of the diffusion-convection equation. This model is didactic, useful and simple. It allows to explore different options for both the fluid and the poroelastic parameters, with different boundary conditions. Graphic results illustrate practical cases with fluid extraction/injection into a reservoir using real data.


Introduction
Thermoporoelasticity is a branch of poromechanics that describes general thermo-hydro-mechanical phenomena (THM) [1], [2] for real-world processes occurring when the reservoir rock is subjected simultaneously to geomechanical, thermal, hydraulic and other physical effects. The linear poroelasticity theory of M. Biot [3], [4], [5], [6] is isothermal; it uses Hookean classic elasticity to describe the mechanical response of the rock, coupling Darcy's law to model the fluid transport within pores and fractures subjected to different types of stresses and boundary conditions. The processes involved in geothermal reservoirs can be isothermal or non-isothermal. In the first case, local thermal equilibrium (LTE) is assumed. The second case occur when the reservoir temperature exhibits changes, which can be assumed in LTE or under local thermal non-equilibrium conditions (LTNE). This occurs, for example, when liquid is injected at lower temperature into a hot reservoir; then liquid and solid phases interact through a volumetric heat transfer mechanism [7], [8]. Darcy's law is influenced by rock deformation because there are changes in porosity and permeability when pressure or temperature changes. Concerning thermal diffusion, it is assumed that rock strains and darcian flow have little effect in pure heat conduction. Both situations, LTE and LTNE are formulated in this paper, but only the uncoupled radial thermoporoelastic model is solved exactly in LTE. The more general     , , , , , , , , , , , , , , Where p f , ζ f , v D , u r , v S , ε r , ε θ , ε B , φ, σ r , σ θ , T f , T S , are fluid (pore) pressure, variation of the liquid content, Darcy velocity, solid displacement, solid deformation velocity, radial, tangential and bulk strains, porosity, radial and tangential stresses, fluid temperature and solid rock temperature, respectively. If only radial coordinates are considered, all variables involved in the isothermal radial model are functions of radius and time f (r, t). In non-isothermal conditions, the rock and liquid temperatures can be obtained as analytic functions of the radius only in a pseudo-stationary state. The vector displacement of the solid rock has only one component acting in the radial direction and, therefore its rotational is zero (figure 1). All symbols, functions and procedures to obtain the solutions are completely defined within the text, all their analytical formulae are in the Appendix.

Governing equations and radial analytic solutions under isothermal conditions (LTE)
The vector displacement of the solid particles u r (r, t) has only one spatial component acting in the radial direction and, therefore, its rotational is zero. The fluid content variation inside the pores ζ f (increment or removal) can be measured experimentally or computed from poroelastic formulae [9], [10]. Under the condition of null rotational of the vector radial displacement, the variation of the fluid content ζ f (r, t) becomes proportional to the fluid pressure p f [9], [2]: Where S r is the uniaxial specific fluid storage in the radial direction, which is calculated in terms of the radial compressibility and other poroelastic coefficients defined below. If the fluid pressure satisfies the geometry and hypothesis of the infinite continuous line source (see figure 1), then it can be computed approximately with the modified Theis solution:  Figure 1. Reservoir radial geometry showing its main elements and thermodynamic functions. Vector e r indicates radial direction, r is the radius, t is current time and z w is the length of the liner.
Where p 0 , Q V ,  f , z w , k r , η f are initial reservoir pore pressure, volumetric fluid flow rate, fluid viscosity, liner's length, absolute radial permeability, and hydraulic diffusivity, respectively. The symbol E 1 is the exponential integral of degree 1. From this classical solution (2) Darcy velocity v D is deduced (gravity is neglected): The uniaxial specific fluid storage S r in the radial direction is calculated in terms of the radial compressibility 1/K v and other poroelastic coefficients defined as follows: Where K B , is the rock drained bulk modulus, G is the shear coefficient, γ e is the poroelastic loading efficiency, ν U is the undrained Poisson's coefficient, B is the Skempton coefficient, and b is the Biot-Willis modulus [6]. The fluid pressure satisfies the geometry and hypothesis of the infinite continuous line source (figure 1), it can be computed with the modified Theis solution, from which Darcy velocity v D is deduced. The model defined by equations (1, 2, 3 and 4) is valid for homogeneous reservoirs of radial geometry, with a fully penetrating well in a very large isotropic porous medium of thickness z w . The main hypothesis is that the reservoir is initially a single-phase system, with uniform pressure and uniform temperature everywhere.

Thermoporoelastic governing equations and temperatures in general form (LTNE)
The basic equations governing the radial behaviour of a linear thermoporoelastic rock are deduced from general principles and physical laws well established [1], [2], [3]. [4], [5], [6], [9], [10], [11]. The main general principle is the equilibrium equation in poroelasticity: Where σ T is the total stress tensor acting in the fluid-rock system. The strains in the solid skeleton are defined in terms of the components of the solid vector displacement u = u i (x i ) e i : Where ε is the strain tensor and x i represents any kind of coordinate, Cartesian, radial, etc. The equation relating stresses and strains is: Where ε B (= ε xx + ε yy + ε zz ) is the volumetric strain, δ ij is the unit tensor, λ and G are the Lamé and shear coefficients respectively for drained conditions, p 0 , T 0 are the reservoir initial pore pressure and initial temperature, b is the Biot-Willis coefficient [6], K B is the bulk modulus, and γ B is the bulk thermal moduli; all coefficients are defined in the next section. Substituting equation (7) into the equilibrium condition (5) and using equation (6), we obtain the first governing linear thermoporoelastic formula: Using the law of mass conservation and Darcy's law, the second governing equation for the fluid flow is obtained: Under the most general conditions of local thermal non-equilibrium state (LTNE, T s ≠ T f ), the heat transfer process between phases is modelled by two partial differential equations, one for the solid phase (s) and one for the fluid phase (f ): where t, φ, c, ρ, T, k, v D , q s and q sf are time, porosity, specific heat capacity, density, temperature, conductivity tensor, Darcy velocity and volumetric heat generation of the solid (s) and fluid (f ) phases respectively. The symbol q sf is the amount of volumetric heat transferred from the solid matrix to the fluid and viceversa [10], [11]; this term is a function of the solid-fluid temperature difference (T s -T f ); the volumetric heat q sf depends also on a thermal coefficient and geometric variables. The microscopic fluid velocity in the pores v f , is related to Darcy flux v D by the Dupuit-Forchheimer relation v D = φ v f [2].
The general model (equations 5 -11) include the energy transfer between the solid and the fluid phases at different temperatures; they correspond to a set of coupled equations governing general thermo-hydro-mechanical phenomena (THM) in porous rocks in LTNE. The main unknowns are , the vector displacement of the solid particles u, the strain tensor ε, the fluid pressure p f and the fluid and rock temperatures T f and T s respectively. However, this linear formulation corresponds to a partial thermally uncoupled model, because both temperatures can be solved independently of the other unknowns if fluid velocity, porosity and initial and boundary conditions are known.
Where γ B is the bulk thermal expansivity defined in next subsection. Porosity is defined as a bilinear function of the fluid pressure and reservoir temperature [1], [10], φ 0 is the initial porosity.
Where N b is the tangent Biot modulus [1], γ φ is the thermal expansion coefficient of pores at constant p f and defined in next subsection. Physically, any applied radial stress σ r generates a strain ε r producing simultaneously a perpendicular tangential strain ε θ , which corresponds to a tangential stress σ θ (figure 1). The strains of the radial displacement are deduced from the volumetric strain, which is equal to the divergence of the vector displacement u r :

Theis model of the diffusion equation for the fluid content variation
The classical Theis model has the following initial, internal and external boundary conditions that are also satisfied by the fluid content variation ζ f : The function ζ f satisfies the same diffusion equation with appropriate initial and boundary conditions [9]: ζ f also satisfies the following relationship in polar coordinates [2], [9]: Where γ e is the coupled poroelastic coefficient defined in equation (4), η f is the hydraulic diffusivity defined in equation (2) and is assumed to be a function of the fluid temperature T f only. The boundary and initial conditions in equation (17) allows to integrate this equation exactly; the corresponding integration constants are eliminated, because u r must be bounded for any value of (r, t ): , The irrotational condition for the solid displacement, implies that the liquid content variation ζ f is proportional to the fluid pressure p f , which can be computed with the Theis solution [9], [10]: Theis model , 4 4 The viscosity of liquid water grows slightly when pressure decreases, but its variability is larger when temperature changes. Therefore, Theis model in equation (2) can be used for every isothermal curve with fixed porosity, and is approximately valid when both hydraulic diffusivity and porosity are updated for different fluid temperatures. It is interesting to mention here that S. Garg [11] derived a simple diffusivity equation for the two-phase flow of water in geothermal systems. His model is valid for reservoirs of radial geometry with the same Theis' hypothesis, which assume a fully penetrating well in a very large homogeneous, isotropic reservoir of thickness z w . The main hypothesis in Garg's model is that the reservoir is initially a two-phase system with uniform pressure and temperature everywhere. The resulting two-phase partial differential equation and its solution are completely analogous to equation (2) for liquid condition.

Computation of the radial displacement, strains and stresses
To obtain the solid radial displacement, the following procedure is used. The sudden injection of fluid at point (0, 0) is mathematically equivalent to the classic problem of an instantaneous line source with injection of heat at time t = 0, which satisfies the diffusion equation (17) The radial displacement produced by an infinite continuous line source, corresponding to the Theis solution, is obtained by integrating equation Once the radial displacement u r is obtained from equation (22), all the other unknowns can be computed directly with the relationships given by equations (12, 13 and 14) 1 e 2

Experimental thermoporoelastic coefficients
The coefficients introduced in previous sections are defined as follows. The variation of the fluid content ζ f can be experimentally measured as [3], [9]: where ρ f is fluid density, φ is current porosity and ρ 0 is the initial fluid density. The thermal expansion coefficients, bulk modulus and compressibility are defined as follows: Where C f , C B , are fluid and bulk compressibilities respectively, and K f is the bulk rock modulus. Biot moduli are defined next: Detailed developments of all these moduli are described in [9], [10].

Analytical solutions for temperatures in the pseudo-stationary LTNE radial model
Assuming a quasi-stationary state (∂/∂t ~ 0), equation (10) in radial coordinates T S (r) becomes: Where η S is the thermal diffusivity of the solid rock. Integrating twice, the analytic solution for this PDE is: Where Q S (°C/s) is the global solid heat transfer and δ S is thermal diffusivity in geothermal rocks. The conduction-convection heat partial differential equation for the liquid when the fluid flow is radial and quasi-stationary, with corresponding boundary conditions is: Where δ f is the thermal diffusivity of the fluid. Integrating twice and replacing the constant boundary conditions, the solution of this equation is: Where u w = u r (r w , 0); this solution is too large to include it here explicitly.

Poroelastic results from the isothermal radial model with fluid extraction
Assuming there are no temperature changes (T = T 0 ), the following graphics were obtained with the radial model herein described.   10 9 Pa γ e (ad)

ThermoPoroelastic results from the non-isothermal radial model with fluid injection
Assuming there are temperature differences during fluid injection, it is necessary to add the effect of the thermal stress. Once both strains are known, stresses and porosity are computed using equations (12) and (13) to include the effect of both pressure and temperature changes in the fluid/rock system. The following graphical results are obtained for an injection rate Q v = + 0.05 m 3 /s.

Discussion of graphical results of the non-isothermal radial model with fluid injection
The results obtained with the radial model herein presented, confirm several interesting experimental facts. First of all, fluid extraction or injection in a well produces poroelastic deformations (figures 3, 4, 6 and 7). Even if they are of small magnitude, the strains and solid displacements are numerically detected at very early times, in the order of seconds. This occurs because the liquid flows through a porous rock whose solid skeleton can be deformed elastically instantaneously. Fluid extraction/ injection in the reservoir causes the reduction/increment of the internal pore pressure (figures 2 and 5 left) affecting both, the liquid content and the effective porosity (figures 3 and 6). Their magnitude is proportional to the volumetric flow rate (equation 22). A wave of non-negligible amplitude appears in the porous rock immediately after the fluid extraction/injection begins (figures 4 and 7), but it is rapidly attenuated during the next simulation times. Therefore, the presence of the moving fluid in the porous rock modifies its mechanical response.
The internal tension produced by the liquid extraction (figure 2), induces a decline of the fluid pressure and of the pore fluid content (figure 3 left). The corresponding reduction of effective porosity (figure 3 right) is the principal source of liquid released from storage. When the poroelastic rock is subjected to internal compression because of injection, the resulting matrix deformation leads to a volumetric increase of the pores containing the fluid ( figure 6). This increment of the pore volume must be bounded by physical poroelastic limits, defining a transitional zone before the rock enters the non- Where P frac is the minimum pressure for the fracture to occur, σ θ is the previously defined tangential stress (equation 12 and figure 8 right) and p I is the extra pressure induced by the injected fluid. This breaking pressure depends on several factors, specifically on the injection rate.

Graphical solutions for temperatures in the radial model (LTNE)
In a pseudo-stationary state (∂/∂t ~ 0) for the temperatures of the solid and liquid phases, a Local Thermal Non-Equilibrium uncoupled radial model is obtained. Assuming that the temperature in the solid skeleton at r = r w is T w , and at distance r = r L >> r w is T L , the analytic solution for the solid T S (r) and fluid temperature T f (r) in radial coordinates for these boundary conditions are given by equations (30) and (32) respectively. The corresponding graphics of both temperatures are:   The pseudo-stationary state is only valid for short injection times measured in minutes, or few hours. The approximation extracted from these analytic solutions (30) and (32) is not valid for larger residence times of the fluid. Full numerical solutions are necessary in that case.

Conclusions
Assuming the hypothesis of the classic Theis solution for pressure in a radial reservoir, an analytic mathematical thermoporoelastic model was constructed. This model includes fluid flow, rock deformation and temperatures distribution. The total number of unknowns of the radial model are thirteen with the same number of equations. The solutions of the radial model presented herein are complete and self-sufficient. Interested readers could use them directly and get their own outcomes.
The results obtained show four interesting key facts: (a) -The transport of a fluid in the reservoir modifies its mechanical response. Fluid extraction or injection in a well produces poroelastic deformations. Even if they are of small magnitude, the strains and solid displacements are detected at very early times, in the order of seconds, because the liquid flows through a porous rock whose solid skeleton can be deformed instantaneously in elastic linear form. (b) -Fluid extraction or injection in the reservoir causes the reduction/increment of the pore pressure affecting both, the liquid content and the effective porosity in direct proportion to the volumetric flow rate. (c) -A wave of non-negligible amplitude appears in the vicinity of the well immediately after the fluid extraction/injection begins. This wave is rapidly attenuated. (d) -The internal tension produced by liquid extraction induces a decline of the fluid pressure and of the pore fluid content. There is a corresponding reduction of porosity that can be the principal source of liquid released from storage. Other conclusions follows: -It is well known that hydraulic fracturing can create high-conductivity paths within a large area of the reservoir. The pressure increment caused by liquid injection induces the rock formation to fracture hydraulically.
-The internal compression produced by fluid injection, causes a skeleton deformation that leads to a volumetric increase of the pores containing the fluid. This increment of porosity must be bounded by physical poroelastic limits, defining a transitional zone before the poroplastic region appears, where the rock can be fractured.
-This model can be useful to explore the start of a hydraulic fracturing process. Unpublished experimental data show that during a stimulation treatment performed in a low-permeability reservoir, which consisted in the injection of fluids at high pressure and flow rate into the formation interval of interest, a vertical fracture of 2.7 cm aperture was created after 300 seconds of continuous injection. The peak value of the injection pressure was 282 bar, which corresponded to the breaking pressure. After this peak, it was a pressure draw-down, which was stabilized at 160 bars during the next 10 hours of continuous injection.
-The breaking pressure required to induce fractures in a rock at a given depth can be estimated with this radial model if the maximum poroelastic radial displacement of the rock is known.
-This work is intended to explore the physical limits of a correct analytical solution of the thermoporoelastic problem in geothermal reservoirs. All the unknowns were successfully coupled in the isothermal poroelastic case. However, under local thermal non-equilibrium conditions the coupling is not possible using the exact model. The development of a numerical solution for the fully coupled non-isothermal case is a current research work in progress.