This manuscript investigates the thermal states of permafrost using a recently developed Bayesian inverse method, Ensemble Kalman Sampling (EKS). The study uses borehole and meteorological data from four sites in different regions with continuous permafrost and performs Bayesian trend analysis and inversion to investigate (1) the long-term trend from 2000 to 2020 of air and soil temperature and their relationship using available meteorological and borehole data; (2) modeled ground temperature and its annual variability compared with borehole data; (3) the relationship between change in latent heat (and thaw depth) and change in temperature; (4) the modeled trend of energy change partitioned into latent heat (permafrost), latent heat (active layer) and sensible heat. The results show the two colder sites (Samoylov and Barrow) have substantially different thermal states from the two warmer sites (Bayelva and Parson’s Lake), which are attributed to different historical climatology and soil types determining the contribution of latent heat. The authors highlight the importance of considering the latent heat (and therefore the soil characteristics) and climatology when interpreting the warming trends of permafrost. They conclude the thawing status of areas with warmer permafrost or higher silt/clay content soils has larger uncertainties in response to the ground temperature trend compared to areas with colder permafrost or low silt/clay content soils.

The research topic falls in the interests of the Cryosphere Journal. The manuscript is well-written and presented. However, the manuscript has three substantial issues that need to be addressed. 

(1) The research design has some issues making it unclear if the main conclusions are attributed to the physical processes or the modeling uncertainties. 

First, the available depths of borehole data are not the same. At the two colder sites (Samoylov and Barrow, Fig 4b and 4c), both sites have deep borehole data although Barrow does not have shallow borehole data. In contrast, neither of the two warmer sites (Fig 4d and 4e) has deep borehole data. This could be the main reason causing the much larger temperature variability (Fig 4d and 4e), more scattered relationships in Fig 5, and more observed uncertainties in Fig 6. Therefore, the majority of conclusions made by comparing colder and warmer sites are not convincing. One or more warm sites with deep borehole data are needed to validate this study’s conclusions. It is also worth performing the inverse modeling again on the Smoylov site excluding its deep borehole data to see if its thermal behavior stays the same or changes toward the warm sites. Line 374 seems to demonstrate depth alone cannot explain the variability. However, the statement is not strong because 82 cm is too small on a 10 m scale. Also, the observations of Bayelva also have less variability than those of Parson’s Lake, which likely explains the less variability in the modeled temperature at Bayelva. The authors do have a full section 5.6 to discuss the limitations. While these limitations are definitely important, the current research design is not strong enough to support the conclusions even neglecting other uncertainties. 

Secondly, section 5.3 discusses the role of surface conditions on ground warming based on the n-factor change before and after 2005. Again, uncertainties can be the main driver because no borehole observations are available to constrain the model before 2005. This is another key point made based on the comparison of two data not having the same conditions.

(2) The manuscript has a large space describing the modeling method but most contents are too technical and not accessible to people who are in the cryosphere community but do not have expertise in numerical modeling, inversion, and Bayesian method, etc. The authors focus too much on the advanced topics of the method but completely missed the information on the basic idea of the applied method. Also, in many cases, the authors only cite some references without explicitly describing the methods, which makes the readers difficult to follow or understand. Below, I list the locations where I think additional explanations and clarifications are necessary.

Section 3.1. The introduction of Bayesian inference involves too many technical terminologies. Please consider adding supporting sentences to make it easier for people not familiar with the Bayesian method to understand it.

Lines 140-144. Need to briefly explain the bias correction procedure.

Line 150. Need to briefly explain the numerical procedures and parameterizations of CryoGrid.

Section 3.6 This section introduces a key methodology EKS. It presents the advantages of EKS over MCMC and EKI without explaining the basic theory/idea of EKS in the first place. Again, this makes researchers not familiar with EKS very difficult to follow and understand it.

Line 250. Need briefly explain what a mean vector from Garbuno-Inigo et al. 2020 is. 

(3) The discussion needs more quantitative and specific analysis. When interpreting the results, the authors only briefly propose possible factors without explaining how would these factors impact the results. Below are some examples.

Paragraph 255. I may miss something but I did not get the purpose of this paragraph. It states that the prior distribution over model parameters is important but does not explain what was done to improve performance.

Line 356. This sentence does not explain why Samoylov has deep soil temperature warming faster than the air temperature. Factors other than air temperature should be included here.

Paragraph 360. Besides only presenting the potential factors impacting the soil thermal states, I would include how they impact the thermal states. For example, how does the ground temperature change with air temperature giving increasing (or decreasing) snow thickness and soil thermal diffusivity?

Line 390. Please explain more about why latent heat is lost so that the temperature is warmer. Please also explain why drainage and evapotranspiration cause latent heat loss.

Section 5.6. It would be helpful if include some discussion about the expected changes after addressing each limitation.

Minor comments:

Line 345 information about site location is needed for Biskaborn et al. 2019. 

Line 386. This may be due ‘to’ the thermal…

Line 445. The second point. Warm permafrost could have slow refreezing when warming due to the effects of latent heat. 

Fig B1. In my opinion, this figure is important as it shows the basic settings of forward modeling. Please consider moving it into the main text.

Line 170. The depth information of each layer is missing.

Dear authors,

I read with interest your paper entitled “Investigating the thermal state of permafrost with Bayesian inverse modeling of heat transfer” which proposes a stochastic inversion of well temperature data to investigate the current trend in permafrost evolution. I don’t have the necessary expertise in permafrost to comment in depth the discussion and interpretation of the results, so I focused my review on the technical aspects linked to Bayesian inversion. The study uses a state-of-the-art technique (EKS) that uses an iterative linearization of the inverse problem (based on the approximation of covariance matrices) to converge towards an ensemble of models of the posterior distribution. Since the method is already published, and that no additional development is made, the description is kept to a minimum (which is fine), but the analysis seems robust and scientifically sound. I regret the absence of a more through sensitivity analysis for the parameters that were not included. I have some comments that could be included in the manuscript to further strengthen the message:

1. L65-73. I regret that other solutions to this issue are not discussed at all. If EKS is one possible avenue, others have been proposed, such as the combination of efficient multiple-chain McMC algorithm with reduced dimension representation of the parameter space (Laloy et al., 2018) , or bypassing the inverse problem by directly predicting the posterior distribution from simulation-based machine learning approaches (Thibaut et al., 2022 for example in hydrogeology dealing with the same type of prediction (temperature field)). Since this aspect is also included in the discussion, it could be interesting to expand the perspectives beyond the technique used in the paper (i.e. EKS).
2. L133-134. The prior does not encode information about Y, it encodes information known about the unobserved quantities (X), before the data Y is actually collected and thus correspond to what we know and don't know about X before the experience.
3. Section 3.2. It was not directly clear to me that a Bayesian approach was applied to the trend analysis. Maybe it could be more explicit.
4. L134-136. Is it ? Since this integral is a constant for a given problem, the comparison of the likelihood (ratio) is sufficient to sample the posterior distribution (see rejection sampling or Metropolis sampling) and the integral is not such a problem. The main problem lies in the computation of the likelihood p(Y|X) which requires to solve the forward problem, generally through numerical approximation of partial derivative equations.
5. L190-191. I am wondering about the effect of this fixed boundary. Given the effort for modelling the uncertainty on the upper BC, why not also considering uncertainty on the bottom one ? This flux is certainly not known for sure and it could contribute to significant uncertainty at depth.
6. Section 3.4. is empty.
7. L218-220. This is a weird formulation. Any Bayesian approach will include some prior uncertainty on model parameters, and if modelling error is often neglected, it is generally included in the observation error from the likelihood. Maybe what is specific to your approach is that the target X and the predictor Y are actually the same (temperature) and that you have first to estimate the distribution of parameter phi?
8. L224-225. It sounds like a classical McMC approach wouldn't work. Any method sampling the posterior can solve the problem, right?
9. Computational time for one forward model?
10. L253-254. Assuming uncorrelated noise in time and space might be one of the unrealistic assumptions, depending of the type of sensors of course. Maybe mention it in the discussion ?
11. L258-259. It is also a requirement for any Bayesian inference. The posterior is directly related to the prior, so the prior should reflect the actual knowledge about the site. McMC is sampling from the prior distribution as well.
12. L284-291. I wonder about the validity of deducing trend with such short data sets. In the discussion, comparison with longer trend is introduced, but I feel this could be strenghtened.
13. What is the reason ? Is this biased visible consistenly throughout the year. If yes, all the models of the posterior must have a high misfit, what could indicate a lack of consistency between the prior and the data set (e.g., Lopez-Alvis et al., 2019).
14. Close to what? I am not sure the sentence makes sense. Please clarify.
15. L374-375. Is this significant enough to state that the difference is not only due to the depth of the sensor ? If both had some sensors deeper, this would largely reduce the uncertainty and has likely nothing to do with the fact that they are warmer, don’t you think ?
16. Why a reference to the median suddenly ? The discrepancy is for all models, not only the median or the mean, all predictions are wrong.
17. L386-387. Have you tried enlarging the prior ? Do you observe a similar bias for other years ?
18. What do you mean by “smooth regularizers”. One of the advantage of Bayesian inversion is to use realistic prior and avoid regularization (such as smoothing). Smoothing is normally only visible in the mean which is not necessarily a sample of the posterior.
19. I think the term “appropriate” is badly chosen. The prior should reflect the uncertainty on model parameters before considering the data set, and what is mentioned should then result in larger prior. The main challenge is maybe to consider correlation between parameters ?
20. L515-516. This is also the case for most McMC approaches.
21. L520-523. See comment 1.
22. Table B5. Display prior distributions behind posterior in figures A2 to A5 to immediately grasp the reduction in parameter uncertainty ?

Laloy, E., Hérault, R., Jacques, D., & Linde, N. (2018). Training-Image Based Geostatistical Inversion Using a Spatial Generative Adversarial Neural Network. Water Resources Research, 54(1), 381–406. https://doi.org/10.1002/2017WR022148

Lopez-Alvis, J., Hermans, T., & Nguyen, F. (2019). A cross-validation framework to extract data features for reducing structural uncertainty in subsurface heterogeneity. Advances in Water Resources, 133, 103427.

Thibaut, R., Compaire, N., Lesparre, N., Ramgraber, M., Laloy, E., and Hermans, T. Comparing Well and Geophysical Data for Temperature Monitoring within a Bayesian Experimental Design Framework. Water Resour Res, 58, e2022WR033045. https://doi.org/10.1029/2022WR033045