1. In the introduction (lines 92–104), the authors clearly outline two key challenges in current research: first, the significant heterogeneity of Ts leads to unclear relationships between variables; second, modeling is hindered by data scarcity and uneven distribution. However, in lines 106–113, when introducing the objectives and scope of this study, the authors do not explain how the study addresses these two challenges. It is also unclear what specific methods are used to overcome them, and why these methods are effective. It is recommended that the authors restructure this section by focusing on the core problems, rather than simply listing the research contents. This would improve the clarity and logical flow of the introduction.

1. In Section 2.1, the authors describe the use of CMA Ts observational data. However, it is unclear how these data were processed. Were the observations directly provided as daily averages, or were they aggregated from hourly data? Was any quality control applied? How were missing data handled, both in the vertical profile and in the time series? Were any filtering or screening steps performed, and if so, what were the specific criteria?

1. In lines 186–190, as well as in Section 4.3, the authors provide a brief discussion of the study’s limitations. However, it is concerning that the missing land surface temperature (LST) data caused by cloud cover were filled using a simple linear interpolation method. This approach may be questionable, as the interpolated values represent a theoretical cloud-free state, while cloud presence can significantly influence radiative transfer and thus impact LST. There are existing interpolation methods that take into account energy transfer and energy balance. It is recommended that the authors investigate these alternatives and consider adopting a more reliable method.

1. In Section 2.3.1, it is suggested to provide further explanation of the Variance Inflation Factor (VIF). Specifically, what is its purpose, how is it calculated, and if possible, a formula should be included to make the description more complete.

1. In Section 2.3.2, a substantial portion is devoted to the spatial partitioning strategy based on a rotated quadtree. I have several questions regarding this part. First, why was the quadtree data structure chosen? The manuscript does not clearly explain this. Is it intended to address the issue of uneven distribution of observation sites? If so, why is the quadtree suitable for this purpose? Second, what was achieved by using the quadtree? Was there an effort to ensure that each node contains a roughly equal number of sites, for example around 30? Why was 30 selected as the threshold, and what is the basis for this value? Lastly, a minor suggestion (optional for consideration): if the goal is to achieve a more balanced spatial distribution of stations, a top-down data structure such as the K-D tree (with K = 2 in this study) may be more effective than the bottom-up quadtree. A K-D tree can ensure the difference in the number of points between leaf nodes does not exceed one, and can also support rotation operations.

1. In Section 2.3.3 (lines 285–295), the authors introduce XGBoost as the core machine learning algorithm used in the study. They present its advantages and compare it with other methods such as SVM, RF, and neural networks. However, the stated advantages are not sufficient to demonstrate that XGBoost is superior to the other listed methods. Machine learning models differ in structure, number of parameters, optimization strategy, and suitability for different tasks. Therefore, the current explanation is not enough to justify the model choice. Considering that the algorithm is not the main focus of this paper, it is suggested to either include a brief comparative experiment to support the claimed superiority or rephrase the section to emphasize the strengths of XGBoost without direct comparison to other models.

1. In lines 296–297, the validation set is twice the size of the test set. Is this split reasonable, and can it effectively evaluate the generalization performance of the model? Why not adopt more common ratios such as 8:1:1 or 6:2:2? In addition, the manuscript later mentions that five-fold cross-validation was used for evaluation. In this context, what are the roles of the two validation sets? Are they used for model selection, parameter tuning, or testing? It is recommended that the authors provide a clearer explanation. It is also suggested to report the specific sample sizes for each dataset.

1. The authors produced data at a 1-kilometer resolution for China. How did the authors account for the spatial scale difference between point observations of soil temperature and the 1-kilometer resolution results? How was it ensured that the dataset constructed through point observation training could represent results at the 1-kilometer spatial scale? Additionally, regarding the dataset production, I am very interested in the subsequent maintenance and updates of the dataset over time. Can the authors' method be extended to produce datasets for subsequent years?

1. How did the authors account for the impact of uneven spatial distribution of observation data points on the development of the national soil temperature dataset? How do factors such as topography, landform, and vegetation cover types influence the results and uncertainties of this dataset?

1. There are also some minor issues that should be addressed. For example, in Figures 6 and 7, it is recommended to include a color bar legend. As it stands, it is difficult to interpret the exact values represented by the orange points. In Equation (2), the variables x and y lack subscripts i. In Equation (4), the variable i used for summation is not defined. In the references, line 667 and 763 include “others” among the authors—what does this mean? It is suggested to carefully check the manuscript for such details, including grammar, figures, and reference formatting.

1. The method used many data (primarily including in-situ observations and indicator variables) to produce soil temperature. By the way, the in-situ in Figure 3 is wrongly spelled as in-suit. Since these data are with varying spatial scales, and many complicated steps are involved in this procedure to produce the Ts at different depths. I just wonder why the outputted Ts is with that good accuracy. Given than even the acknowledged MODIS LST (nearly Ts at 0 cm) is 1-2K, and it has been taken as an input in this study.
2. To my knowledge, LST changed very quickly and is seriously affected by cloud. The local observation time differ across China, and most regions in the South are covered by cloud at most time. How you process these data, and whether the accuracy can be guaranteed in your study?
3. On the other hand, the seemingly good accuracy is not that strange. Because the authors used the same ground measurement to validate the estimated values. Although the entire data has been divided into two sections of training and validation. They are actually homologous with the similar schemes by CMA. How about validating the estimated results with data collected from different sources
4. The authors used on XGBoost, why not try other machine learning algorithms. It is not sure that XGBoost perform best. Maybe a balance of multiple algorithms is more convincible.
5. Like many other overabundant pure machine learning articles, the present study lacks of innovation, but eligible as a data description paper. Reanalysis data such as ERA5-Land also have Ts at multiple layers, except for the finer spatial resolution (they can also do that if they want), what are the advantages of your data? Why do you think a user should consider your data?

This study generates a high-resolution (1 km, daily) multi-layer (0~40 cm) soil temperature (Ts) dataset for China covering the period 2010-2020. The authors developed a "spatially adaptive layer-cascading" framework based on the XGBoost machine learning algorithm. The spatial adaptiveness was achieved through a rotated quadtree partitioning scheme to handle the non-uniform distribution of in-situ observation sites, while the layer-cascading approach used the temperature of an upper soil layer to predict the temperature of the layer below. However, I have several major concerns about the robustness of the methods and reliability of the results, as detailed below. Based on these, I cannot recommend its publication, at least not in its current form.

1. The authors state that the sites were randomly split into training (70%), validation (20%), and test (10%) sets. For geospatial data like soil temperature, which exhibits strong spatial autocorrelation, this random splitting is a critical methodological flaw. It almost certainly leads to "data leakage", where test sites are geographically close to training sites. Consequently, the model can achieve high performance on the test set even though it did not learn the true underlying relationships between predictors and Ts. This means the model's ability to generalize to new, un-sampled areas is not being properly evaluated. The reported performance metrics (e.g., R² > 0.93 in Fig. 5) are therefore very likely to be significantly inflated and overly optimistic.

The authors should implement a more rigorous validation scheme that accounts for spatial autocorrelation. A spatial block cross-validation approach is strongly recommended.

2. The manuscript claims that the generated dataset accurately captures the spatial distribution of Ts. However, the evidence provided is the high R² (and low RMSE) of the daily time series at individual stations. These temporal variations are heavily dominated by the seasonal cycle, which is easy for any model to capture using predictors like air temperature. A high temporal R² does not prove that the model correctly reproduces the spatial gradients across China. I suggest the authors conduct a spatial-only validation, using mean Ts (for the whole year and for specific seasons) across the sites.

3. The results indicate that model performance is worse at the surface (0 cm) and improves at intermediate depths (e.g., 5-20 cm), as shown in Figures 4-7. This is a counter-intuitive result given the layer-cascading methodology, where the prediction for a deeper layer depends on the prediction from the layer above. This structure implies that errors from the surface prediction should propagate downwards, theoretically leading to a degradation of performance with depth. This apparent paradox should be discussed. 

4. In the VIF analysis (Fig. S1), sand, silt, and clay percentages were included. As these three variables are compositional and should sum to a constant (100%), they are perfectly collinear by definition. This should result in an infinite (or extremely large) VIF values. However, the reported VIFs are relatively low (5.6 to 10). This discrepancy is concerning and suggests a methodological error. 

5. The model uses solar radiation as a predictor but omits downward longwave radiation (LWD). Considering that LWD is a critical driver of the surface energy balance (particularly for nighttime and winter temperatures) and that LWD has been identified as a main driver of Ts trends in process-based models (Peng et al., 2016, https://doi.org/10.5194/tc-10-179-2016), I suggest the authors include LWD as a predictor, or provide a strong justification for its exclusion.

Minor comments:

Line 490 “Notably, RMSE at the surface (0 cm) is slightly lower than at 40 cm, possibly due to stronger direct influences from surface cover and meteorological conditions.” – This is not the case for Fig. 12 cd. Furthermore, making this statement based on only a few sites is not adequate.

Line 535 “Figure. S5 demonstrates that LST is more effective than air temperature in detecting spatial variations in surface Ts in sparsely vegetated areas” – I do not see how this conclusion can be derived from Fig. S5.