Summary: Past global surface temperatures over the past few centuries can be estimated from present borehole temperature profiles applying inversion techniques based on the solution of the heat transfer equation. Form a set of sites where temperature profiles have been measured, large-scale temperature reconstructions can be derived by averaging the local retrievals. The manuscript presents a method to estimate the uncertainties  in the large-scale average

 estimations based on a bootstrap approach. the authors  conclude that this new method provides better and more realistic uncertainty estimates that previous methods. Those previous methods simply calculated the average of the high and low ends of the local uncertainty ranges.

Recommendation: I think that in general the manuscript is valuable  and should be published after some revisions. However, I am afraid that one of the motivations of the present study, namely that the previous estimations of uncertainness was unrealistic, contains a conceptual misconception, although it has been previously published. Therefore, the motivation of the present manuscript should be amended to present a correct statistical case. I explain below in more detail my main  concern.

1) The manuscript presents a base method to estimate global or large-scale uncertainties that has been published previously. This method just constructs the high-end (and low-end) uncertain range of the global average by calculating the average of the high-end (or low-end) range of the local estimations. This is, however, not correct , as it can be illustrated in a short counter-example. The interpretation of a 5-95% uncertain range in a frequentist approach is that the range covers the true value with 90% probability (technically, it means that a putative infinite number of realizations of the measurements and their corresponding uncertain estimations will contain the true value 90% of the time). For the sake of this reasoning , we can a bit sloppily say that that the probability that the true value is within the estimated uncertainty range is 90%. However, if the uncertainty ranges are constructed by averaging the  5% and the 95% local ranges, this probability is much much larger than 90%. Let us focus on the high end (95%). The probability for the average  to be outside that 95%  range is not 0.05, but  actually 0.05 to the Nth power , where N is the number of profiles (sites).This results because each profile from which that average is constructed, has a probabibily 0f 0.05.  If N=100, this number is very small,, much smaller than 0.05. 

The authors realize in the discussion that indeed this estimation is not correct i. There, they apply a much more correct estimation assuming that the global profile is the average of N random variables, and therefore, assuming that these N random variables are independent, the error  in the average amounts to the sqrt of the average squared error. If all individual errors are equal, this amounts to that individual error divided by the sqrt(N).

There is one  important underlying assumption:  the errors should be independent across space. But even if this assuption is not completely fulfilled, this estimation is much more realistic that simply the average of the upper and lower local percentiles, which is clearly incorrect.

Thus, to some extent, the manuscript corrects a previous statistical  misconception. In this sense, it is useful, but the motivation of the manuscript should be cast differently,as the reader will be really surprised to see, without any caveat,  a clearly wrong method as a benchmark.

On the other hand, the question of the spatial correlation of uncertainties, which is critical for the validity of both methods (bootstrap and error propagation) is not mentioned at all.

The bootstrap approach is definitively better - and I could not see any clear error in this application of bootstrapping. However, this approach does not take into account the possible spatial correlation of the local errors. I do not know how significant these correlations might be, but if they are, then the bootstrap estimation of the uncertainty will be too narrow - in the same as the error propagation proposed by the authors in the discussion - since the effective number of degrees of freedom will not be N, but smaller. 

If these correlations are relevant, the bootstrap should take it into account, e.g. by block-bootstrap, in which correlated regions are first arevaged together, and then bootstrapped. I think that this problem is technically very difficult to solve satisfactorily, but again, I believe that the presented bootstrap approach is indeed useful.

Particular points

2) A definition of the quasi-equilibrium temperature will help some readers

3). Section 3.5, perhaps the most important section,  is not very clearly written (I needed to read it several times). For instance , line 257: ‘ named Sampling and Bootstrapping ensembles (S and B ensembles in Figure 2). The Sampling ensemble consists ...’. and the reader expect the following sentence to explain what the Bootstrapping ensemble is. However, the text goes on with ‘  The BTI method considers the uncertainty arising from ...’. This is quite confusing. Actually, the bootstrapping ensemble is a typical bootstrap sampling from the set of individual local profiles, where each profile has been derived from one value of the uncertain parameters (T0, Gamma0, and thermal conductivity). The only restriction is that each sites contributes with one member to the ensemble. 

All in all, I found the technical description unnecessarily too cumbersome.

In this article, the authors introduce a new technique to estimate ground surface temperature and ground heat flux histories.  Using an artificial temperature profile and real data from the Xibalbá subsurface temperature profiles, they show how this new technique performs against the tested singular value decomposition and perturbed parameter inversions methods.  Little difference is noted when examining individual profiles.  However, notable differences in the uncertainty range for the inversion of large sets of subsurface temperature profiles are observed with the new bootstrap technique producing smaller uncertainty ranges.  

The authors clearly explain the importance of these various methods and shedding light on the thermal state of the land surface.  They elaborate the three different methodologies and provide useful figures to illustrate their points.  I found Figure 2 to be especially helpful.  I believe that the outlined new methodology is a useful new tool for the community and would recommend publication after some minor changes.

Comments:

There is little discussion in the data section as to the uncertainty associated with the proxy-based temperatures and the Xibalbá subsurface temperature profiles.  Are there errors associated with converting the PAGES2k global temperatures to a land temperatures?  With respect to the Xibalbá subsurface temperature profiles, they haven’t all been measured at the same time.  This will influence the reconstruction.  I think a couple of lines here elaborating these uncertainties could be helpful.  It would also help the reader further understand the complexity of undertaking these inversions as the data is not perfect and that any tool that can minimize uncertainty is important.


On L153 and later in Section 3.8, the authors state how the continental subsurface is considered homogeneous (i.e. thermal properties are constant) and give the example of the Arctic.  However, I cannot believe that the Arctic could be considered representative of the globe.  How can one consider a constant thermal diffusivity and conductivity for all regions?  Is this an unfortunate trade off due to the lack of subsurface thermal data?

On L418, the authors state: “Another remarkable results is the agreement between the PAGES2k land temperatures and the ground surface temperature histories for most of the period.”  While this appears to be true for a quick glance at Figure 7, notable differences can be seen when examining Table 2.  The three methods reconstruct a warming about two times greater in 1950-2000, 1900-1950, 1850-1900, and 1800-1850.  The authors do note this following L418.  However, I believe this sentence should be rephrased to emphasize the how the tendency is captured but not the magnitude and/or highlight the excellent job the methods do in reconstructing the period of 1600-1800.

In Figure 8b, PPI and BTI with varying κ and λ show an increase in heat flux as of about 1970 that is not observed in SVD nor PPI and BTI with constant κ and λ.  This is not elaborated in the text.  It would be helpful to a reader to have a couple lines clarifying this.

  

Technical Points:

L51-52: Should read: “Nevertheless, several sources of uncertainty arise in the inversion process, the most important being…”

Overall, watch the use of “consists in” in the text. For the majority of time used, it should actually be “consists of”

L59-60: Should read: “..the deepest part of the observed profile, then providing a best estimate…”

L101: “These experiments also allow to identify…” should read “These experiments allow for the identification of…”

L153: Should read: “…typical vary by a relatively small…”

L196: Should read: “Finally, small eigenvalues are from S-1 in order to stabilize the solution…”

L204-205: Should read: “The errors in the estimates of the long-term surface temperature (T0) and the geothermal gradient (Γ0) have …”

L251: Should read: “As explained in the Introduction, the SVD and PPI techniques do not provide a comprehensive…”

L263: “..which are retrieved form…” should be “which are retrieved from…”

The units of thermal diffusivity are shown as m2s but they should be m2s-1

L394: sensible should be sensitive

L433: “This three factors…” should be “These three factors…”

In Table 1, please define Ns and NB

In Figure 3, the shading associated with the purple line isn’t purple but blue.  For consistency with the other figures, I recommend using the same purple shading colour as the other figures.

In the title of Figure 4, I would clarify that B in the title stands for Bootstrapping.