The manuscript presents a massive dataset of δ²H values of inclusion-hosted water for five speleothems, whose age segments are overlapping. The stable isotope composition records of the speleothems fit each other, proving that the signals are reproducible. The speleothems collectively cover the last 16.5 ky, their age precisions are very good. The paper is written clearly in most parts. All these together would support publication in a well known journal like Climate of the Past. The datasets, their descriptions, and the comparisons with independent records are all fine. However, the reviewer misses the discussion of consequences, the explanations how and why these data modify the knowledge of late Pleistocene and Holocene climate conditions and their governing processes in Iberia and the wider region. Both the Abstract and the Conclusions sections are confined to the descriptions of analytical results, the Discussions contain comparisons with other paleoclimate records, but not more. What is the significance of temperature differences between GS and GI periods and the Holocene? I suggest a major revision to address the governing processes of climate changes, most probably initiated in the North Atlantic.

Specific comments:

Fig. 4. B axis titles are mixed.

line 113: it is not clear how the isotopically equilibrium precipitation of carbonate would affect the stable hydrogen isotope composition of inclusion-hosted water. Please explain.

line 207: the method description is too weak. At least the TC/EA technique and the instrument should be mentioned.

line 284: these age periods are mentioned first and last time here.

line 292: precipitation temperature is one of the most important factors of δ¹⁸O values of carbonates.

line 301: plotting the record of seawater oxygen isotope composition would be informative.

line 306: there are too many pieces of information in earlier publications. Rainfall effect is mentioned here without detailed description why the given location is affected.

lines 317-346: I miss the δ²H/T gradient numbers. It is not clear how the δ²H/T relationship was obtained. At line 337 the δ²H/T gradient seems to obtained from a two-point linear regression using MAAT ad δ²H of drip water at two sites. On the other hand, at line 370, it is written that the δ²H/T gradient is obtained by adjusting the δ¹⁸O/T gradient with a factor of 8. Somehow the entire description of gradient calculation is confusing. This should be carefully revised in order to make the process more clear.

line 364: it would be easier to follow if the term δ²Hd was defined here as the δ²H value of drip water.

I suggest to mention here that the equation’s form expresses the fact that the pre-Holocene temperatures are lower than today’s, hence the calculated temperature difference should be subtracted from T(modern).

line 372: line 105 may suggest that the δ²H/T gradient was obtained using direct monitoring data.  Since both the δ¹⁸O and δ²H values monitored along with surface temperature change, this relationship can be directly calculated and it is not necessary to multiply the δ¹⁸O/T gradient by a factor of 8.

line 388: the negative excursion is shown by a single point, it might derive from stochastic scatter. It can be safely written that GS-2.1a had generally cooler than GS-1.

This manuscript provides a rather interesting fluid inclusion isotope record of two caves in the western Pyrenees, that cover the time interval from Heinrich-1 to today. The isotope record is converted to a paleo-T record, based on the present-day relation of rainfall oxygen isotope data with temperature. The resulting temperatures show a pattern of lowest T in H-1, higher T’s in GI-1, lower T’s again in the Younger Dryas (GS-1), and near modern temperatures for the Holocene. This is an interesting high-resolution record, providing valuable information for the reconstruction of paleotemperatures through the deglaciation in Northern Spain

I still have some comments and suggestions, however, that may improve this ms.

• FI isotope analysis:

The technique used for the FI isotope analysis is well-established, and an extensive dataset has been produced for this ms. This has been a formidable achievement, and provides valuable insights. The analytical uncertainty as presented, is based on the replicates, and results in a ~ 2.7 permille precision at 1SD level. Routinely analysed replicates of a standard also suggest the uncertainty is around 2.7 permille (I assume that is a 1SD statistic, but please add that to line 210).  

I believe it would improve the ms to further report the d18Ofi data that were presumably collected alongside the d2Hfi data. I understand and accept that the technique used does not yield top results for d18Ofi (particularly at lower yields), but since you probably have the data, I believe it is better to show that, than to say that. Further, if part of the d18O data are robust, then d2Hfi vs d18Ofi cross-plots can help check the robustness of the analysis as a whole, and help identify analytical artefacts, like those described by Fernandez et al. (2023). I would not ask to do that in the main text, but availability of such data (and plots?) in the supplementary materials would really bolster the quality of your dataset. My stance here would be that one should not write off the d18Ofi data without trying. Quite a few previous publications show that d18Ofi data often preserve rather well in speleothem fluid inclusions, and yield meaningful information.

In relation to that, I would strongly suggest to add a section on (isotope equilibrium) temperature calculations based on d18O calcite - d18Ofi pairs. This is a paleothermometer, used in several previous studies (e.g. Fernandez et al 2023) and has a different (isotope equilibrium) approach to temperature calculation than the one you use in this ms. Even if the d18Ofi data are not good enough, you can recalculate the d2Hfi data to d18Ofi data using the global or local MWL. I expect it may give valuable insights to have two different techniques to calculate paleotemperatures from the same set of data.

• Paleothermometer calibration:

The description on how you calibrated the paleothermometer (with a modern rainfall record further E and higher up in the Pyrenees) is not so clear to me. You use a dataset of individual rain shower datapoints to make a d18O vs T relationship, and then calculate that back to d2H vs T relationship using “the factor eight” which I presume is the slope of the global MWL (or are you using a local MWL?). Several other corrections are applied (ice volume, elevation), which is correct I believe, but makes the entire calculation process a bit difficult to follow for the non-expert reader. Even when there is more info in Giménez et al (2021), I would ask to clarify this in more detail in the text of this paper, because the calculation is essential to your interpretation of the data.

Another specific question is: You calculate via d18O data now (Fig 2B), but why did you not use the d2H data of the rainwater samples, from Giménez et al 2021, straight away (because that saves you an unnecessary conversion from d18O to d2H values)?

• Propagation of uncertainties:

The authors state the uncertainties of the analyses and calibration are fully propagated but I could find no details on how that is done. If I look at the data underlying the d18O vs T relation (Fig 2B), then I see a significant, but seemingly not so strong correlation. The plot suggests a lot of d18O variability is not controlled by T, and I guess that should be very similar if you use the d2H data. You probably need to report the r² for the Villanúa d2H vs T plot to quantify the eventual effect on the uncertainties of the calculated T’s. What I don’t quite understand from the text, in that context, is how you calculated the rather good (0.03 permille) uncertainty of mean annual d18Or and MAAT. 

It would be good to provide more error propagation details to underpin the uncertainties on the eventual T-data that you produce, and the choices you made to get there.

Fernandez et al., 2023: Characterization and Correction of Evaporative Artifacts in Speleothem Fluid Inclusion Isotope Analyses as Applied to a Stalagmite From Borneo (G³)

More specific comments and suggestions

286: Is that really significant? I do agree it is at the “right” point in time, but as a result this does not really stand out for me.

324-329: I don’t quite understand what you have done here. This may need further explanation.

351: you mention Younger Dryas here, but that is not used in your figures (you use GS-1 there)

.

453-454 I don’t immediately understand what a more smoothed T signal is? Less amplitude?

477 the word “record” is once too many in this sentence.

482-483: masking? How does that work? Is that centennial-scale d2H variability that is not T-related?

Line 515: That 2.7 degrees C is awfully close to your T uncertainty. I would certainly not dare to say that that is more pronounced than elsewhere (like you do in line 522) based on the uncertainties you are dealing with.

Table A1: reporting reproducibility with a 1SD metric on only two replicates is statistically a bit on the edge, perhaps.

Fig 4 panel B: Y axis labels are in the wrong place (swapped).