Thermo-hydro-mechanical Research Articles

Thermo-hydro-mechanical (THM) coupled simulation of the land subsidence due to aquifer thermal energy storage (ATES) system in soft soils

Aquifer thermal energy storage (ATES) system has received attention for heating or cooling buildings. However, it is well known that land subsidence becomes a major environmental concern for ATES projects. Yet, the effect of temperature on land subsidence has received practically no attention in the past. This paper presents a thermo-hydro-mechanical (THM) coupled numerical study on an ATES system in Shanghai, China. Four water wells were installed for seasonal heating and cooling of an agriculture greenhouse. The target aquifer at a depth of 74–104.5 m consisted of alternating layers of sand and silty sand and was covered with clay. Groundwater level, temperature, and land subsidence data from 2015 to 2017 were collected using field monitoring instruments. Constrained by data, we constructed a field scale three-dimensional (3D) model using TOUGH (Transport of Unsaturated Groundwater and Heat) and FLAC3D (Fast Lagrangian Analysis of Continua) equipped with a thermo-elastoplastic constitutive model. The effectiveness of the numerical model was validated by field data. The model was used to reproduce groundwater flow, heat transfer, and mechanical responses in porous media over three years and capture the thermo- and pressure-induced land subsidence. The results show that the maximum thermo-induced land subsidence accounts for about 60% of the total subsidence. The thermo-induced subsidence is slightly greater in winter than that in summer, and more pronounced near the cold well area than the hot well area. This study provides some valuable guidelines for controlling land subsidence caused by ATES systems installed in soft soils.

Coupled thermo-hydro-mechanical model for simulating thermal history-induced voids in buffer material

The thermo-hydro-mechanical processes in a deep geological repository (DGR) influence the bentonite buffer’s long-term performance. The buffer’s temperature evolution over a repository’s service life imposes a thermal history in the DGR. The developed thermal history cause shrinkage of the buffer and generates the technological voids (i.e., annular gaps) in a DGR. The void formed compromises the intimate contact between the dissimilar materials within the DGR and thus impedes the containment ability of the buffer in the long run. Therefore, to investigate the complex coupled processes involved, the present study developed a coupled thermo-hydro-mechanical (THM) model to evaluate thermal history’s influence on buffer’s hydro-mechanical performance. The developed THM model was validated using data from an infiltration test conducted under isothermal conditions. Furthermore, the THM simulations on thermal history indicated radial and vertical shrinkage within bentonite. The formed voids were sealed due to rehydration occurring through the bottom boundary. The hydro-mechanical parameters, such as swell deformation in radial and vertical directions, porosity, permeability, water intake, and swell pressure, were assessed during the hydration process. A sudden increase in the swell pressure was observed when the void gets sealed completely due to rehydration. The rate and magnitude of swell pressure generated were proportional to the peak temperature imposed at the non-yielding boundaries and the consequent void size formed. The increased swell pressures could have impaired the structural integrity of the waste contained canister at a DGR scale. Given this, the implications were comprehended based on the outcome of numerical investigations for the plausible adverse effects inside a DGR due to thermal history.

A semi-implicit material point method for coupled thermo-hydro-mechanical simulation of saturated porous media in large deformation

The coupled thermo-hydro-mechanical (THM) response of liquid-infiltrated porous media underpins the safe operation and maintenance of key engineering infrastructure. Challenges remain in modeling and understanding the complicated multiphysics processes of porous media subjected to THM loads and undergoing large deformation. In this study, we develop a stabilized material point method (MPM) for modeling the THM responses of large deformation problems in biphasic solid–fluid mixtures. A novel and efficient staggered solution scheme is proposed to solve the governing equations of the coupled system formulated in terms of four primary variables: solid displacement (u), liquid velocity (v), pore pressure (p), and temperature (T). The scheme solves the energy balance equation first and employs the resulting temperature to further advance the calculation of the momentum and mass balance equations using a semi-implicit fractional step method to facilitate equal-order interpolations. Both the incompressible and weakly compressible fluid are considered in the presented fractional step formulations. We also develop the axisymmetric form of the coupled MPM to increase the applicability and efficiency of the method in THM problems. The validity, stability, and robustness of the proposed method are demonstrated through three benchmark problems, including the heating of a saturated half-space, the non-isothermal consolidation of a soil column, and a three-dimensional axisymmetric problem pertaining to the thermoelastic response around a point heat source. The predictive capability of the proposed method for large deformation problems is further showcased by the simulation of the progressive failure process of a thermal-sensitive slope.

Tensile strength degradations of mineral grain interfaces (MGIs) of granite after thermo-hydro-mechanical (THM) treatment

Granite for hosting the deep geological repository (DGR) is subjected to thermo-hydro-mechanical (THM) treatment for decades of years, while the variations in the tensile strengths of mineral grain interfaces (MGIs) with THM treatment are still unknown. Served as an important parameter in determining the macro mechanical properties of granite, quantifying the tensile strength degradation of MGI with THM treatment is crucial for ensuring the long-term safety of DGR. In this study, the varied tensile strengths of MGIs in the granite after different THM treatments were first directly measured based on a self-developed single-MGI mechanical test system. Three-bending tests were conducted on six-hundred MGIs, which were submitted to various THM treatments with a temperature of 25–175 °C and a pore pressure of 0–15 MPa. The tensile strengths of the interfaces of quartz, feldspar and the interfaces between quartz and feldspar were derived based on the Dugdale–Barenblatt (D–B) model. Results show that the average tensile strength of the interfaces of quartz is 1.06 times larger than that of interfaces between quartz and feldspar, and 1.13 times larger than that of the interfaces of feldspar. Besides, the tensile strength of MGIs seems more susceptible to the pore pressure than to the temperature, and the degradation of tensile strength is the most obvious in the interfaces between quartz and feldspar, followed by the interfaces of feldspar and finally the interfaces of quartz after THM treatments. Based on the characteristics of the tensile strength distributions of different MGIs after THM treatments, the tensile strength degradation models of MGIs related to THM treatments were proposed in meso scale.

Transient temperature and water distributions in compacted MX80 bentonite under high temperature gradients

Numerous element-scale and full-scale tests have been carried out on compacted bentonite to investigate its coupled thermo-hydro-mechanical (THM) behavior when used as a buffer material in high-level nuclear waste geological repositories, such as the Full-scale Engineered Barriers Experiment (FEBEX) and its forerunner laboratory mock-up test [5]. In these past studies the applied temperatures were below 100°C, while new repository conditions being explored in the recently started HotBENT project at the Grimsel Test Site involve waste canister temperatures up to 200°C [2]. Recently, Lu and McCartney [3] conducted a tank scale test where a temperature of 200°C was maintained by a cylindrical heating element having a diameter of 12.5 mm at the center of a compacted MX80 bentonite layer. A notable temperate drop in temperature was observed with distance from the small diameter heating element. A similar phenomenon was also reported by Gens et al. (2020) for heating tests on MX80 pellets during heating to 140°C. Measurements from the FEBEX project indicate a more distributed temperature decay away from the heater [3]. This discrepancy between laboratory and field tests provided the motivation for better understanding the temperature and water distribution near a heater under high temperature gradients to better guide the design of laboratory-scale experiments in the future.
 A one-dimensional heating test was performed in this study to investigate the coupled THM response of a MX80 bentonite layer in unconstrained conditions during basal heating. Measurements include the transient temperature and volumetric water content redistribution at different vertical distances from the flat plate heater under high-temperature gradients as well as global volumetric strains. Three dielectric sensors (Decagon 5TM) were placed into a compacted MX80 layer having an initial gravimetric water content of 12.2% and a dry density of 1.15 Mg/m3 within a PVC modified Proctor mold (inner diameter of 152 mm and height of 178 mm), as shown in Figure 1(a). Low-density polyethylene sheeting was stretched across the top and bottom of the soil layers to minimize any global loss of water from the compacted bentonite layer during heating. A temperature-controlled hot plate was used as the basal heat source, and a dial indicator was mounted on the top of the soil layer to measure volume changes. Fiberglass insulation was wrapped around the cell to minimize lateral heat loss from the side boundary. A basal temperature of 128°C was maintained during a heating stage lasting 4000 h to ensure THM equilibrium, after which three short-term (96 h for each) cooling-heating cycles were applied. Results indicate that soil temperatures increase rapidly in the first 20 h of heating then gradually stabilize (Figure 1(b)). The volumetric water content increases at the beginning of heating then slowly decreases (Figure 1(c)). This increasing-decreasing trend is due to the upward movement of a wetting front away from the heater due to coupled heat transfer and water flow. An initial increase in volumetric strain (swelling) was followed by a decrease (contraction) as heating progressed (Figure 1(d)). This is consistent with the observed thermal response of normally consolidated expansive soils (e.g., [4]). A sharp monotonic decrease in volumetric strain occurred during the three short-term cooling-heating cycles. The temperature and water content follow nearly identical trajectories in each cycle. The temperature distribution in the bentonite with distance from the heating element after 4000 h of heating is smoother than the sharp drop-off in temperature from the small-diameter heater observed by Lu and McCartney [3], as shown in (Figure 1(f)). The main implication of the results presented in this study is that the surface area of the heater in contact with the surrounding bentonite plays a major role in the temperature distribution and the associated coupled water flow. This is a critical finding that may help guide the geometric design of small-scale laboratory mock-up experiments used to study THM processes in controlled laboratory conditions.

Numerical simulations of thermal fracturing in nuclear waste disposal

Regarding high-level nuclear waste (HLW) and spent fuel (SF) disposal, underground geological disposal is likely to be the most feasible and effective solution [1]. The coupled Thermo-Hydro-Mechanical (THM) behaviour of host rocks has been widely studied by laboratory tests, in-situ experiments and numerical analyses. In field experiments performed in underground laboratories, pore pressure generation and damage evolution have been observed. In this context, coupled THM numerical analyses can be an essential tool to better understand the coupled multiphysics behaviour of host rocks, buffer materials and canisters. In addition, large-scale field experiments provide a crucial opportunity to demonstrate understanding and to assess the predictive/modelling capabilities under realistic conditions.
 A number of hydromechanical (HM) numerical simulations of underground excavations have been carried out and the mechanisms of pore pressure and damage evolution have been studied [2]. Also, numerical analyses have been used to predict the observations of THM behaviour of in-situ heating tests, and there is some research on the damage evolution due to thermal loading [3], but few studies focus on the analysis of strain localization problems, in part due to specific numerical difficulties. However, thermal fracturing behaviour has indeed been observed in experiments and a better understanding of the phenomenon and of shear band evolution is crucial for a proper assessment of the stability of HLW/SF disposal.
 In this study, the coupled THM response of a generic HLW disposal facility under high-temperature has been analysed; in particular, the mechanism of thermal pressurization and the corresponding thermal fracturing behaviour have been examined. For this purpose, an advanced elastoplastic constitutive model with damage has been used to represent the mechanical behaviour of the host rock; it incorporates the following features: (1) a hyperbolic approximation to the Mohr-Coulomb yield surface to represent shear and tensile failure modes, (2) stiffness and strength anisotropy, (3) hardening and softening behaviour, (4) viscous effects and (5) a non-local integral type approach for the plasticity component of the model. The adopted overstress theory and the nonlocal formulation are capable of regularising the spatial distribution of plastic strains, eliminating mesh dependency and improving numerical efficiency. The coupling between accumulated plastic strain and hydraulic permeability has been considered in the analyses. Note that the term damage used in this study, refers to a state of the material, instead of the concept associated with damage mechanics theory; in addition, fracture represents localization phenomena in the damaged area.
 Numerical simulations have been performed using the finite element method software CODE_BRIGHT. Two-dimensional generic disposal configurations have been carried out for the two cases of unsupported and supported tunnels. The host rock is assumed to be the Callovo-Oxfordian claystone (COx). Three modelling stages have been considered: the excavation stage, the waiting stage and the heating stage, as shown in Figure 1 (a). In the heating stage, the thermal flow at the borehole wall is set as a constant power (200 W/m) and a no hydraulic flow boundary condition is prescribed. A maximum temperature of around 80 °C has been obtained on the drift wall at heating time t = 10 years.
 In addition, Figure 1 (b) shows the contours of plastic shear strains and liquid pressures at heating time t = 10 years. Near-field damage is caused by both the excavation of galleries and the effect of thermal pressurization. In the simulations presented, shear bands develop significantly during the heating stage. By comparing two cases in Figure 1 (b), it can be concluded that the installation of liners can efficiently limit the damage evolution and reduce the likelihood of thermal fracturing, significantly improving the stability of HLW disposal. Note that the scales of the two cases in Figure 1 (b) are different. Figure 1(c) shows that a lower pore pressure close to the gallery wall for the unsupported tunnel due to the higher level of damage in this case that results in a larger permeability of the rock in the damaged area.
 Apart from the analysis of the mechanism of THM behaviour in a high-temperature environment, numerical simulations have also been carried out of the large-scale PRACLAY Heater test (performed in the HADEs underground laboratory) and of the CRQ heating experiment (performed in the Meuse-Haute Marne underground laboratory). A good agreement has been observed between numerical predictions and experimental data, in terms of both temperature and pore pressure. These numerical validation exercises are, however, outside the scope of this abstract.

Thermo-hydro-mechanical modelling of geothermal energy extraction in deep mines in spatially heterogeneous settings

With the increasing demand for mineral and alternative energy resources, as well as the gradual depletion of shallow resources, the exploitation and utilization of mineral resources and geothermal energy in deep strata is an effective way to solve the problem of resource shortage [1]. In recent years, as a new type of resource mining mode, the co-mining of deep mineral and geothermal energy has developed rapidly [2, 3]. This method can make use of the original equipment of the mine for geothermal exploitation. However, the deep co-mining system faces two significant challenges: the first is the significant uncertainty inherent in subsurface properties, while the second is the high levels of geostress and temperature associated with deep mining. These challenges are adding some constraints on the practicality of exploiting such systems and limit the feasibility of deep resource co-mining, so that modelling efforts are needed for actual risk assessment.
 Consequently, we developed a Thermo-Hydro-Mechanical (THM) coupling framework for geothermal energy exploitation in deep mines using COMSOL to quantitatively characterize the temperature field of the geothermal system and predict the stress field of the mining system, considering the joint effects of large uncertainties and THM coupling. Through SGeMS, the uncertainty and spatial heterogeneity distribution of porosity are first generated. Then, the uncertainty of the hydraulic parameter [4] (permeability), mechanical parameter [5] (elastic modulus), and thermal parameter [6] (heat capacity and heat conductivity) was derived from the porosity. 500 samples were generated within a given uncertainty range, by means of Monte Carlo simulations. The spatial and temporal distributions of the temperature field of the geothermal system, and the stress field of the mining system were simulated, for each sample with COMSOL. Using the distance-based global sensitivity analysis [6], the most sensitive parameters for deep mining are identified, the heat storage capacity of the system and evolution of the maximum stress ratio are evaluated, including uncertainty.

Effect of climate change on the behaviour of thermo-active diaphragm walls

Energy geo-structures are becoming more common as a renewable energy solution which utilises shallow geothermal energy to provide heating and cooling to buildings and civil infrastructure projects. There has been an upward trend in the use of thermo-active geo-structures in the United Kingdom, from 150 in 2005 to over 5000 in 2012 [2]. Studies so far have shown that diaphragm walls subjected to combined thermo-mechanical loading show overall increases in lateral displacements, bending moments, shear forces, axial forces, and settlements on the retained side with thermal cycles [6, 7, 8].
 This study uses a variation of a validated numerical model [4] to predict the behaviour of thermo-active diaphragm walls in the longer-term including accounting for the influence of climate change. This numerical model also assesses the impact of different modelling assumptions on the model output by comparing a simplified model with a more complex model where atmospheric temperatures affecting ground surface temperatures are included. Research from the IPCC shows increases in ambient temperatures and in maximum extreme temperatures as a result of climate change. Two models are compared in the IPCC report [5], Representative Concentration Pathway (RCP) 2.6 where governments across the world work together to meet the recommendations in the report to reduce overall warming, and RCP8.5 where the current trends of heating continue without any governmental interventions. According to climate model RCP2.6, global mean surface air temperatures will increase by around 1.5°C and extreme daily temperatures (i.e. days where temperatures exceed expected temperatures for the time of year) will increase by around 2°C over the 50-year period chosen for this study, whilst using RCP8.5 global mean surface air temperature increases by approximately 3°C and extreme daily temperatures increase by approximately 5°C.
 This study uses the same model geometry as a model validated through centrifuge testing [4] which is a 5m deep cantilever diaphragm wall, like those that would be used in deep excavations for open-topped tunnels in large civil infrastructure projects. The numerical model was developed using the semi-coupled thermo-hydro-mechanical (THM) thermal module available with PLAXIS 2D and uses the HS-Small soil model [3] to allow for variation in shear modulus at small strains (i.e. those associated with retaining walls) to be captured more accurately. Thermal loads were applied by applying heat flows through a plate element with the same mechanical properties of the model diaphragm wall used in centrifuge testing. This model attempts to capture the influence of atmospheric temperatures by applying thermal boundary conditions to the ground surface as well as to the thermo-active diaphragm wall. Studies have demonstrated that a more complex model capturing atmospheric temperatures increases the impact of thermal loading [4], whilst also capturing a more realistic temperature profile throughout the entire depth of the soil body both adjacent to the wall and in the free field condition. In this study, atmospheric temperatures were captured in RCP2.6 by cycling between 30°C and -5°C over the course of the cycle to represent a typical ‘hottest’ and ‘coldest’ day of the year in cycle 1, increasing to 32°C and -5°C in line with the temperature increases described in the IPCC report [5]. Ambient mean temperatures between heating and cooling cycles increased from 20°C to 21.5°C over the 50-year period. In the same model, the temperature within the diaphragm
 wall cycled between 35°C and 10°C in cycle 1, representing the typical operation of a ground source heat pump (GSHP) in heating and cooling. These temperature inputs increased to 37°C and 10°C over the 50-year period. The RCP8.5 atmospheric temperatures cycled between 30°C and -5°C to represent summer and winter in cycle 1 and increase to 35°C and -5°C over the 50-year period. Thermal flows in the diaphragm wall cycled between 35°C and 10°C in cycle 1 and increased to 40°C and 10°C over the 50-year period. The results from this study show increases in lateral displacement, maximum bending moments, positive and negative shear forces and axial forces (compressive (-) and tensile (+)). These results are shown in Figure 1.
 Significantly, the RCP2.6 model shows that these increases begin to stabilise over the 50-year period, with the increases in the final cycles becoming negligible. With the RCP8.5 model these increases are increasing linearly at the end of the modelling period which can potentially have significant implications not only in terms of serviceability of the structure, but potentially on ultimate limit state as linear increases suggest that design limits will eventually be reached. The findings from this study generally agree with other similar studies, such as a study that analysed the geotechnical performance of energy piles with additional environmental loading as a result of climate change which found that additional compressive and tensile stresses are generated during heating and cooling phases, respectively [1].

Comparison of two approaches for modelling fracture opening due to cold water injection in geothermal reservoir

Heat and electricity production from deep hot reservoirs through hydrothermal or petrothermal system requires to inject cold fluid in a naturally or artificially fractured medium. Cold water injection in a hot reservoirs causes thermo-hydro-mechanical (THM) coupled processes that may have several influences on operations and long term production [1]. Cold water has a higher viscosity, which means injection can become difficult due to higher flow impedance, while thermal diffusion in the rock matrix causes shrinkage and thus an increasing fracture aperture over time that can lead to flow channeling and a reduction of operation time. A deep understanding of the contribution of coupled THM processes to injection into fractured media is thus important to predict the long-term performance of a geothermal power plant. In this study, coupled THM processes in a single fracture are numerically investigated with two different approaches to model the discontinuity.
 The first numerical approach proposes modelling an implicit fracture in a fracture zone using solid elements (Fig.1 (a)). The impact of the fracture is introduced by means of a stress dependent fracture aperture, which in turn determines the fracture transmissivity. This function allows recovery of the opening, and thus the model replicates opening and closure of the fracture. Then, an equivalent continuum permeability [m2] is obtained as follows:
 where [m] is the elastic opening, br and bmax [m] are the residual and maximum apertures, respectively, a is the stress dependency coefficient; [Pa] is effective normal stress, and [m]the fracture spacing. The combination of the continuum elements to model an implicit fracture embedded in a fracture zone and the evolving permeability has been used successfully used to model elastic opening in volcanic rock [2] and granodiorite rock [3, 4].
 In the second approach, zero-thickness interface elements are used to explicitly model the opening of an existing fracture (Fig.1 (b)), following the work by Liaudat et al. [5]. A constitutive law that is capable of describing the fracture initiation, propagation, closing and opening is implemented in this approach, considering heat transfer in the discontinuity. In this approach, the cubic law is used to estimate the longitudinal transmissivity [m3] as a function of the normal fracture aperture:
 where [m] is the normal separation and [m3] is a constant value which makes it possible to assign an initial longitudinal transmissivity even if the fracture is mechanically closed.
 In both approaches, the simulation of cold water injection requires to consider the heat and hydraulic flows in the fracture alongside mechanical behaviour, since changes in pore water pressure and temperature influence the fracture aperture, thus modifying the fracture transmissivity.
 In order to compare both approaches, a simple synthetic model is studied. The model with its boundary and initial conditions are shown in Fig.2. Water with a temperature of 20 °C and pressure varying from 4 MPa – 10 MPa is injected at the centering point of the fracture. Figure 3 shows the injection pressure profile. Fracture permeability evolution is compared in Fig.4, which shows the two different methods capture different characteristics of the natural fracture behaviour. In the continuum method, the fracture is modelled as a continuum element, thus any pressure or stress change will result in corresponding changes in the fracture permeability, which is reflected in Fig.4 showing that the permeability changes with the injection pressure changing, In contrast, in the interface element, Fig.4 shows only when the injection pressure reaches 10 MPa, the permeability increases sharply from the initial value to the peak value. This is because penalty method is used to allow the small overlapping of the two faces of the discontinuity, thus keeping negligible permeability variation, before it opens. The result in Fig.4 also illustrates the effect of thermal stress on the fracture opening, considering the total vertical stress is 10 MPa while injection pressure is below or equal to 10 MPa.
 In conclusion, though the peak permeability predicted from two methods are consistent, both around 1.5 x 10 -12 D, the behaviour before fracture opens in these two approaches are quite different. For the interface method, permeability variations keep negligible before fracture opens, while for the continuum method, permeability changes with corresponding pressure or stress changes without dependence on the opening of the fracture.

Direct shear setup for studying the THM shear response of the soil-structure interface in the context of energy geostructures

Developing pathways for climate-resilient development involves integrating mitigation and adaptation actions, ensuring sustainable development [1]. A climate-resilient development can be achieved through the inclusion of effective mitigation approaches into development planning, reducing vulnerability, conserving ecosystems, and restoring ecosystems [1]. In this regard, energy geostructures are introduced as an effective mitigating approach, providing renewable energy while limiting the emission of greenhouse gases [2-4]. The operation of the system is closely linked to daily and seasonal cycles, which leads to cyclic temperature and water content fluctuations at the soil and the soil-structure interface [4, 5]. Thus, the study of the shear response of the soil–structure interface subjected to different thermo-hydro-mechanical (THM) conditions is of importance. Testing techniques used to study the THM behaviour of unsaturated soils require advanced laboratory equipment, as well as protocols for correcting measured data due to errors in the test conditions and apparatus calibration. This paper presents the development of a new direct shear setup to measure the non-isothermal shear strength of the partially saturated soils and soil-structure interface. The modified setup, a unique one to the authors' knowledge, enables simultaneous control of temperature, matric suction, and mechanical stress state within the soil specimens. The operational temperature range of energy geostructures (i.e., 5°C to 50) is applied through a thermal plate, developed from corrosion-resistant stainless steel with high thermal conductivity, placed at the base of the soil specimen. Matric suction (i.e., in the range of 0 to 100 kPa) is controlled using the axis-translation technique and measured using a pressure transducer connected to the back of the top cap, incorporating the HAE disk, facilitating the measurement (Figure 1) [2, 4-6]. The direct shear setup is modified to accurately measure the shear strength and deformation characteristics of soil samples under controlled laboratory conditions. The design of the device is based on previous direct shear devices but includes several improvements to enhance its accuracy and ease of use [4, 5]. The device has been tested to measure the shear response of both soils and soil-structure interfaces, and the results were compared to those obtained using conventional direct shear devices [2, 4, 7-10]. The results indicate that the new device is accurate and reliable and represents a significant advancement in the field of soil testing.
 In this study, the interface is formed by kaolin clay, a temperature-sensitive clay, and concrete, a structural material widely used for energy geostructures [4, 11]. The soil samples were prepared by static compaction at 30% initial water content () and an initial void ratio () of 1.2. Initially, all samples were inundated with distilled water at room temperature (24°C) [4, 12], followed by sequential hydraulic (i.e., matric suction, = 0 or 70 kPa), mechanical (i.e., net normal load, =100 or 300 kPa), and thermal (i.e., temperature, =24°C or 45°C) loading of the interface. The shearing was initiated after ensuring equilibrium criteria were met. It was necessary to limit the shearing rate to 0.005mm/min and the thermal loading rate to 3°C/hr to maintain the drained condition [4, 13]. Additionally, the vertical displacements associated with hydraulic load were subjected to an equilibrium criterion of 0.025%/day strain rate [14]. As presented in Table 1, the results revealed that the apparent interface friction angle was not significantly affected by matric suction at varying temperatures, but a slight decrease has been observed upon heating at all matric suctions. The apparent adhesion increased in response to temperature increase/decrease, with a decreasing rate as suction increased, while the interface desaturation led to higher apparent adhesion at all temperatures, corresponding to peak and residual values. In interfaces subjected to identical normal stress and temperature but different matric suction values, the shear stress-shear displacement curves showed greater peak stress with increasing matric suction. Furthermore, despite the same net normal stress and matric suction, higher temperatures resulted in slightly lower peak shear stress and less contractive volume change behaviour at the interface. Non-isothermal volumetric behaviour and soil dilatancy play a significant role in governing the THM shear response of the interface [15].

Numerical simulation of frost heave of saturated soil considering thermo‐hydro‐mechanical coupling

AbstractFrost heave can lead to both the ground uplifting and frost heave pressure under different circumstances, and cause many engineering problems. To describe the characteristics of frost heave under various thermo‐hydro‐mechanical (THM) coupling conditions and calculate both the frost heave amount and the frost heave pressure, the coupled THM process, as well as the phase change of the pore water, should be considered for freezing soil. In this paper, thermodynamic equilibrium conditions in saturated freezing soil were derived to account for the mechanical effect on the cryogenic suction and unfrozen saturation of pore water during phase transition. Then the governing equations were developed considering the water flow, heat transfer, stress equilibrium, phase change, ice segregation, and their complex interactions. Based on that, a numerical model considering fully THM coupling is presented within the framework of finite element provided by COMSOL. Both frost heave amount tests and frost heave pressure tests were simulated to verify the proposed model. Then the model was applied to subgrades with different types of soil and different boundary conditions, to reveal the characteristics of frost heave amount and frost heave pressure under different conditions, and at the same time to demonstrate the model's applicative prospect.

Thermo-hydro-mechanical coupled model of unsaturated frozen soil considering frost heave and thaw settlement

Frost heave and thaw settlement of subgrade is the most typical engineering disease in permafrost regions, which seriously threatens the safe construction and operation of transportation infrastructure. On the basis of the Modified Cam Clay model, a thermoelastoplastic constitutive model of unsaturated frozen soil is proposed. This model can effectively describe two key issues of the mechanical behavior of frozen soils, that is, soil hardening caused by temperature reduction and volume compression during thawing. The proposed constitutive model is realized in the established thermo-hydro-mechanical (THM) coupled mathematical model, and the reliability of the coupled model is verified by comparing with the experimental results. In addition, the proposed THM framework is applied to describe the long-term freezing and thawing behavior of typical railway subgrade sections in the permafrost region of western China. The results show that water-heat transfer and deformation frequently occur on the subgrade surface as the seasons change. Moreover, the variation of thaw settlement is greater than that of frost heave, which seriously threatens the operation safety of railways in permafrost regions.

Response Prediction of Asphalt Pavement in Cold Region with Thermo-Hydro-Mechanical Coupling Simulation

Although the theoretical pavement structure design method (TPSDM) is widely used for designing asphalt pavements in Japan, it still exhibits certain limitations, such as not considering the variation in moduli of the base and subgrade layers due to water contents, freeze–thaw action, and stress states. This study aims to enhance Japanese TPSDM’s accuracy by considering variations in the resilient modulus of environmental impacts, pavement materials, pavement structure, and traffic load actions to accurately calculate the mechanical responses and predict pavement fatigue life. Firstly, the study develops a 3D Thermo-Hydro-Mechanical (THM) model using the finite element method (FEM) to investigate temperature and moisture distributions of the pavement with time. Then, based on the numerical results of the moisture, temperature, and stress state obtained from the THM analysis, the constant resilient modulus of the base and subgrade layers in the Japanese TPSDM is replaced with a resilient modulus that considers the stress state and the combined effects of water content fluctuations and freeze–thaw action. Finally, the fatigue life of the pavement is calculated based on the obtained mechanical response in THM analysis. The reliability and validity of the proposed fatigue life prediction method are well verified by comparing the calculated with the actual pavement fatigue life. Results indicate that the modifications improve the Japanese TPSDM by considering the environmental impacts, traffic load actions, pavement materials, and pavement structure, thereby improving the accuracy of predicting the fatigue life of asphalt pavements, particularly in cold regions.

Numerical study on the effects of thermoelastic and poroelastic parameters on the geomechanical behaviour of Hot Dry Rock geothermal reservoirs

ABSTRACT During the last two decades, energy production has been directed towards renewable resources, one of the most important of which is geothermal energy. Despite the reliability of geothermal energy as a resource, its effects on the surrounding environment have not been investigated in detail. This study focuses on reservoir deformation as one of the environmental concerns of geothermal energy extraction. A coupled thermo-hydro-mechanical (THM) model is employed for simulating the reservoir deformations. The validity of the model results is examined by comparing the numerical results with the analytical solution of a benchmark THM problem. The model is then utilised to study the behaviour of a Hot Dry Rock reservoir consisting of two wells. The obtained results show that in some cases the reservoir deformation is significant, making it an important factor in design considerations. Also, it is found that the deformation caused by thermal volume change is up to 6 times more significant compared to that of poroelastic effects. The conducted parametric study demonstrated that the coefficient of thermal expansion of rock and the production rate severely influence the reservoir deformation, while rock elasticity modulus, porosity and Biot-Willis coefficient only affect the behaviour of the geothermal reservoir to some extent.

An equivalent thermo-hydro-mechanical model for gas migration in saturated rocks

Gas transport process in saturated rocks is gaining increasing attention in the environmental/energy geomechanical engineering fields, e.g., shale gas exploitation, carbon dioxide sequestration and nuclear waste disposal, etc., as it is closely associated with gas production rate, sealing capability of caprocks. To investigate gas migration behavior through saturated rocks which induces microcracking inside the material, an equivalent thermo-hydro-mechanical (THM) model is developed that considers coupled behavior of fractures. The fractured rock material is considered as the equivalent continuum, which consists of isotropic porous matrix and parallel fracture set. The mathematical formulations are constructed within the thermodynamic framework, where the THM constitutive relations are proposed. Compared to the model without considering fracture behavior, the equivalent model improves the calculation accuracy of about 60% against laboratory experiments (e.g., observed gas flow rate), which proves the model capability. Furthermore, the sensitivity analysis shows that highly-pressurized gas induced microcracking has significant influence on the coupled process, which greatly contributes to the observed gas flow phenomenon.

A framework to model the hydraulic fracturing with thermo-hydro-mechanical coupling based on the variational phase-field approach

In this research, a numerical framework for hydraulic fracturing has been formulated, incorporating thermo-hydro-mechanical (THM) coupled effects. While previous studies have reported various hydraulic fracturing models based on the phase-field method, a THM coupling scheme grounded in the variational phase-field approach remains unexplored. The THM coupling is of paramount importance for understanding underground fracture propagation. This study’s principal innovation lies in its extension of the variational phase field model, introducing it to a thermoporoelastic medium, thereby integrating fluid flow and heat transfer elements in fractured materials. In our advanced model, the fluid flow and heat transfer in the matrix and fracture are modeled independently. Subsequently, employing a set of reasonable assumptions, we introduce unified flow and heat transfer equations derived from phase-field calculus. A pivotal contribution of our work is the development of an iterative solution algorithm based on the fixed stress method, pivotal for tackling the THM coupled system. The proposed equations are discretized employing the FEM. Additionally, to address the numerical oscillations in advection-dominant issues that manifest in the heat transfer model, we have applied a generalized Streamline Upwind Petrov Galerkin method for 8-node elements. The accuracy of fracture propagation was assessed by comparing the numerical results of the evolution of fracture length, fracture width at the injection point, and the injection pressure under the K regime with the corresponding analytical solution. Furthermore, the process of heat transfer within a fractured domain, often referred to as the Barends problem, has been studied by proposed model. Additionally, some numerical experiments were conducted to gain insights into the intricate interplay between the thermal–hydraulic–mechanical system during hydraulic fracturing, thereby further emphasizing the significance of our advanced model.

Landslide susceptibility evaluation in Alpine environment: 2. Thermo-hydro-mechanical modeling for the response to climate-related variables

This paper is Part 2 of two companion papers, proposing a multidisciplinary approach to assess stability and velocity evolution of a large landslide located in the Central Italian Alps (upper Valtellina region): the Ruinon landslide. Part 1 of this work presented a 3D stress–strain finite element analysis, which assessed the morphological and geomechanical predisposition of the slope to gravitational instabilities and defined the current stress state along the slope. In this paper, a thermo-hydro-mechanical (THM) numerical analysis is applied to the landslide shear zone, to assess the link between landslide driving factors and the shear band material response. Data used as input for the model were pore pressure, reference stresses and initial temperature at the sliding surface, as well as the monitored velocity of the landslide body, assumed to move as a rigid block. The shear band material was modeled as a visco-plastic medium with thermal softening and velocity hardening, thus thermal- and load-rate sensitivity of the material were estimated through laboratory testing. To this end, triaxial compression tests with thermal control were performed on rock samples representative of the shear band. To constrain the model, results of the analysis presented in Part 1 were used to define the stress state at the sliding surface and the relationship between pore pressure and shear stresses. Then, pore pressure data from in-situ piezometers relevant to the period 2014–2018 were introduced and a best fitting between modeled and monitored landslide velocities was obtained. Finally, velocities were forecasted for the period 2018–2020​ and a process of validation was performed using field displacement data. The outputs of the model adequately simulate the measured landslide velocity, reproducing the sliding behavior and its relationship with pore pressure. The presented approach may be applied to further case studies, aimed at defining a novel physics based early warning strategy for landslides.

Numerical modeling on nuclide transport around a nuclear waste repository under coupled thermo-hydro-mechanical condition

Coupled thermo-hydro-mechanical (THM) processes may raise the uncertainty of nuclide transport behavior in long-term operation of nuclear waste repository. A far-field modeling framework on nuclide transport around a nuclear waste repository under coupled THM conditions was developed and verified. Base on conceptual models of high-level waste repository in fractured rocks, a three-dimensional repository model including a vertical fault and a horizontal fracture was built, and influence of coupled THM processes on nuclide transport is numerically investigated. Effects of repository depth, fracture inclination, and internal aperture heterogeneity were investigated using the developed model. The results show that the coupled THM processes significantly affect nuclide transport behavior, especially in the early stages of nuclide leaching. The thermal stress induced by heat release of nuclear waste causes the reduction in fracture aperture around the repository. The nuclide transport under THM conditions is retarded before 103 years than that under hydraulic condition. The deeper locations of the candidate repository exhibit the lower rate of nuclide transport. The inclined fracture induces opposite flow directions between two fractures, which causes the delaying of nuclide transport. Fractures with heterogeneous apertures have preferential flow paths, which results in early arrivals of nuclide transport than fractures with homogenous apertures.

Open-source Python library for modeling coupled thermo-hydro-mechanical (THM) processes

Abstract. The interactions between temperature, fluids, and mechanical properties in a repository system are essentially described scientifically by coupled thermo-hydro-mechanical (THM) processes. THM modeling, i.e., the prediction of the behavior of materials under different conditions, is the fundamental numerical tool here. The confident handling and deep understanding of numerical computation methods is thus the prerequisite for performing and evaluating preliminary safety analyses in the site selection process. At the Federal Office for the Safety of Nuclear Waste Management (BASE), a software library for the simulation of coupled THM processes is currently being developed. The main goal of the in-house development is an open-source toolbox in the scripting language Python and is motivated by several long-term sub-goals: targeted development of expertise within BASE regarding numerical modeling of safety-relevant aspects in the long-term safety analyses; diversification of the testing capabilities regarding the preliminary safety investigations by means of an in-house, independent modeling software; foundation of a library of known benchmarks and evaluation examples for the comparison of different software tools; documentation and processing of basic THM scenarios for internal or, if necessary, public technical training. The focus of this development is on creating a toolbox that is easy to use and at the same time highly flexible. The main methodical aspects are as follows: building a new library based on the pyGIMLi pre- and postprocessing framework (Rücker et al., 2017); creating a finite-element reference implementation in the Python scripting language for maximal transparency; creating an easy-to-use interface to the solution of the weak formulation for the finite-element theory with expressions of a symbolic manner allowing maximal flexibility; defining an interface to allow for the integration of alternative, third-party high-performance libraries; creating a collection of Jupyter notebooks of well-documented test cases and benchmarks. Choosing the open-source approach ensures the best possible transparency and, in the long term, also allows the provision of appropriately quality-assured and documented simulation tools to the public. The presented poster shows the current development status of the software library and the currently implemented quality assurance concepts and gives an outline of the potential applications of the library.

Landslide susceptibility evaluation in Alpine environment: 1. 3D Finite Element modeling of the Ruinon (IT) case study

This manuscript is Part 1 of two companion papers that explore a multidisciplinary approach to predict velocity and stability of a large landslide located in the Central Italian Alps: the Ruinon landslide. The area is of high geological interest due to the presence of numerous shallow and deep gravitational instability processes that affect valley flanks, mainly driven by unfavorable morphological conditions and geomechanical properties of rock masses. In this manuscript, a 3D finite element model (FEM) was implemented in order to analyze the stress–strain distribution along the Ruinon rock-slope. Goals are to define the relation between morphological factors, mechanical parameters and the development of irreversible strains. The model was defined based on morphological features and mechanical properties detected along the slope, as well as on piezometric data from the landslide monitoring system. In a first step of the analysis, a static simulation was carried out under dry conditions. Then, a validation process was performed by comparing numerical outputs with geomorphological field observations. Finally, a parametric analysis was carried out where different piezometric level scenarios were evaluated in order to assess the influence of both mechanical parameters and pore pressure on the distribution of high sliding susceptibility areas. By overlaying satellite images with the model outputs, results were shown to accurately reproduce the extent of the slope areas subject to active gravitational instability. Parametric analyses showed a clear relationship between the input factors and the magnitude of strain, while the extension of areas subject to irreversible deformation did not change significantly. Stress distribution and stress–strain relations defined in this article are subsequently introduced in a thermo-hydro-mechanical (THM) numerical model presented in Part 2 of this work, where the evolution of the Ruinon landside is simulated.