# Full Information Numerical Simulation of Two-Dimensional Steady Ground Temperature Field
## 4 Iterative Method for Solving  
Inverse Fourier transform is applied to the temperature field distribution in the space-wavenumber mixed domain obtained by solving Equation (18) to obtain the temperature field distribution in two-dimensional space. Inverse Fourier transform is applied to the anomalous temperature field value in the space-wavenumber domain to obtain the anomalous temperature field in the spatial domain. Since the differential equation satisfied by the temperature field belongs to the elliptic partial differential equation routine [15; 16], so the algorithm in this paper adopts the iterative method to solve the anomalous temperature field in the space domain. By analogy with convergent Born series (compact operator) in electromagnetic field, a stable and convergent compact operator for numerical simulation of ground temperature field is constructed [17; 18; 19; 20; 21]. The iteration format is  
\[T^{(n+1)}=\alpha_{T}T^{(n+1)}+\beta_{T}T^{(n)}. \tag{25}\]  
In the formula, the right end item \(T^{(n+1)}\) and \(T^{(n)}\) represent the total field obtained by the \(n+1\) and \(n\) calculation, respectively \((n\geq 0)\), and the left end item \(T^{(n+1)}\) is the total field of the \(n+1\) iteration updated by the compact operator, \(\alpha_{T}\), \(\beta_{T}\) is the coefficient related to the background thermal conductivity \(\lambda_{p}\) and anomalous thermal conductivity \(\lambda_{s}\), namely  
\[\begin{array}{l}\alpha_{T}=\frac{2\lambda_{p}}{2\lambda_{p}+\lambda_{s}}\\ \beta_{T}=\frac{\lambda_{s}}{2\lambda_{p}+\lambda_{s}}\end{array} \tag{26}\]  
Extensive test results demonstrate that the compact operator achieves stable convergence when applied to the iterative solution of the ground temperature field in non-adiabatic media. The computational results are rigorously validated, confirming the operator's correctness and effectiveness in numerical simulations of geothermal fields. Detailed test cases and convergence analysis are presented in Section 7.  
The specific steps of the proposed algorithm are described as follows:  
Step 1: The model parameters are input, 2-D mesh generation is performed, sampling points in space domain \(N_{x},N_{z}\) are discretized, and sampling points in the wavenumber domain \(k_{x}\) are obtained using the sampling theorem.  
Step 2: Equation (2) is calculated using the finite element method, background temperature field \(T_{0}\) is obtained, and \(T_{0}\) is considered as the initial value of the total temperature field in the first iteration.  
Step 3: Anomalous temperature field in the mixed space-wavenumber domain \(\widetilde{T}_{s}^{n}\) is calculated by solving Equation (3).  
Step 4: Anomalous temperature field in the space domain \(T_{s}^{n}\) is obtained using a 1-D inverse Fourier transform.  
Step 5: A compact operator is applied, and total temperature field \(T^{n+1}\) is updated.  
Step 6: The relative root-mean-square error is determined. If \(RMSE=\sqrt{\frac{\sum\limits_{i}^{N_{k}}\sum\limits_{j}^{N_{k}}\left(T^{(n+1) }\right)_{1-t_{1}}-T^{(n)}\left|\,{}_{t-t_{1}}\right)^{2}}{N_{x}N_{y}N_{z}}}<\varepsilon\), then the results are output. Otherwise, Step 3 is repeated until the accuracy requirement is satisfied. The error in this study is less than 0.1%.  
The algorithm flowchart is shown in Figure 2.