A mathematical framework for modelling 3D coupled THM phenomena within saturated porous media undergoing finite strains

Abstract A mathematical formulation of a coupled thermo-hydro-mechanical model for saturated porous media undergoing finite deformations is presented. The model, developed by unifying some of the most relevant works on poromechanics, introduces a reformulation of the stored energy function, the inclusion of thermo-osmotic effects as well as the adoption of a permeability tensor depending on deformation. Analytical solutions are derived to qualitatively evidence the features of the proposed approach, showing a good agreement with similar results reported in literature. Different thermal osmosis scenarios are additionally evaluated, particularly referring to negative thermal osmosis, positive thermal osmosis and no osmosis, so predicting the relevance of such a parameter when analysing water transport and mechanical deformation of porous media.


Introduction
Modelling coupled heat-moisture flows within deformable saturated/partially saturated porous media is a widely investigated field since the pioneer works by Biot [1][2][3] for establishing complete theoretical frameworks able to represent multi-component systems consisting of air, vapour, water and solid skeleton in non-isothermal conditions. Correspondingly, it has been possible to reconstruct coupled processes of heat-fluid flow and deformation within saturated/partially saturated media, as e.g. in [4][5][6] to name a few. Anyway, these works generally neglect [7] the influences of pressure gradients on heat flow (thermal-filtration which is analogous to Dufour effects in solutions) and the temperature gradient on fluid flow (thermo-osmosis effect which is analogous to Soret effect in solutions) [8,9]. Particularly the latter is referred as the so-called mechano-caloric effect which has appeared to be of relevance in argillaceous media [10,11]: the counterpart of the mechano-caloric effect derives from Onsager's reciprocity principle [12]. This thermodynamic principle requires that the reciprocal process exists, i.e. a fluid flow caused by a temperature gradient. Such a thermo-osmotic process has been observed in clay-rich materials (see e.g. [13][14][15]). Thermal and even chemical osmotic flows are in fact two of the main coupled flows caused by heat and ion transport in clay-rich media [16]; the overall direction of water flow by thermo-osmosis depends on the entropy difference between water in the clay pores and that of water trapped in the diffuse interface layer or clay interlayers [12,17,18].
Thermo-osmosis might contribute to convective transport of unsorbed radionuclides in the context of nuclear waste confinement in argillaceous media [12]. In general, to ensure nuclear safety, the waste packages are designed to delay the radionuclide release beyond the maximum heat period and no radionuclides are expected to reach the boundaries of a confining shale unit in less than 100 ky [19]. Thermo-osmosis becomes significant with very low hydraulic conductivity (10 -10 -10 -14 m/s, typical of clay barriers for waste disposal) in presence of temperature gradients [7]; the effect can generally be neglected in traditional engineering applications since the hydraulic conductivities are usually higher and the temperature gradient is often negligible. Particularly, by assembling the most exhaustive data including the parameters necessary for the analysis on thermo-osmosis available in literature, Gonçalvès et al. [12] confirmed that thermo-osmosis may prevail over the classical Darcy flow driven by pressure gradients, in clay-rich formations under natural or artificial thermal gradients [10,20]. Hence, given the contribution of this flow to convective transport of dissolved species, thermo-osmosis has been understood of being fundamental to describe and explain the pressure distribution and flows in shale layers, as well as to predict transport phenomena in such media, especially for nuclear-waste storage projects.
Among the most complete thermal-hydro-mechanical (THM) models accounting for thermoosmosis we recall Zhou et al. [7] for saturated, linear-elastic porous media in small strains, Ekbote & Abousleiman [21] presented a coupled poro-chemo-thermo-elastic model, Chen et al. [22] and Chen et al. [23] for unsaturated, elastic (the former) and elastic-plastic porous media in small strains, whereas Gonçalvès et al. [12] developed a mathematical expression for thermo-osmotic permeability. Chen et al. [24] extended the modified mixture theory to non-isothermal conditions, deriving a fully coupled thermal-hydro-mechanical model within the framework of finite deformations, and Yang et al. [25] presented a thermo-poroelasticity theory fulfilling local thermal non-equilibrium condition with thermo-osmosis effect for a saturated medium. A simplified formulation of the THM model by Thomas and He [26] for unsaturated soils is presented in Zagorščak et al. [27] including the effects of thermo-osmosis on hydraulic behaviour for saturated clays neglecting the mechanical formulation.
In this work a coupled three-dimensional THM model in finite strains [28] based on the modified mixture theory [24,[29][30][31] is developed from a previous HM version [32] for saturated porous media. Particularly, the upgraded model includes a permeability tensor dependent on deformation of the solid skeleton, it accounts for a modified expression of the stored energy function and a thermoosmosis contribution within the Darcy's law for fluid flows.

Balance equations for the porous medium
By assuming incompressibility of the solid and liquid phases, the balance of linear momentum for a saturated porous medium can be written as [24,29,30] where ρ = ρs(1 -φ) + ρwφ is the saturated mass density of the solid-fluid mixture (φ being the porosity) and σ the Cauchy total stress tensor, with being σ' effective stress, p pore pressure and B = Biot's coefficient [33] with κ and Ks drained bulk modulus of the solid skeleton and of the solid grain, respectively; the Biot's coefficient is generally assumed equal to 1 being Ks >> κ.
In the total Lagrangian formulation, Eq.(1) is rewritten in the reference configuration via the Piola transformation [34], i.e.
where P is the total first Piola-Kirchhoff stress tensor obtained from pull-back of σ and J = det(F), being F the deformation gradient of the solid skeleton (see Section 2).
Introducing the symmetric total Kirchhoff stress tensor τ = Jσ and the effective Kirchhoff stress tensor τ' = Jσ', it follows that [35] which represents the effective stress equation (2) in Lagrangian form.
Correspondingly, in saturated non-isothermal conditions, the mass balance equation for the solid can be directly written as which is identical to the expression derived by Borja & Alarcón [29] and may be used to describe the material evolution of porosity as a function of the Jacobian when the solid constituent is incompressible.
The local fluid content continuity equation in the current configuration can be written as [36][37][38] the modified Darcy's velocity ( § 2.2) in non-isothermal condition and incompressibility of the skeleton material (Eq.(6)), Eq. (7) can be rewritten as The Piola transform of w is with the Piola identity so the balance of fluid content (7) in the Lagrangian configuration reads which becomes in case of incompressible constituents, negligible thermal expansion of the solid grains and thermal expansion coefficient of the mixture coinciding with the standard thermal expansion coefficient in infinitesimal range. Additionally, as reported by Sun [31], if 0   w  , by applying the Piola transform, Eq.(15) reduces in isothermal conditions to

Balance equations for heat transport
It is here considered that all phases of the saturated porous medium are locally in thermal equilibrium, and hence the temperature of both solid and fluid constituents are identical locally (in a homogenized sense of each elementary representative volume), i.e. Ts = Tw = T. Except the additional advection term, the local balance of energy is analogous to that of the single-phase thermo-plastic material [31,40,41], hence with R the heat source term and c being the specific heat capacity per unit volume of the porous where J is the determinant of the deformation gradient (see below), ρw0 and ρs0 the initial densities, and cs and cw the specific heat capacities (per unit mass) of the solid and fluid constituents, respectively.
It is to be noticed that in Eq. (18) the contribution from the structural heating and dissipation is neglected, as commonly done for many petroleum and geotechnical engineering applications (see e.g. [42][43][44]).
The Cauchy heat flux q in Eq. (18) is often written as the dot product of the volume averaged heat conductivity tensor and the gradient of temperature [45], i.e.
with kT the volume averaged heat conductivity tensor.
For sake of brevity the balance of energy in the reference configuration is here omitted.

Solid model
The multiplicative decomposition of the deformation gradient F(X,t) is here employed [34,46,47] so that , the former occurring at a constant (fixed) temperature T along with the stress producing deformation gradient FM and the volume change JM, the latter characterizing the intermediate configuration.
In case of a mechanically incompressible material (i.e. no volume change during an isothermal process), JM = 1.
We assume a thermally isotropic material, so that the deformation gradient FT may be given by an isotropic tensor, e.g., where F(T) is a scalar function determining the volume change relative to the reference configuration and α = α(T) is the temperature dependent expansion coefficient.
It is to be noticed that effective stress-based constitutive models, such as the Cam-clay model and its enhancements, are developed without regard for the presence or absence of water in the voids of the soil. Models of this type typically replicate the phenomenological responses of dry soils, or those of water-saturated soils deforming under fully drained conditions. Additionally, together with material frame indifference, the approach is restricted to isotropy [29].
Recalling the Helmoltz free energy function [40] and following the lines of [48][49][50], we reformulate the strain energy function as where b e is the left elastic Cauchy-Green strain tensor, , and ξ is an internal plastic variable defined such that characterizes the hardening response of the soil; m  represents the purely mechanical response and can be identified as the conventional isothermal strain energy density function associated e.g. to the Cam-Clay model; t  is the purely thermal portion of the Helmoltz free energy density, whereas tm  represents thermo-mechanical effects, based on the assumption that the primary coupling is through volumetric expansion. The stored energy function in the reference state for simplicity is here assumed to vanish [49].
If now in Eq. (24) the dependence of m  on ξ is dropped to allow for the introduction of an empirical hardening law such those used in Cam-Clay models [51], we can write and, together with the assumptions above, being e A  the elastic logarithmic principal stretch in the A-th direction.
The thermal contribution reads being T0 the initial uniform temperature experienced by the medium.
On the assumption that thermal effects in shear (i.e. deviatoric effects) can be neglected relative to thermal effects in dilatation [49], the purely mechanical effect is here equated with the deviatoric term in Eq.(29), hence The above form of m  is valid only within a certain range of elastic deformations, and where the stored energy is a convex function of the elastic logarithmic principal stretches. Now with "elastic" entropy denoted by η e , ψ satisfies the state equations directly coming from the dissipation inequality [40].
The two-invariant yield function has the form where Pc is the Kirchhoff pre-consolidation pressure, M the slope of the critical state line, and P and Q the mean normal and deviatoric Kirchhoff stress invariants, respectively with β the vector of principal Kirchhoff effective stresses whose components A β provide for the spectral decomposition of the Kirchhoff effective stress tensor in Eq. (32) and e ε vector of the elastic logarithmic principal stretches (Eq. (26) (40) so that the yield function assumes dependency on thermal loads via P, i.e. ) , ), ( ( The evolution equations coming from the principle of maximum dissipation read being v L the Lie derivative,  the Lagrange multiplier satisfying the Kuhn-Tucker conditions and ηp representing an internal variable associated to the entropy production p  p e      (42) contributing to the thermal dissipation p T . It is to be noticed that the plastic entropy and correspondent thermal dissipation does not explicitly appear in the temperature evolution equation Eq.(18).

Fluid flow constitutive model
It has been found [36] that Darcy's law can be formulated as a constitutive representation of the "drag"-interaction between the phases; the upgraded form of Darcy's velocity (Eq.(2)) accounting for the thermo-osmotic potential [12, 22-25, 27, 31] results in with rl k relative permeability, μ dynamic viscosity of water, kT thermo-osmosis coefficient and k = k(J) saturated permeability tensor [35] being D the effective diameters of the grains, φ0 the initial porosity and 1 the second order identity tensor.
It is to be noticed that in Eq.(43) another coupled flow called thermo-diffusion or Soret flow, by which solute will diffuse with the effect of temperature gradient, should be added [53]; Soret flow is different from Fickian flow induced by concentration gradient and it is to be accounted for when studying contaminant transport in underground water. Anyway, such an effect is neglected in the present study.

One-dimensional analytical solution
 given by the constitutive equation (32) Correspondingly, neglecting the spatial variations for the fluid density, the balance of fluid content (12) reads in which the thermo-osmotic conductivity has been referred with Th k to avoid confusion with heat conductivity.
In absence of heat sources and assuming negligible contributions of pore pressure gradients on the thermal field, Eq. (18) reduces to With the above simplifications Eq. (48) becomes the typical parabolic heat conduction equation with constant nonhomogeneous boundary conditions, whose solution is given by being n a the Fourier sine series coefficients of ) ( ) ( x with the characteristic (diffusive) time scale in the problem given by t * = L 2 /α ; it can be proved that the first term of the series dominates the sum of the rest of the terms when t' = t/t * ≥ 1/π 2 i.e., with the data of Table 1, t ≥ 24 days. By posing T* = (T-T1)/(T0-T1) the solution of the heat balance equation simplified into Eq. (48) gives the spatial and temporal distributions of Fig. 1 and 2, respectively.
Following such an scheme, the heat and fluid balance equations appear nearly decoupled, being just the latter depending on the solution of the former; hence two different approaches have been followed: through the first one the deformation of the solid phase has been maintained as variable via the standard infinitesimal meaning and the contribution of the fluid pressure neglected; through the second one the pressure term was kept and v   rewritten, as reported below.
Hence, the first rearrangement of Eq. (47) leads to solve a thermo-mechanical problem in which infinitesimal deformation are driven by temperature, i.e.    Recalling that J x v ln     and being possible to demonstrate that, for the assumed temperature range, 1 ) (  T  , the constitutive relationship (46) can be rewritten as Eq. (55) in incremental form is introduced within the equilibrium equation (1) and, by which can be treated as a nonhomogenous problem of the type, being u(x,t) = p(x,t), with Q(x,t) representing a source term, and initial and boundary conditions given by u(0,t) = u(L,t) = 0, u(x,0) = Φ(x) = p0. The solution of the problem (58) is given by By comparing the pressure results coming from the two approaches (the first reducing the problem to a thermo-mechanical scenario in small strains, the second one treating the fluid balance equation in the most general way roughly accounting for a finite strains problem structure), Fig. 4 was obtained in which the underestimation of the pressure contribution from the infinitesimal approach clearly appears. Evidently, such a result is primarily given by the fact that pressure is fictitiously obtained via an equilibrium with stresses, but the main result fits satisfactorily with that from [38].
Additionally, the trend of the analytical solution accounting for the proposed approach agrees similarly with that reported within the cited reference (Fig. 5).  Additionally, the effect of thermal osmosis within the developed model has been analysed; starting from the value reported in Table 1  in absence of thermal osmosis the driving force of water flow is only the water pressure gradient, whereas in the latter situation water flows from lower temperatures to higher temperatures if there is no water pressure gradients (negative osmosis) [24].
When there is no thermal osmosis ( Figure 6), water flows from higher pore pressures to lower pore pressures, showing a monotonic decreasing pressure profile, typical of a diffusion curve. When the coefficient is set to positive, the thermal osmotic effect comes into play and causes influx of water into the porous medium, so producing a pore pressure rise. Differently, under a negative value water flows from lower to higher temperatures, so that the porous medium experiences a water efflux and a faster pore pressure decrease in the beginning of the analysis. The same phenomenon has been evidenced in [24]. Figure 6. Effect of thermo-osmosis conductivity on pore pressure distribution.

Conclusions
A unified mathematical framework for modelling saturated porous media in finite strains has been presented, characterized by a reformulation of the stored energy function, the inclusion of thermoosmotic effects as well as the adoption of a permeability tensor depending on deformation. The model incorporates the influences of the geometrical non-linearity on the full coupled solid deformation, pore-fluid diffusion, and heat transfer processes. Analytical solutions have been derived to show, although qualitatively, the relevance of a nonlinear approach, even demonstrating a satisfactory agreement with similar results reported in literature. Additionally, the contribution of thermal osmosis has been evaluated when accounting for negative thermal osmosis, positive thermal osmosis and no osmosis, so predicting the relevance of such a parameter when analysing water transport and mechanical deformation of porous media. A three-dimensional variational formulation is currently being prepared for subsequent implementation within a Finite Element research code.