## 3 Numerical Simulation  
### Mathematical Models  
In tunnel construction, the freezing temperature field is a transient heat transmission issue with a water-ice phase change. The main form of heat transfer in DEM is contact heat transfer. The heat conduction law is discretized to derive its difference equation to enable the solution of material-specific heat conduction problems under operating conditions [25]. Neglecting the effect of deformation on temperature, the heat transfer equation is shown in Equation (1):  
\[-\frac{\partial q_{i}}{\partial x_{i}}+q_{v}=pC_{v}\frac{\partial T}{\partial t} \tag{1}\]  
where \(q_{i}\) is the heat flux, \(\mathrm{W/m^{2}}\); \(q_{v}\) is the heat source intensity per unit volume, \(\mathrm{W/m^{3}}\); \(\rho\) is the material density, \(\mathrm{kg/m^{3}}\); \(C_{v}\) is the specific heat capacity, \(\mathrm{J/kg\cdot^{\circ}C}\); \(T\) is the temperature, \({}^{\circ}C\). The relationship between heat flux and temperature gradient satisfies [PERSON]'s law, as shown in Equation (2):  
\[q_{i}=-k_{\mathrm{ij}}\frac{\partial T}{\partial x_{\mathrm{j}}} \tag{2}\]  
Figure 3: Cooling curve of brine. where \(k_{\rm ij}\) is the thermal conductivity tensor, \(W\)/m\(\cdot\)\(\circ\)C.  
According to the sensible heat capacity method [26], the expressions for its equivalent volumetric specific heat C\({}^{*}\) and equivalent thermal conductivity k\({}^{*}\) when simulated are as follows:  
\[\text{C}^{*}= \begin{cases}\begin{array}{ccc}C_{f}&,&T<T_{d}\\ \frac{C_{f}+C_{u}}{2}+\frac{L}{T_{r}-T_{d}},&T_{d}\leq T\leq T_{r}\\ C_{u}&,&T>T_{r}\end{array}\end{cases}  \tag{3}\]  
\[\text{k}^{*}= \begin{cases}\begin{array}{ccc}k_{f}&,&T<T_{d}\\ k_{f}+\frac{k_{u}-\bar{k}_{f}}{T_{r}-T_{d}}(T-T_{d}),&T_{d}\leq T\leq T_{r}\\ k_{u}&,&T>T_{r}\end{array}\end{cases}  \tag{4}\]  
where  
\(C_{f}\), \(C_{u}\)--volumetric specific heat of frost and thawed soils, \(\text{kcal}/\left(\text{m}^{3\cdot\circ}\text{C}\right)\);  
\(k_{f}\), \(k_{u}\)--thermal conductivity of frost and thawed soils, \(\text{kcal}/\left(\text{m}\cdot\text{d}\cdot\text{C}\right)\);  
\(T_{r}\), \(T_{d}\)--melting and freezing temperatures of the soil, \(\circ\)C.  
Ice plays an important role in the frozen soil wall. The DEM uses a parallel bond constitutive model to simulate the cohesive action within ice. It can achieve the transfer of forces and moments between particles when frozen soil is stressed. It can also characterize the damage form of frozen soil [27; 28]. The parallel bond constitutive model corresponds to alienated particles, which are bonded into a single unit by contact bonds between the particles. These contact bonds are similar to springs with constant stiffness on the contact surface, which act together with the sliding model to counteract the forces and moments generated. The specific intrinsic structure model is shown in Figure 4.  
The forces and moments generated by the parallel bond constitutive model contact are shown in Equations (5) and (6):  
\[\overline{F}_{i}=\overline{F}_{i}^{n}+\overline{F}_{i}^{s} \tag{5}\]  
\[\overline{M}_{i}=\overline{M}_{i}^{n}+\overline{M}_{i}^{s} \tag{6}\]  
where \(\overline{F}_{i}^{n}\), \(\overline{F}_{i}^{s}\) denote the normal and tangential components of the force (N); \(\overline{M}_{i}^{n}\), \(\overline{M}_{i}^{s}\) denote the normal and tangential components of the moment (N\(\cdot\)m). Figure 4: Parallel bond constitutive model. Soil freezing and thawing deformation can be regarded as a temperature-dependent transient volume strain problem [29; 30], and the DEM can realize this transient change process of soil by particle radius expansion. DEM thermal coupling is achieved through the existence of temperature differences between particles or between particles and walls, where heat is transferred from one particle or wall to another through contact, causing a change in particle temperature, which in turn causes the particles to expand and generate thermal strain.