Coupled Heat and Mass Transfer and Thermodynamic Response Characteristics of Surrounding Rock in Underground Tunnels

The development and construction of deep underground engineering projects have been progressively extended into the complex stratigraphic formations and deep geological environments of western China. Consequently, the issue of high-geothermal hazards in deep-buried tunnels has become increasingly prominent, representing a core challenge to the safe construction and long-term serviceability of deep geotechnical engineering [1, 2]. The in-situ high-geothermal environment of deep strata, combined with the disturbances induced by tunnel excavation unloading and groundwater seepage, gives rise to intense heat transfer and moisture migration behavior within the surrounding rock, thereby triggering compound engineering thermal hazards [3, 4]. In numerous high-altitude deep-buried tunnel projects across China, conditions involving the coexistence of extremely high-temperature surrounding rock and high-temperature groundwater have been observed. This severe thermal-humid environment not only significantly compromises the operational safety and work efficiency of construction personnel, but also persistently erodes the tunnel lining structure, accelerating its cracking, aging, and failure. Simultaneously, it induces internal stress redistribution and fracture propagation within the surrounding rock, substantially degrading its overall stability [5, 6]. The heat and mass transfer process in tunnel surrounding rock is a quintessential nonlinear, multi-physical field coupling evolution process. Thermal driving forces govern internal heat conduction and fluid convection, modifying the fundamental thermal properties of the rock mass and the occurrence state of groundwater. Seepage flow transports heat across spatial domains, altering the effective stress state of the surrounding rock through variations in pore pressure. The dynamic adjustment of the stress field further reshapes the fracture structure and permeability characteristics of the rock mass, thereby exerting a feedback control on the transport pathways for both heat and fluid [7, 8]. The coupled effects arising from this dynamic multi-physical field interaction directly govern the damage evolution and instability processes in high-geothermal tunnel surrounding rock. A systematic elucidation of the coupled heat and mass transfer mechanisms and the associated thermodynamic response characteristics in deep-buried tunnel surrounding rock is therefore essential. Such an understanding serves not only to advance the fundamental theory of multi-field coupling in deep rock masses, but also to provide a reliable theoretical foundation and engineering basis for the prevention and control of thermal hazards, ventilation optimization, and support system design in high-geothermal tunnels [9, 10].

Extensive theoretical, experimental, and numerical investigations have been conducted on the coupled thermal-hydraulic-mechanical behavior of fractured rock masses, progressively establishing a multi-field coupling analysis framework suitable for underground rock engineering [11, 12]. In early studies, the local thermal equilibrium assumption was commonly adopted, wherein the fractured rock mass was treated as an equivalent homogeneous continuum, with the temperature of the rock skeleton and the internal fluid assumed to remain consistent. This approach effectively simplified the solution of multi-field coupling problems and proved applicable to tunnel engineering analyses under conventional temperature conditions [13, 14]. To accommodate the complex conditions inherent in deep high-geothermal engineering, the local thermal non-equilibrium theory has been progressively applied to heat transfer studies in rock masses. By distinguishing the temperature difference between the rock skeleton and the fracture fluid, this theory enables an accurate representation of the heat exchange process between the solid and liquid phases under high-temperature conditions, effectively remedying the limitations of the traditional thermal equilibrium theory [15, 16]. Through sustained development, substantial progress has been achieved in the theoretical modeling, numerical simulation, and engineering application of multi-field coupling in rock masses [17]. Nevertheless, significant systematic deficiencies remain within the existing research framework. Most thermal non-equilibrium theoretical models are restricted to the analysis of a single heat transfer process, without incorporating the mechanical response mechanism into a unified coupling framework. Consequently, these models are incapable of elucidating the stress evolution and surrounding rock damage mechanisms induced by solid-liquid temperature differentials, and a complete mathematical characterization system for fully coupled thermal-hydraulic-mechanical processes remains absent. Furthermore, prevailing numerical simulation studies commonly prescribe core thermal properties, such as thermal conductivity and specific heat capacity of the rock mass, as fixed constants [18], while neglecting the nonlinear evolution characteristics of these parameters under coupled thermal and stress effects. However, laboratory experiments have confirmed that variations in temperature and pressure conditions significantly modify the thermophysical properties of rock masses. The adoption of constant parameter values can introduce considerable deviations in the computed thermodynamic responses, severely compromising the accuracy of numerical analyses [19, 20]. In addition, current thermal hazard assessment approaches for high-geothermal tunnels remain relatively simplistic, relying solely on dry-bulb temperature for safety evaluations. This practice neglects the synergistic effects of multiple factors—including humidity, airflow velocity, and thermal radiation—under high-temperature and high-humidity conditions. As a result, an accurate characterization of the true thermodynamic state of the surrounding rock and the on-site thermal hazard level cannot be achieved, impeding the implementation of refined hazard classification and scientifically based prevention strategies [21, 22].

To address the deficiencies in theoretical modeling, parameter characterization, and engineering evaluation prevalent in existing studies, a systematic and integrated research framework is established. An in-depth investigation is conducted on the coupled heat and mass transfer and thermodynamic response characteristics of surrounding rock in high-geothermal tunnels, yielding multiple core innovative contributions. A dual-medium fully coupled thermal-hydraulic-mechanical numerical model is developed based on the local thermal non-equilibrium theory, enabling the synergistic solution of multi-physical field variables and elucidating the microscopic thermodynamic mechanisms of heat and mass transfer in surrounding rock under non-equilibrium conditions. The coupling effects of temperature variation and stress have been embedded into the model, establishing a dynamic bidirectional feedback mechanism between thermophysical properties of the rock mass and the multi-physical fields, thereby overcoming the limitations of conventional constant-property modeling. Furthermore, a comprehensive classification and evaluation system for tunnel thermal-humid environments has been established by integrating multi-dimensional thermodynamic evaluation indicators, facilitating a scientific and refined assessment of thermal hazards in high-geothermal tunnels.

This study is organized into six chapters, with the overall research logic progressing in a hierarchical and coherent manner. Chapter 1 provides a review of the engineering background and the current state of domestic and international research, identifying the deficiencies in existing studies, and delineating the core research content and innovative contributions of this work. In Chapter 2, key thermophysical parameters of the surrounding rock are obtained through laboratory testing, with the evolution of these parameters as functions of temperature and stress systematically revealed and corresponding constitutive relationships formulated. Chapter 3 presents the derivation of the multi-field coupled governing equations and the construction of a complete theoretical model for coupled heat and mass thermodynamics. Chapter 4 establishes the numerical solution methodology, with model reliability validated through multiple sets of comparative tests and engineering monitoring data. In Chapter 5, parametric analyses are conducted under various working conditions using the validated model, identifying the dominant controlling factors and evolution patterns of the thermodynamic response of the surrounding rock. Optimization recommendations for engineering prevention and control are proposed in conjunction with the classification and evaluation system. Chapter 6 synthesizes the core conclusions of the study, summarizes the application value of the research findings, and outlines directions for future in-depth investigations.

Based on Biot's poroelasticity theory, with theoretical extensions made to accommodate the structural characteristics of the dual-medium rock mass, a coupled thermal-hydraulic-mechanical stress solution framework was established, suitable for fractured surrounding rock in high-geothermal environments, incorporating both the pore pressure effects of the dual-medium and the mechanical effects of thermal deformation of the rock mass. The total stress field of the model is governed by the combined action of the skeleton effective stress and the equivalent pore pressure, expressed as:

$\sigma=\sigma^{\prime}-\alpha_B \bar{p} I$             (15)

where, $\text{ }\!\!\sigma\!\!\text{  }\!\!'\!\!\text{ }$ is the effective stress tensor of the rock mass; ${{\text{ }\!\!\alpha\!\!\text{ }}_{\text{B}}}$ is the Biot effective stress coefficient; $\text{I}$ is the unit tensor; and $\bar{p}$ is the weighted mean pore pressure of the dual-medium, calculated from the fracture volume fraction, pore pressure, and fracture pressure through weighted averaging. In strict adherence to the physical nature of thermal deformation in rock masses, thermal expansion strain is generated exclusively in the solid skeleton, with the fracture fluid not participating in the thermal deformation process. The corresponding thermal strain expression is given as:

${{\text{ }\!\!\varepsilon\!\!\text{ }}_{\text{T}}}\text{=}{{\text{ }\!\!\beta\!\!\text{ }}_{\text{T}}}\text{(}{{\text{T}}_{\text{s}}}\text{-}{{\text{T}}_{\text{0}}}\text{)I}$                    (16)

where, ${{\text{ }\!\!\beta\!\!\text{ }}_{\text{T}}}$ is the coefficient of thermal expansion of the rock mass. By combining the linear elastic constitutive relationship with the static equilibrium conditions, the complete governing equations for the mechanical field were established:

$\text{ }\!\!\sigma\!\!\text{  }\!\!'\!\!\text{ =D:( }\!\!\varepsilon\!\!\text{ -}{{\text{ }\!\!\varepsilon\!\!\text{ }}_{\text{T}}}\text{)}$                                    (17)

$\nabla \cdot \text{ }\!\!\sigma\!\!\text{ +}{{\text{ }\!\!\rho\!\!\text{ }}_{\text{b}}}\text{g=0}$         (18)

where, $\text{D}$ is the elastic stiffness matrix of the rock mass; $\text{ }\!\!\varepsilon\!\!\text{ }$ is the total strain tensor of the rock mass; and ${{\text{ }\!\!\rho\!\!\text{ }}_{\text{b}}}$ is the bulk density of the rock mass. This theoretical framework forms a closed-loop bidirectional multi-field coupling mechanism, in which temperature elevation in the rock mass induces thermal strain and triggers global stress redistribution; the alteration of effective stress regulates the degree of fracture closure and aperture evolution, thereby modifying fracture permeability and seepage field distribution; the convective heat transfer conditions subsequently vary and inversely modify the temperature field distribution. Thus, the complete evolution process of dynamic multi-physical field interaction in the surrounding rock is fully reproduced.

To overcome the idealized limitation of constant thermophysical parameters in conventional coupled models, the temperature-stress coupled constitutive relationships calibrated through laboratory tests in Chapter 2 were embedded into the thermal-hydraulic-mechanical governing equations, enabling real-time dynamic iterative updating of core thermophysical parameters in accordance with the multi-field state. The heat conduction and thermal storage performance of the rock mass are governed by the synergistic regulation of temperature and stress, with the dynamic evolution equations expressed as:

$\text{ }\!\!\lambda\!\!\text{ }_{\text{s}}^{\text{eff}}\text{(}{{\text{T}}_{\text{s}}}\text{,}{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}\text{)=}{{\text{ }\!\!\lambda\!\!\text{ }}_{\text{s0}}}\cdot {{\text{f}}_{\text{ }\!\!\lambda\!\!\text{ }}}\text{(}{{\text{T}}_{\text{s}}}\text{)}\cdot {{\text{g}}_{\text{ }\!\!\lambda\!\!\text{ }}}\text{(}{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}\text{)}$                (19)

${{\text{(}{{\text{c}}_{\text{p}}}\text{)}}_{\text{s}}}\text{(}{{\text{T}}_{\text{s}}}\text{,}{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}\text{)=(}{{\text{c}}_{\text{p}}}{{\text{)}}_{\text{s0}}}\cdot {{\text{f}}_{\text{c}}}\text{(}{{\text{T}}_{\text{s}}}\text{)}\cdot {{\text{g}}_{\text{c}}}\text{(}{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}\text{)}$                 (20)

where, $\text{ }\!\!\lambda\!\!\text{ }_{\text{s}}^{\text{eff}}$ is the effective thermal conductivity of the rock mass under coupled stress and temperature effects; ${{\text{ }\!\!\lambda\!\!\text{ }}_{\text{s0}}}$ is the initial thermal conductivity of the rock mass under ambient temperature and pressure reference conditions; ${{\text{(}{{\text{c}}_{\text{p}}}\text{)}}_{\text{s0}}}$ is the equivalent specific heat capacity of the rock mass under variable temperature and stress conditions; ${{\text{(}{{\text{c}}_{\text{p}}}\text{)}}_{\text{s}}}$ is the reference specific heat capacity; ${{\text{f}}_{\text{ }\!\!\lambda\!\!\text{ }}}\text{(}{{\text{T}}_{\text{s}}}\text{)}$ and ${{\text{f}}_{\text{c}}}\text{(}{{\text{T}}_{\text{s}}}\text{)}$ are the temperature correction functions for thermal conductivity and specific heat capacity, respectively; ${{\text{g}}_{\text{ }\!\!\lambda\!\!\text{ }}}\text{(}{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}\text{)}$ and ${{\text{g}}_{\text{c}}}\text{(}{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}\text{)}$ are the mean stress correction functions for thermal conductivity and specific heat capacity, respectively; ${{\text{T}}_{\text{s}}}$ is the temperature of the rock skeleton; and ${{\text{ }\!\!\sigma\!\!\text{ }}_{\text{m}}}$ is the mean confining pressure of the rock mass.

Fracture permeability characteristics simultaneously couple stress-induced fracture deformation and temperature-induced fluid viscosity variation, with a multi-dimensional coupled permeability model constructed as:

${{\text{k}}_{\text{f}}}\text{=}{{\text{k}}_{\text{f0}}}{{\left( \frac{\text{b}}{{{\text{b}}_{\text{0}}}} \right)}^{\text{3}}}\frac{\text{ }\!\!\mu\!\!\text{ (}{{\text{T}}_{\text{0}}}\text{)}}{\text{ }\!\!\mu\!\!\text{ (}{{\text{T}}_{\text{f}}}\text{)}}$                       (21)

where, ${{\text{k}}_{\text{f}}}$ is the real-time fracture permeability; ${{\text{k}}_{\text{f0}}}$ is the initial reference permeability; $\text{b}$ is the real-time fracture aperture under coupled stress and temperature effects; ${{\text{b}}_{\text{0}}}$ is the initial fracture aperture; and $\text{ }\!\!\mu\!\!\text{ (}{{\text{T}}_{\text{0}}}\text{)}$ and $\text{ }\!\!\mu\!\!\text{ (}{{\text{T}}_{\text{f}}}\text{)}$ are the hydrodynamic viscosities of fracture water corresponding to the reference temperature and the real-time fluid temperature, respectively.

This set of equations achieves deep coupling among thermophysical properties, seepage properties, temperature field, stress field, and fracture structure, accurately reproducing the nonlinear evolution of physical parameters of surrounding rock under deep high-geothermal and high-stress environments, and completely eliminating the systematic computational errors introduced by constant parameter assumptions, thereby ensuring the fidelity and accuracy of multi-field coupled evolution simulations.

In the present model, rock displacement, pore pressure, fracture pressure, skeleton temperature, and fluid temperature were selected as the five independent fundamental solution variables. A complete solution system was constituted by three scalar stress equilibrium equations, two seepage governing equations, and two solid-liquid energy conservation equations, with the number of equations strictly matching the degrees of freedom of the unknowns, achieving full mathematical closure of the model. To conform to actual engineering ventilation and heat exchange conditions, a third-type Robin boundary condition was adopted at the tunnel wall, with the Dittus-Boelter correlation employed for dynamic calculation of the wall convective heat transfer coefficient, replacing the simplified treatment of a constant heat transfer coefficient:

$\text{h=(}{{\text{ }\!\!\lambda\!\!\text{ }}_{\text{air}}}\text{/D)}\cdot \text{0}\text{.023Re}_{\text{air}}^{\text{0}\text{.8}}\text{Pr}_{\text{air}}^{\text{0}\text{.4}}$                  (22)

where, ${{\text{ }\!\!\lambda\!\!\text{ }}_{\text{air}}}$ is the thermal conductivity of air; D is the equivalent hydraulic diameter of the tunnel; and $\text{R}{{\text{e}}_{\text{air}}}$ and $\text{P}{{\text{r}}_{\text{air}}}$ are the Reynolds number and Prandtl number of air flow, respectively. This dynamic boundary condition enables accurate characterization of the regulatory mechanism of ventilation velocity on heat exchange efficiency of the surrounding rock, establishing a quantitative relationship between the thermodynamic response characteristics of the surrounding rock and engineering ventilation parameters. The far-field boundaries of the model were set with fixed temperature and stress constraints, with the computational boundary located at a distance from the tunnel center of no less than ten times the tunnel diameter, thereby minimizing the interference of boundary truncation effects on the computed multi-field coupled responses in the near-field surrounding rock.

To elucidate the nonlinear regulatory mechanism of seepage field variations in surrounding rock, induced by excavation unloading of underground tunnels, on the redistribution of the pristine high-geothermal field, and to accurately assess the potential impact of heat and mass interactions in porous media on the safe serviceability of tunnel support structures, a spatial evolution numerical comparative analysis was conducted in Figure 7, focusing on the competitive equilibrium between heat conduction in deep rock masses and convective heat exchange with groundwater. Figure 7 clearly presents the spatial morphological characteristics of the surrounding rock temperature field under two differentiated hydrogeological conditions. Under the low-seepage static heat transfer condition, where pore water migration is minimal and the thermal conductivity remains quasi-isotropic, the isotherms of the surrounding rock strictly maintain a circularly symmetric distribution conforming to the tunnel geometry, indicating that when the seepage velocity falls below the critical threshold, near-field heat migration is predominantly governed by solid-phase conduction through the rock mass, with convective heat transfer effects being negligible. However, once the seepage velocity is elevated to the convection-dominated threshold representing high-seepage conditions, the initial thermal equilibrium state of the surrounding rock is thoroughly disrupted by the vigorous transport of fluid particles, and the high-temperature zone rapidly exhibits a pronounced directional shift along the groundwater seepage vector, resulting in isotherms that display a significantly asymmetric water-drop-shaped or elliptical non-equilibrium distribution at the macroscopic scale. The quantitative results from this sub-figure demonstrate that high-velocity groundwater migration, serving as the core controlling factor in coupled heat and mass transfer in deep tunnels, completely reshapes the non-steady-state thermodynamic response pathways in the near-field surrounding rock through the induction of intense convective heat exchange. This key finding provides compelling evidence that conventional linear thermal field analysis models, which neglect the dynamic feedback of seepage, are incapable of accurately reproducing the true state of surrounding rock in high-geothermal tunnels, thereby furnishing indisputable physical evidence and theoretical support for the asymmetric support structure optimization design and long-term operational safety control strategies proposed in this study, which account for the spatiotemporal evolution of heat and mass processes.

To overcome the technical limitation of the conventional single dry-bulb temperature evaluation, which fails to account for the synergistic thermal hazard effects of temperature, humidity, and thermal radiation, a multi-indicator collaborative evaluation system with rigorous thermodynamic physical foundations was established based on environmental parameters output from the multi-field coupled simulations, enabling refined quantitative assessment of the thermal-humid environment in tunnels. The system incorporates three complementary evaluation indices: wet-bulb globe temperature, heat stress index, and discomfort index. The calculation formulas and physical definitions of each index are provided below. The wet-bulb globe temperature comprehensively characterizes the coupled effects of air temperature, humidity, and environmental radiation, and serves as the core index for engineering thermal hazard classification, with the calculation formula given as:

$\text{WBGT=0}\text{.7}{{\text{T}}_{\text{nwb}}}\text{+0}\text{.2}{{\text{T}}_{\text{g}}}\text{+0}\text{.1}{{\text{T}}_{\text{a}}}$             (23)

where, ${{\text{T}}_{\text{nwb}}}$ is the natural wet-bulb temperature, ${{\text{T}}_{\text{g}}}$ is the globe temperature, and ${{\text{T}}_{\text{a}}}$ is the in-tunnel dry-bulb temperature. The heat stress index primarily reflects the thermal load imposed on the human body by coupled high-temperature and high-humidity environments, and is quantified using the Rothfusz modified formula:

$\text{HI=-42}\text{.379+2}\text{.04901523}{{\text{T}}_{\text{a}}}\text{+10}\text{.14333127RH-}\\ \text{0}\text{.22475541}{{\text{T}}_{\text{a}}}\cdot \text{RH-6}\text{.83783 }\!\!\times\!\!\text{ 1}{{\text{0}}^{\text{-3}}}\text{T}_{\text{a}}^{\text{2}}\text{-5}\text{.481717}\\ \text{ }\!\!\times\!\!\text{ 1}{{\text{0}}^{\text{-2}}}\text{R}{{\text{H}}^{\text{2}}}\text{+1}\text{.22874 }\!\!\times\!\!\text{ 1}{{\text{0}}^{\text{-3}}}\text{T}_{\text{a}}^{\text{2}}\cdot \text{RH+}  \\ \text{8}\text{.5282 }\!\!\times\!\!\text{ 1}{{\text{0}}^{\text{-4}}}{{\text{T}}_{\text{a}}}\cdot \text{R}{{\text{H}}^{\text{2}}}\text{-1}\text{.99 }\!\!\times\!\!\text{ 1}{{\text{0}}^{\text{-6}}}\text{T}_{\text{a}}^{\text{2}}\cdot \text{R}{{\text{H}}^{\text{2}}}$               (24)

where, RH is the ambient relative humidity. The discomfort index is employed for rapid determination of the stuffiness level of the environment, suitable for quick inspection and assessment scenarios in engineering, with the calculation formula given as:

$\text{DI=0}\text{.5}{{\text{T}}_{\text{a}}}\text{+0}\text{.5}{{\text{T}}_{\text{wb}}}$               (25)

where, ${{\text{T}}_{\text{wb}}}$ is the ambient wet-bulb temperature. These three indices form a differentiated and complementary evaluation logic, and through multi-index cross-validation, the one-sidedness of any single evaluation index is effectively avoided, substantially enhancing the accuracy and reliability of thermal hazard risk assessment in high-geothermal tunnels.

Based on the global numerical simulation results and field engineering measurement data, a five-level classification standard for tunnel thermal-humid environments was established in this study, dividing the working environment into five grades: comfortable zone, generally comfortable zone, safe zone, thermal tolerance zone, and hazardous zone. The critical threshold values for temperature, humidity, and thermodynamic indices corresponding to each grade were explicitly defined. The critical threshold for safe engineering operations was set at a dry-bulb temperature of 37 ℃, relative humidity of 64%, and wet-bulb globe temperature of 32 ℃, while the critical threshold for hazardous operations was set at a dry-bulb temperature of 43 ℃, relative humidity of 81%, and wet-bulb globe temperature of 40 ℃. Corresponding grading thresholds for the heat stress index and discomfort index were also established, forming a multi-dimensional closed-loop evaluation system. Based on the combined characteristics of surrounding rock temperature and seepage velocity, differentiated engineering thermal hazard prevention and control strategies were formulated. Under low-temperature and low-seepage conditions, conventional ventilation suffices to meet operational safety requirements; under moderate-temperature and moderate-seepage conditions, increased ventilation rates combined with thermal insulation support structures are required; under extreme high-temperature and high-seepage conditions, ventilation alone is insufficient to mitigate thermal hazard risks, necessitating the simultaneous implementation of seepage sealing and active cooling measures. Fracture seepage significantly deteriorates the thermal-humid working environment of the tunnel, capable of increasing in-tunnel temperature, heat stress index, and wet-bulb globe temperature by 18%, 24%, and 12%, respectively. Moreover, the response sensitivities of different thermodynamic indices to extreme thermal hazard conditions exhibit notable differences, with the wet-bulb globe temperature showing the most sensitive response and dry-bulb temperature exhibiting significant hysteresis, further validating the engineering applicability and superiority of the multi-index joint evaluation system. This system provides theoretical support and a technical basis for graded thermal hazard management and long-term operational maintenance of high-geothermal tunnels.