## 2 Materials and Methods
### Temperature Field Governing Equation  
The thermal behaviour of an energy pile is characterized as a transient process. The differential equation that governs heat transfer in the saturated ground is:  
\[\rho c_{\mathrm{P}}\frac{\partial T}{\partial t}-\rho_{w}c_{w}k\
abla\Psi\
abla T =\
abla\bullet(\lambda\
abla T) \tag{1}\]  
Figure 1: Schematic diagram of heat transfer process of energy pile. Low permeable and hydraulic gradients are more favorable to ensuring the seasonal energy balance of the energy pile [3]. Therefore, groundwater-induced thermal convection is neglected, and the heat transfer can be regarded as a pure heat conduction process. Equation (1) can be simplified as follows:  
\[\rho c_{\rm p}\frac{\partial T}{\partial t}=\
abla\bullet(\lambda\
abla T) \tag{2}\]  
In the three-dimensional Cartesian coordinate system and at a constant thermal conductivity, Equation (2) can be transformed into the following Equation (3):  
\[\frac{\partial T}{\partial t}=\frac{\lambda}{\rho c_{\rm p}}\bigg{(}\frac{ \partial^{2}T}{\partial x^{2}}+\frac{\partial^{2}T}{\partial y^{2}}+\frac{ \partial^{2}T}{\partial z^{2}}\bigg{)} \tag{3}\]  
where \(T\) is the transient temperature, \({}^{\circ}\)C; \(t\) is time, s; \(\rho\) is the concrete density, kg/m\({}^{3}\); \(c_{\rm p}\) is the mass-specific heat capacity at constant pressure, J/(kg\({}^{\circ}\)C). The concrete thermal conductivity \(k\), W/(m\(\cdot\)K), can be assumed constant because it is not substantially affected by the temperature variation within the studied range. In order to solve the differential equation above, boundary conditions need to be determined. There are three types of boundary conditions: (1) Dirichlet boundary condition, i.e., the temperature on the boundary is known or is a function of time. It is suitable for the situation if the specific temperature on the periphery of the surrounding ground can be determined after running for a period of time. (2) Neumann boundary condition, i.e., the heat flux on the boundary is known or is a function of time. For the heat convection between the heat exchange liquid and the inner wall of the pipe, it can be simplified (as shown in Equation (4)). (3) Cauchy boundary conditions, i.e., the temperature and heat transfer coefficient of the fluid medium in contact with the boundary, are known, as shown in Equation (5). \[-k\left.\frac{\partial T}{\partial n}\right|_{\Gamma}=q_{2} \tag{4}\]  
where \(q_{2}\) is the heat flux, W/m\({}^{2}\); the outflow boundary \(q_{2}\) positive, which means the direction is from hot to cold; and the outflow boundary \(q_{2}\) is negative. \[-k\left.\frac{\partial T}{\partial n}\right|_{\Gamma}=\alpha(T-T_{t}) \tag{5}\]  
where \(\alpha\) is the heat transfer coefficient, W/(m\({}^{2}\)C); \(T_{t}\) is the fluid medium temperature, \({}^{\circ}\)C.  
The transient temperature field under such complex geometric boundary conditions is unlikely to be solved using the analytical method; the finite element method is an alternative choice [20; 22; 24; 25], which adopts the weighted margin to solve the mathematical problem:  
\[\frac{\partial J^{D}}{\partial T_{l}}=\iint\limits_{D}\bigg{[}k(\frac{ \partial w_{l}}{\partial x}\frac{\partial T}{\partial x}+\frac{\partial w_{l} }{\partial y}\frac{\partial T}{\partial y})+\rho c_{\rm p}w_{l}\frac{\partial T }{\partial t}\bigg{]}dxdy-\oint_{\Gamma}kw_{l}\frac{\partial T}{\partial y}d_{ \rm s}=0 \tag{6}\]  
where \(w_{1}\) is the weighted function, \(n\) is the direction vector of the normal outer line at any boundary surface. The boundary size can be estimated according to the far boundary theory:  
\[R_{\infty}\geq 4\sqrt{\alpha_{\rm s}t} \tag{7}\]  
where \(R_{\infty}\) is the influence radius, m; \(\alpha_{\rm s}\) is the average heat transfer coefficient, W/(m\({}^{2}\)\({}^{\circ}\)C).