A simple computational framework for THM processes in saturated geomaterials

Investigations of thermo-hydro-mechanical (THM) processes in geomaterials have become a continuing concern due to the increasing interests in the application of geothermal energy. In the last two decades, a number of THM coupled constitutive models were proposed for soils and rocks based on either the elastoplastic theory or the thermodynamics. The practical application of the constitutive models inevitably requires computational approaches such as finite element or finite difference methods. In this paper, a simple computational framework for THM processes in saturated geomaterials is proposed. An effort is made to explain the details and major concerns related to the governing equations of three physical fields. The computational framework is validated through the benchmark problem of 1D THM processes in saturated linear elastic material by comparing with the analytical solution. We then further investigate some key factors of the multiphysical processes and their influence on the time-dependent evolution of excess pore pressure, effective stress and strain. The flexibility of the framework allows the effect of some coupling components, such as the temperature-dependent dynamic viscosity of pore fluids, which are often ignored in previous studies, to be examined. The framework serves as the basis for the next-step in-depth analysis of theoretical multiphysical problems and practical design of energy foundations.


Introduction
Saturated geomaterials, such as clays, sands and rocks, are porous media composed of solid and liquid phases. A change of boundary temperature or internal body temperature will trigger a series of complex internal mechanical and seepage mechanisms. Such scientific problems are usually called the thermal-hydraulic-mechanical (THM) coupled problems. The systematic study of the THM coupled problems started from 1980s [1]. Recently, due to the increased application of energy geostructures, nuclear waste disposal and CO2 storage, many related experimental and numerical work have been carried out [2].
To numerically investigate THM coupled problems, it is necessary to first develop coupled governing equations based on the laws of each physical field and then assume one constitutive relation for each simulated material. In recent years, Many studies focused on the development of coupled 2 constitutive equations to investigate elemental thermo-mechanical soil behaviors (e.g., [3,4]). Meanwhile, only a few studies considered both the governing equations and the constitutive models. For example, Di Donna and Laloui [5] applied ACMEG-T model in the FE software LAGAMINE to investigate the behavior of energy pile foundations during heating-cooling cycles. Rui and Soga [6] used an advanced critical state model in a FE code developed at the University of Cambridge to simulate the thermal response of a 23 m long thermal pile installed in heavily overconsolidated London Clay. Tamizdoust and Ghasemi-Fare [7] applied an thermo-elastic model in FE software COMSOL to investigate the heat and mass transfer in deformable porous media in a transient and quasi-steady state conditions.
In fact, some temperature-related factors in the governing equations may have a fundamental impact on the THM coupled soil behaviors. These factors include transient or steady state heat transfer, nonlinear temperature-dependent behavior of soil permeability, thermal expansion of pore fluids, etc. For instance, from 20 to 90°C, the thermal expansion coefficient of water triples its value, and the dynamic viscosity of water drops to about one third [8]. The former significantly increases the generation of excess pore pressure upon heating, and the latter leads to an increase in seepage velocity at elevated temperature. Such concerns, although important, are quite often simplified or overlooked in previous studies, especially if the constitutive models used are complicated.
The goal of the present paper is to introduce a simple computational framework and to numerically investigate the THM processes in saturated geomaterials. Unlike several previous studies focusing on the performance of constitutive models, this study tries to explain the important factors related to the governing equations. The paper first introduces the general computational framework, with special focus on temperature-dependent parameters and state variables. Then, the framework is validated through the benchmark problem of 1D THM processes in saturated linear elastic material by comparing with the analytical solution. Finally, parametric studies are conducted to highlight the importance of temperature-dependent parameters and state variables for the accurate prediction of stress-strain and excess pore pressure evolution within saturated geomaterials.

Computational framework for THM processes
In this section, the derivation of the governing equations of the three physical fields will be illustrated separately. The forms of the governing equations, if applied to single-field analyses, are classic and universal. For multiphysical field problems, however, the governing equations consist of coupling terms, which control the accuracy and rigor of performed analyses.
This study assumes instantaneous temperature equilibrium between local solids and pore fluids, and no phase change is assumed for pore fluids under elevated temperature. The formulation takes the compression positive sign convention.

Mechanical governing formulation
The fundamental equation of the mechanical field for saturated geomaterials is the equilibrium equation.
where ij is the total stress tensor,  the bulk density and gi the gravity tensor.  can further be expressed as  = nl + (1 − n)s, where n is the porosity, l the density of liquid phase and s the density of solid phase.
Equation (1) is only suitable for extreme soil conditions, such as fully drained or undrained conditions. For other cases, effective stress principle should be introduced. where ij' is the effective stress tensor, p the pore water pressure, i the Biot's coefficient affecting the plane perpendicular to the vector i, and ij the Kronecker delta. A Biot's coefficient of unity can be applied to sands and clays, and values of less than unity are more suitable for hard rocks.
For thermo-mechanical problems, the incremental form of the total strain tensor can be expressed as where ∆ij  and ∆ij T are the incremental strain tensors due to mechanical stress change and temperature change, respectively. Both strain tensors can further be separated into elastic and plastic parts. The thermally induced elastic strain, ∆ij T,e , is generated by the principle of thermal expansion and related to the volumetric thermal expansion coefficients of the solid phase, s. It is generally believed that the thermally induced plastic strain, ∆ij T,p , is caused by the physicochemical interactions between clay particles and the adsorbed water [9][10][11]. Recently, it is discovered that another mechanism contributes to the thermally induced plastic strain of both clays and sands, which is the increase in granular packing caused by the thermal expansion of granular particles [12,13].
Here we should introduce a constitutive model to generate a relationship between the effective stress and the strain tensors. A general form of the incremental relationship should be where Eijkl is the constitutive matrix that depends on the specific constitutive model. For modern nonlinear elastoplastic or thermodynamic models, Eijkl is generally a function of both the strain tensor and the density. Under small strain assumption, the connection between the density and the total strain tensor is expressed by the mass continuity equation where n0 is the initial porosity.

Hydraulic governing formulation
For the hydraulic field, the vector of seepage velocity, vi, can be expressed by the Darcy's law where vi t is the relative velocity vector of liquids to solids,  the intrinsic permeability and  the dynamic viscosity of the pore fluid. The intrinsic permeability is a function of the porosity and independent of temperature. In the meantime, both l and are very sensitive to temperature change. Therefore, it is very important to take into consideration the temperature dependence of l and if the expected temperature variance is larger than 10 º C. Only a few recent studies have taken such consideration in their THM coupled analyses [7,14,15].
The classic continuity equation of the liquid phase is expressed as where vi l is the velocity vector of the liquid phase and t stands for time.
Similarly, the continuity equation of the solid phase is expressed as where vi s is the velocity vector of the solid phase. It should be mentioned that s is also temperature dependent following the principle of thermal expansion.
Note that div(vi s ) = − ∂kk / ∂t, equations (6～8) can be summarized into Equation (9) includes the deformation rate of the solid skeleton and the liquid phase due to thermal expansion, forced convection and free convection. Further derivation of equation (9) requires explicit expressions of the temperature-dependent state variables, such as s and l, and those of the temperature-dependent parameters, such as . Some reference data can be found in [15,8].

Thermal governing formulation
Heat transfer in geomaterials usually involves three mechanisms: conduction, convection and radiation. For geotechnical applications, the following two assumptions are adopted: 1) the effect of radiation is assumed negligible, and 2) free thermal expansion or contraction of a porous material do not change its porosity. Thus, the energy conservation equation reads as where T is the temperature; ceq = (1 − n)scs + nlcl is the averaged specific heat; cs and cl are the specific heat capacity of the solid and the liquid phase, respectively;  = (1 − n)s + nl is the averaged thermal conductivity; s and l are the thermal conductivity of the solid and the liquid phase, respectively. The first term in equation (10) is the increment of internal energy, the second term is the heat transfer by conduction (Fourier's law), and the third term is the heat transfer by convection, which shows thermal-hydraulic coupling. The thermal conductivity is independent of temperature. But it may vary with the change of external stresses and the internal microstructure. It should be kept in mind to check that the unit of the thermal conductivity is compatible with other parameters, otherwise the calculation result may be unrealistic.

THM coupled governing formulation
Adopting appropriate boundary conditions and standard Galerkin finite element discretization method, the incremental form of the fully coupled THM formulation can be expressed as where K1 and K2 represent the governing matrices, the first vector is the rate of nodal unknowns, the second the nodal unknowns, and the third the nodal external forces. The formulation can be solved using proper time integration schemes, based on the complexity of the problem.
It should be noted that, for THM coupled analyses, the temperature-dependent state variables are: ij T,e , ij T,p , s, l. The temperature-dependent parameters are: the thermal expansion coefficients of the solid and the liquid phase, and the dynamic viscosity of pore fluids. If the pore fluid is further separated into adsorbed water and free water, the above parameters should be modified accordingly. Since the adsorbed water is often considered to be completely adsorbed on the surface of the solid particles, it can be treated as part of the solid skeleton. As a result, the modified thermal expansion coefficient of the solid skeleton should be different from the intrinsic thermal expansion coefficient of the solid material.

Comparison with analytical solution
To verify the proposed computational framework, a viable way is to compare the results of a case study with those of analytical solutions or commercial software. In this study, the example of onedimensional THM processes in saturated linear elastic geomaterial is adopted because a rigorous analytical solution for this case has been proposed by Suvorov and Selvadurai [16].  Figure 1 presents the FE mesh and the boundary conditions of the present case. Consider a 1D soil column of 10 m in height. The bottom end is fixed, and the top end has no displacement restraints. The initial temperature in the soil column is 0°C. When thermal loading begins, the whole column is instantaneously heated to 100°C without drainage. As a result, the top surface rises, and excess pore pressure builds up in the column. Then, a constant compression stress of 10 MPa is applied on the top surface, and the temperature and the pore pressure at the top boundary are instantaneously reduced to zero. Table 1 shows the list of all the parameters used for the analysis. The proposed framework is simplified to be compatible with the assumptions of the analytical solution. The heat convection is not considered herein. Figure 2 presents the evolution of temperature, excess pore pressure and displacement in the 365 days after loading. Results show good agreement between the finite element and the analytical solution.
Since the heat convection is not considered, the temperature evolution with time and depth in figure 2a is controlled solely by the heat conduction in equation (10). An increase in the thermal conductivity or a decrease in the specific heat will accelerate the process towards temperature equilibrium. Due to the similarities in the thermal and the hydraulic governing equations, the temperature evolution in this case is similar to the excess pore pressure evolution in the classic 1D consolidation theory considering one-way drainage.
In figure 2b, after the 10 MPa compressive stress is applied to the top, the maximum value of the excess pore water pressure is about 36.3 MPa. This suggests that 26.3 MPa of the total excess pore pressure is caused by the 100 ºC temperature rise, which can be explained in the following way. In the current problem, the solid phase has a non-zero thermal expansion coefficient, and the thermal expansion coefficient of the liquid phase is set to be zero. Assuming that displacement is allowed on the side boundary of the column (which means the problem becomes two-dimensional), the temperature rise will cause the liquid phase to be stretched, and negative pore pressure will be generated. Meanwhile, due to the lateral constraints in 1D condition, the liquid phase is squeezed rather than stretched, so the excess pore water pressure is in fact positive. Note that the total vertical stress of the top surface remains zero throughout the process. Therefore, when the temperature is raised to 100ºC, the vertical effective stress of the top surface becomes -26. 3 MPa. This suggests that the solid phase is vertically stretched, which is one of the causes of the upward displacement (as shown in figure 2c). In the meantime, due to the lateral constraints, the solid phase is compressed in the radial direction, and the radial effective stress is about 5.5 MPa.
In figure 2c, the maximum displacement of the top surface caused by the thermal expansion of the solid skeleton is about 0.018m. Later, with the decrease of temperature and the consolidation process, the displacement gradually decreases. In this example, the temperature change has a major influence on the displacement evolution. Due to the linear elasticity assumption of the constitutive model, the temperature-induced displacement is reversible. After the overall equilibrium of temperature and pore pressure, the final displacement of the column will be determined by the soil stiffness and the compressive stress applied.

Effect of dynamic viscosity of pore fluids and intrinsic permeability
The above numerical analysis assumes a constant hydraulic permeability to be in consistence with the analytical solution. This means that the dynamic viscosity of pore fluids, , and the intrinsic permeability,, in equation (6) are assumed constant. In fact, it is widely accepted that  is a function of temperature and  is a function of porosity. Since the variation of  and  was not considered in many previous THM coupled studies, it should be noted herein that the consideration can be important for the accurate prediction of thermal consolidation and heat convection processes. For the present case, the variation of  is negligible, because the maximum increase in porosity is less than 1%. Meanwhile, the variation of  can be significant within the range of 0～100°C. This study adopts the relationship between  and T proposed by Guvanasen Based on equation (12), the value of  at 100°C drops to about one-sixth of the value at 0°C, and its influence on the THM processes is shown in figure 3. b Figure 3. Influence of temperature dependence of dynamic viscosity on (a) pore pressure and (b) displacement evolution Since the heat convection is not considered, a change of  will not affect the temperature evolution. Meanwhile, the decreased  at elevated temperature greatly enhances the hydraulic permeability. As shown in figure 3a, compared with the analysis using constant , the consideration of temperaturedependent  accelerates the excess pore pressure dissipation at the bottom by 8% and 82% at day 1 and day 10, respectively. Accordingly, as shown in figure 3b, the thermally induced primary consolidation at the top accelerates by 6% and 13% at day 1 and day 10, respectively. This is why the thermally enhanced prefabricated vertical drains (PVD) can be used to accelerate the consolidation of reclaimed lands [18]. It is recommended that the temperature dependence of  be always considered in analyses with temperature variation larger than 10°C.

Effect of thermal expansion of liquid phase
The thermal expansion of the liquid phase is not considered in the above analyses. In practice, however, the excess pore pressure cannot be accurately predicted without a proper assumption for l. In this section, a same analysis is repeated using the parameters from    Figure 4 shows the numerical results at the time after imposing the 10 MPa compression stress but before releasing the top thermal and hydraulic boundaries. l ranges from -5.4×10 -5 at 0°C to 7.8× 10 -4 at 100°C, which is on average an order of magnitude larger than s. Therefore, the excess pore pressure in figure 4b reaches almost 1000 MPa, which is 28 times the previous result. The thermal expansion of the liquid phase also leads to an 0.14 m elevation of the top surface, which is 8 times the previous result. Due to the extremely high thermally induced excess pore pressure, the vertical effective stress is large and in tension, and the impact of the 10 MPa compression stress becomes insignificant. This is why heating a saturated rock with low-permeable pores may easily result in hydraulic fractures.

Conclusions
This paper presents a simple and generalized computational framework for THM coupled processes in saturated geomaterials. The purpose is not to put forward particular theories for particular materials, but to illustrate the major concerns in the governing equations regarding the THM coupled processes. The applicability of the framework is demonstrated by the benchmark problem of 1D THM processes in saturated linear elastic material and by comparing with analytical solutions. Results from the current study confirm that the temperature dependence of the thermal expansion coefficient and the dynamic viscosity of pore fluids highly affect the overall thermal performance of the saturated porous medium. It is therefore recommended that changes in parameters and state variables with temperature and pressure must be carefully considered whilst performing THM coupled modeling.