We are grateful to the reviewer for constructive comments and suggestions. We hope the following response addresses the reviewer’s concerns outlined in the general comment.

Specific comments

line 48: “atmosphere-alone simulations”. This implies to me that SSTs were fixed, and the usual assumption would be that were fixed for modern conditions. I think it would be important early on in the paper to describe that while the LMDZOR6A simulations are atmosphere only, the SSTs used as boundary conditions (along with other forcing) were ultimately derived from coupled simulations.

We thank the reviewer for pointing out this confusion. These precisions will be added earlier in the paper, in the introduction.

line 69: “Then, the climate described by the global simulations cannot faithfully represent the magnitude of the differences between the different paleoclimate periods, nor microclimate variations at EH2.” There are two issues here: GCM simulations obviously cannot simulate microclimatic variations, but if they can’t “represent the magnitude of the differences between two different paleoclimatic periods” then what good are they? I think climate models passed the test of being able to simulate climates different from present long ago.

Here, by saying that GCM simulations cannot "represent the magnitude of differences between two paleoclimate periods", we refer to the fact that the model cannot represent the exact magnitude at the site, due its resolution, which is slightly different from the second part which refers to the fact that in addition small local scale phenomena are not represented. However, we agree that the model is useful and represent well the large scale regional features as it is discussed in the text in section 4.1.

line 132: “actualist approach” A little more explanation would be good. I was unable to locate Stoetzel et al. (2009), Stoetzel et al. (2014) is basically a review, and that leaves Stoetzel et al. 2011 as the main reference to the approach. But the approach is not really described in depth there, nor do we have a good idea of how well it works in practice.

By "actualistic approach", we mean that in order to calculate the THI, it is assumed that the species' ecological preferences (habitat) were the same in the past as they are today. This assumption will be better explained in the paper. References relate to THI values, which have been retrieved from the literature.

line 295: “There are also important changes in the magnitude of the seasonal temperature variation from June to October and of the seasonal precipitation variation from October to May between these two periods.” If by “these two periods” is meant 115 ka and present, then at some part of the lower temperatures at 115 ka could be related to the calendar effect. At 115 and 116 ka, perihelion occurred in early January and late December respectively, similar to present, but eccentricity was much greater, at its maximum over the last glacial cycle at 115 ka. This led to months in the second half of the year (July – December) beginning and ending farther from the June solstice, and so the calendar effect could have resulted in apparently colder-than-present conditions during that same interval (Bartlein and Shafer, 2019, Fig. 2 and Fig. 11d) that could have reversed the sign of the “real” temperature difference. I don’t think that’s the case here, but it might be good to indicate that the potential role of the calendar effect was considered and discounted.

Potential biases related to calendar effect were indeed thought about. First, the experimental protocol has been designed with care to make sure we do not bias the simulation when applying the correction method to SSTs coming from the different coupled simulations. All simulations have March 21. The way we performed the correction and the daily SST boundary conditions is consistent with the daily SST evolution of the paleoclimate period. Insolation is computed for each longitude, latitude and hour in the day from the orbital configuration of the paleo period.

The second concern with the calendar is the change in the length of seasons (and of month) between the period in the analyses and mainly the statistical analysis. Our concern was to keep the consistency between the different variables in the analyses. Unfortunately, daily values were not saved, so that correcting the monthly values from the calendar effect has to be done variable by variable using the Bartlein and Shafer (2019) method, and not by recomputing the celestial monthly means from daily values. This could introduce some small inconsistencies for variables that should be interdependent as temperature, specific humidity and relative humidity, or evaporation and temperature and humidity that are not linearly related.

This is why, since we are looking at large difference between periods we performed our analyses using the means and standard deviations which are slightly, if at all, affected by the calendar effect, allowing us to bypass it. This requires that event with present day calendar all these values are properly computed by waiting each month by the right number of days in the month. Thanks to this approach there is no calendar effect (or only extremely small on the standard deviation) on the analysis. These comments about the calendar effect and the reason of choosing the mean and standard deviation in the analyses are indeed missing and will be added in the methodological section.

The place where we should have provided figures using variables corrected from the calendar effect are the maps and curves. For the latter the effect is not that much visible because we present the absolute values. For the maps depending on the period it would affect the magnitude, but not the sign for the period and the seasonal averages we present here.

line 331: As described in its legend or in the text, Figure 8 is incomprehensible. How were the climate variables “centered and reduced”?

Variables were normalized between periods, in order to provide a cross correlation analysis. It is these relative variations between periods that are shown. Below the colored tables, maximums and minimums of variables with their original units are shown for information purpose, to indicate which range of value is explored by each variable. This explanation will be clarified in the methodological section of the article, and figure caption will be revised completed.

What was the “common scaling” (z-scores, as implied by line 278)? If the variables are centered, why is a progressive as opposed to diverging color scheme used?

The color scale was chosen so that the intensity of the color reflects the intensity of the variables. This presentation was chosen to made climate variation over the period perceived intuitively: an intense color indicates a maximum while a clear (white) color a minimum. Precision will be added in the figure caption

What observations are the means and standard deviations calculated over?

The means and standard deviations are computed on the annual mean cycle estimated from the last 30 years of each LMDZOR6A simulations, thus they represent respectively the annual mean and the seasonal amplitude of the variables (as indicated in section 2.2.3 and figure captions). This statement will be better explained and further highlighted in the methods section.

Why is the top half of the figure brown, and the bottom blue?

The brown/blue color code is used throughout the article to distinguish between analyses based on DH1 and DH2 (DH for Dating Hypothesis). This was done to easily identify on which DH the presented results are based on from figure to figure. This information will be added to figure captions where needed.

Why do the variable labels alternate in color? What does “Actual” mean?

The color labels of the variables differ simply to make it easier to read, especially on the "m sd" line which would otherwise be confusing. “Actual” refers to current days. It will be replaced by 0k to avoid confusion.

Following the reviewer advice, special care will be given to the description of the figures in the next version of the article, and all these complements will be added in the captions. However part of the confusion comes from the fact that some explanations were missing in the method section and will be added in the text. The figure captions will only better reflect some of the choices in the figure designed to better guide the reader.

line 338: (Top) What do the axes represent, and what are the numbers in parentheses? Why do the points plotted here differ from those on the biplot for DH1?

In Fig 9, we present the two principal component analyses (PCA) performed on the climatic states by layers, based on DH1 (brown) and DH2 (blue), respectively. The axes represent the two leading principal components, along which data are visualized. The numbers in parentheses are the percentage of variance carried by each axis. Most of these methodological details will be expanded and added to figure captions as well. In the biplot of DH1 (presented in SM Fig2) there is indeed an error in the labelling (the point indicated as L6/L7 is in fact L5), we thank the reviewer for pointing it out.

(Bottom) How are “contributions” calculated? Typically, component loadings (correlations between the original variables and components) would be plotted in this context, and ordering the variables by contribution makes it difficult to compare the analyses.

The contributions of original climate variables are computed from the PCA using the fviz_contrib function, using the factoMineR (Lê et al., 2008) and factoextra (Kassambara & Mundt 2020) libraries. Variables are ordered by contribution in order to make easily identifiable which climate variables contribute the most in each case.

Line 385 (and Fig. 10): The motivation for, design, and interpretation of the 2B-PLS analysis needs to be improved.

We will work on the readability of Fig 10 for the next version of the paper. In particular we will make the entire layout more symmetrical for the two dating hypothesis, which seems to be part of the confusion.

What are the two blocks of variables? (It looks like the isotopes and climate PC’s in one analysis, and THI categories and climate PC’s in the other, but why use the components and not the original climate variables?).

Yes, for the 2B-pls the blocks are on the left the isotopes and successively "all climate variables", PC1 and PC2; and on the right the THI categories and successively "all climate variables", PC1 and PC2. The presentation and designations will be changed for clarity. Working on principal components instead of variables alone allowed us to remain in a multivariate framework and thus to take into account the interrelationships between climate variables. This while considering separately the uncorrelated axes of variance, which could behave differently since they are uncorrelated.

Are there really no significant pair-wise correlations between the principal components and environmental variables, except for a few for DH1 isotopes?

As surprising as it is, there is no significant pairwise correlations except between the PC2 of DH1 (climate) and isotopes.

What are the little scatter diagrams?

2B-pls plots are presented when the tests results are statistically significant, so that the reader can get an idea of the dispersion of the points.

What is the basis for the statement that “there is a statistically significant covariation between THI values and all climate variables”? (line 386).

This statement is based on the statistically significant result of the 2B-pls between the THI categories and all climate variables (p-value=0.04*, r-pls=0.9; Fig 10, top, right panel).

I could not find code- or data-availability statements. Several R packages are mentioned in the acknowledgments, but I was unable to find a function for 2B-PLS among them. (There is a function in the geomorph package, but that’s not listed.)

Indeed, for the 2B-pls the two.b.pls function from the geomorph package (Adams & Collyer 2016) was used. The list of the libraries used will be completed. The raw data used in the analysis and a synthetic code will be made available on an online repository according to the request of the journal.

Technical comments

We thank the reviewer for all technical comments that will help bring clarity to the paper. All will be addressed in the next version of the paper. In particular:

• Some precisions will be added in the abstract
• Wording suggestions will be applied

line 196: “They translate…” Does this mean that the biases are propagated into the paleo simulations of the seasonal cycle of surface air temperature and precipitation?

Yes

line 227: Were the SSTs specified as 50-year time series, or long-term means?

The SSTs were specified as long-term means.

line 265: Why would 2B-PLS be preferred to canonical correlation? (e.g. Rohlf and Corti, 2000, Syst. Biol. 49:740-753).

We preferred 2B-pls to canonical correlation analyses because, according to McIntosh (2022, arXiv:2107.06867v2), high correlations within either block can compromise the reliability of canonical correlation results.

line 287: How were these dissimilarities calculated? (What metric/which variables?).

Dissimilarities were assessed on means and sd of the same variables that the ones used in our analyses. We computed the Euclidean distances between EH2 layers and each cell of the area.

I’m surprised that the grid cell that EH2 lies in wasn’t lit up for the midH simulation because the climate is so similar to that at present (Fig. 6). The key for the different shades of blue isn’t specified.

The mean temperature and precipitation are indeed very similar between present days and midH (Fig 6), however other variables display more variation, and there is also differences in seasonality (Fig 8), which might explain why the cell grid containing EH2 is quite light on Fig 7. In the next version of the paper we will provide a key for the blue shades to easy the reading.

lines 343-355: I guess I see as much difference in the shading on Fig. 8 between L2 and L1 as between some of the other major transitions. I see 15 color changes between L3 and L2 and 17 between L2 and L1. Is this a significant difference?

It is indeed a significant difference. However we defined the “major climate transition” based on the inversion of dominant variables (for example on L1 and L2 hydric stress is relatively more important than the precipitation, while in L3 and L4a it is reversed). These transitions are also based on the proximity between layers on the PCA.

We thank the reviewer for his recommendations. Below are our responses the different points that have been raised.

At present the title suggests that you are considering all of the possible uses of paleoclimate simulations. I think the applicability of this research needs to be further developed in the discussion: for example by outlining the conditions under which this kind of approach would work. I would also like a little sentence about the necessity of the performing bespoke simulations for future efforts. I guess you could alter the title to be “An example of…”; and that might mean you could avoid this additional discussion. But I actually think the discussion is quite helpful, so I'd prefer you not to do that.

We can do both. We can adapt the title to best reflect the exact topic of the article, with something like “Refinement of the environmental and chronological context of the archaeological site El Harhoura 2 (Rabat, Morocco) using paleoclimatic simulations.” We will also add some additional development in the discussion to further explore the possibilities of this approach. A fourth section “4.4 Perspectives and limitations of the interdisciplinary approach” will be added to the discussion:

“The interdisciplinary approach between archeology and paleoclimatology presented in this study opens new avenues on how to test the consistency of paleoclimate or paleoenvironment reconstructions in different regions. The contextualization of archaeological and paleontological sites could greatly benefit from this approach. Environmental and/or chronological uncertainties such as the ones encountered at EH2 are unfortunately common in archeology, since the observed differences depend primarily on the methods-specific biases, not especially on the site. However, while the results of this approach are promising in the case of EH2, extending it to other sites is only possible under certain conditions. 1) The concerned site must have been well studied and data on the chronology of the sequence must be available from different methods. 2) Paleoenvironmental inferences must also to be available, and preferably from different sources. 3) The stratigraphic sequence must be composed of sufficient levels to allow the application of statistical methods. 4) Fully coupled climate simulations of the periods of interest must be available, otherwise their complete production would represent a considerable amount of work.

While not all sites meet these criteria, a large number of them have been heavily studied and could benefit from our approach, such as other Moroccan sites (Ben Arous et al., 2020b) or sites from the cradle of humankind (Hanon et al., 2019; Pickering et al., 2019). Moreover, with the development of more powerful statistical tools, this approach could even be extended to other sites with less referenced context. Paleoclimate simulations such as the ensemble we used here are becoming more common and distributed, so that their availability should be less a of a concern in the coming years. The most crucial point is to encourage collaboration between the fields of archeology and paleoclimatology, as expertise in both disciplines is needed to properly combine the different types of information.”

I didn't really get Figure 10, and I was wondering if this could be presented in a more intuitive fashion.

Following this comment, and those of the other reviewer, we will rework the readability of the figures and their captions, especially for figure 10. In particular we will make it more symmetrical between the two datation hypothesis.

I found the EH2 and DH2 acronyms confusing. Is there a way that you can get rid of them?

Acronyms will be modified to ease the reading. For example by replacing “DH2” by “D2” for Dating Hypothesis 2,  and “EH2” by “EH” so that there is less confusion.