## 2 Temperature Field Coupling Theory (Airflow to Lining/Heat-Insulating Layer to Surrounding Rock)
### Analytical Theory for Temperature Fields of Tunnel Lining and Surrounding Rock  
The governing differential equation for the temperature fields of the solid domain (tunnel lining/heat-insulating layer and surrounding rock) is expressed as:  
\[C_{eq}\frac{\partial T_{s}}{\partial t}+\
abla\left[-\lambda_{eq}\
abla T_{s }\right]=Q, \tag{7}\]  
where \(C_{eq}\), \(T_{s}\), \(\lambda_{eq}\) and \(Q\) are the equivalent heat capacity, temperature, equivalent thermal conductivity and internal heat source for the tunnel lining/heat-insulating layer or the surrounding rock, respectively. At present, the enthalpy method and the apparent heat capacity method are widely used in finite element analyses of phase change problems in rocks and soil. Both of these two methods are able to provide an accurate prediction for the phase change process. In the enthalpy method, it is usually required to carry out an experiment to obtain an accurate enthalpy function. In the current research, the experiment was not carried out, which means that the accurate enthalpy function is not available. Thus, it is more convenient to use the apparent heat capacity method instead of the enthalpy method. Besides, the accuracy of the apparent heat capacity method in analyzing phase change problems in surrounding rocks and soil has been verified by past studies [9; 10; 19]. However, it should be noted that the equivalent heat capacity increases rapidly in the apparent heat capacity method, which may lead to difficulties and errors in the finite element analysis of phase change problems, especially when the phase transition range is small. The apparent molar heat capacity method proposed by [PERSON] et al. [22] was employed to simulate phase change in the solid domain (tunnel lining/heat-insulating layer and surrounding rock) in a tunnel. Furthermore, the heat capacity of frozen zone (\(C_{\mu}\)), the heat capacity of unfrozen zone (\(C_{u}\)), the thermal conductivity in frozen zone (\(\lambda_{\mu}\)), and the thermal conductivity in unfrozen zone (\(\lambda_{u}\)) are assumed to not depend on the temperature of different domains (\(T\)), while the phase change process is usually assumed to occur within a certain temperature range, (\(T_{m}\)\(\pm\)\(\Delta T\)). The expressions of the equivalent heat capacity and the equivalent thermal conductivity for tunnel linings/heat-insulating layers and surrounding rocks are given as:  
\[C_{\alpha}=\left\{ \begin{array}{c}C_{f}\ T<T_{m}-\Delta T\\ \frac{L}{2\Delta T}+\frac{C_{f}+C_{u}}{2}\ T_{m}-\Delta T\ \leq T\leq T_{m}+\Delta T\\ C_{u}\ T>T_{m}+\Delta T\end{array} \right., \tag{8}\]  
\[\lambda_{eq}=\left\{ \begin{array}{c}\lambda_{f}\ T<T_{m}-\Delta T\\ \lambda_{f}+\frac{\lambda_{u}-\lambda_{f}}{2\Delta T}[T-(T_{m}-\Delta T)]\ T_{m}- \Delta T\ \leq T\leq T_{m}+\Delta T\\ \lambda_{u}\ T>T_{m}+\Delta T\end{array} \right., \tag{9}\]  
where \(C_{u}\) and \(C_{f}\) are the respective heat capacities of the frozen and unfrozen zone; \(\lambda_{f}\) and \(\lambda_{u}\) are the respective thermal conductivities of the frozen and unfrozen zone, and \(L\) is the latent heat of phase change of the water-bearing surrounding rock. Based on the above discussions, it is found that the convective heat transfer coefficient is the key factor to analyze the coupling heat transfer between different physical domains. It should be known that the computed convective heat transfer coefficient is not only related to flow parameters, but also affected by mesh resolutions. Thus, the employed numerical model requires a proper meshing scheme to ensure an appropriate grid size in the near-wall regions. Firstly, the heat transfer coefficient (\(h\)) is calculated based on the turbulent kinetic energy (\(k\)), the dimensionless temperature (\(T^{*}\)), and the normalized distance from the first point of the fluid domain to the wall surface (\(y_{p}\)\({}^{+}\)).