## 3 Results and Discussion
### Comparison of Stability of Various Machine Learning Models  
To further compare the performance of different machine learning models, the absolute errors (AE) for each model was also adopted as a statistical metric. The frequency distribution of AE between measured and estimated daily T\({}_{\rm s}\) derived from the three machine learning models under various input combinations during the testing stage was plotted (Figure 10). As shown in the histograms, the AE of the three models mainly distributed at the ranks of \(<\)0.3 \({}^{\circ}\)C, while the proportions of AE greater than 0.5 \({}^{\circ}\)C for the models only ranged from 1.0% to 3.9%, which together suggested that the three types of models all performed well in T\({}_{\rm s}\) estimation. The proportions of AE less than 0.3 \({}^{\circ}\)C differed among model types and input combinations. For the input combination of T\({}_{\rm mean}\) and R\({}_{\rm s}\) (i.e., combination 4), the proportions of M5P and XGBoost models with AE less than 0.3 \({}^{\circ}\)C ranged from 87.4% to 90.9%, higher than the proportions of their corresponding counterparts with the input combination of T\({}_{\rm max}\) and T\({}_{\rm min}\) (ranging from 81.6% to 85.1%). [PERSON] et al. [39] also concluded that the input combination of T\({}_{\rm mean}\) and R\({}_{\rm s}\) provided a more effective estimation of T\({}_{\rm s}\) in the Turkish region. For the complete input combination, the proportions of M5P and XGBoost models with AE less than 0.3 \({}^{\circ}\)C ranged from 87.0% to 93.1%. Moreover, the XGBoost model had the highest proportion of AE less than 0.1 (up to 55.2%) with this input combination among all models, indicating that the XGBoost model had stronger overall estimating performance than both other models. In addition, the XGBoost model outperformed the RF and M5P models in data processing speed, with less computing time required, consistent with results obtained by former research [29, 87, 88].  
Figure 9: Combined indicator values for each model at all stations in the NEC region. (Model 1–18 in the Y axis represents M5P1, RF1, XGBoost1, M5P2, RF2, XGBoost2, M5P3, RF3, XGBoost3, M5P4, RF4, XGBoost4, M5P5, RF5, XGBoost5, M5P6, RF6, and XGBoost6, respectively, which also applies for Figures S8–S13).  
Taking all the previous results together, it can be concluded that, among the different climate zones of China in this study, the best-performed T\({}_{\rm s}\) estimations were observed in the CC region, followed by the SC region. In contrast, models constructed in the study generally showed poor performance in the NEC and IM regions. Models with the combination of T\({}_{\rm mean}\) and R\({}_{\rm s}\) as inputs could archive relatively high performance in estimation, approaching the performance of models in the fully parametric input pattern. Moreover, the former required fewer datasets and was more generalizable. In terms of different types of models, the XGBoost model would outperform the other two models. Moreover, soil organic matter content, soil water content, T\({}_{\rm mean}\), RH, R\({}_{\rm s}\), and U\({}_{2}\) were significantly related to the accuracy of the model estimation of T\({}_{\rm s}\) and were the main factors that caused errors in the model estimation of T\({}_{\rm s}\).