The use of paleoclimatic simulations to refine the environmental and chronological context of archaeological/paleontological sites.

. To reconstruct the paleoenvironmental and chronological context of archaeological/paleontological sites is a key step to understand the evolutionary history of past organisms. Commonly used method to infer paleoenvironments rely on varied proxies such as faunal assemblages and isotopes. However, those proxies often show some inconsistencies. Regarding estimated ages of stratigraphic layers, they can vary depending on the dating method used. In this paper, we tested the potential of paleoclimate simulations to address this issue and contribute to the description of the environmental and chronological context of archaeological/paleontological sites. We produced a set of paleoclimate simulations corresponding to the stratigraphy of a Late-Pleistocene Holocene site, El Harhoura 2 (Morocco), and compared the climatic sequence described by these simulations to environmental inferences made from isotopes and faunal assemblages. Our results showed that in the studied site combined US-ESR ages were much more congruent with paleoenvironmental inferences than OSL ages. In addition, climatic variations were found to be more consistent with isotopic studies than faunal assemblages, allowing us to discuss unresolved discrepancies to date. This study illustrates the strong potential of our approach to refine the paleoenvironmental and chronological context of archaeological and paleontological sites. 2015). Conversely to paleoenvironment proxies, paleoclimate simulations are much influenced by large scale climate processes because they rely on the global dynamic of the Earth’s climate system. 25 We propose to combine those two fields and to use climate simulations to discuss and refine the environmental and chronological context of El Harhoura 2 (EH2) cave, an archeological site located on the north-Atlantic coast of Morocco, whose microvertebrate assemblages have been extensively studied. and L4a. The second principal component, explaining 42.35% of the observed variance, separates the group of L2, L1 and Act from others. L8, L7, L6 and L5 are characterized by a hot and dry climate and a low seasonal variation. L3 and L4a are defined by a wet and windy climate, with important precipitations and high seasonal variability. Finally, L2, L1 and Act 320 present a hot environment associated with an important water stress.


Abstract.
To reconstruct the paleoenvironmental and chronological context of archaeological/paleontological sites is a key step to understand the evolutionary history of past organisms. Commonly used method to infer paleoenvironments rely on varied proxies such as faunal assemblages and isotopes. However, those proxies often show some inconsistencies. Regarding estimated ages of stratigraphic layers, they can vary depending on the dating method used. In this paper, we tested the potential of paleoclimate simulations to address this issue and contribute to the description of the environmental and chronological context of archaeological/paleontological sites. We produced a set of paleoclimate simulations corresponding to the stratigraphy of a Late-Pleistocene Holocene site, El Harhoura 2 (Morocco), and compared the climatic sequence described by these simulations to environmental inferences made from isotopes and faunal assemblages. Our results showed that in the studied site combined US-ESR ages were much more congruent with paleoenvironmental inferences than OSL ages. In addition, climatic variations were found to be more consistent with isotopic studies than faunal assemblages, allowing us to discuss unresolved discrepancies to date. This study illustrates the strong potential of our approach to refine the paleoenvironmental and chronological context of archaeological and paleontological sites.

Introduction
Reconstruct paleoenvironments is a key step to understand the ecological context and evolutionary history of past organisms, and, hence, the link between their phenotype and the environment they lived in. Commonly used method to infer paleoenvironments rely on the organisms themselves, such as the 5 qualitative presence/absence of species (e.g., Tchernov, 2002;Avery, 2007;Denys et al., 2018), the association of certain taxa with their preferred modern habitat (e.g., Fernandez-Jalvo et al., 1998;Stoetzel et al., 2011), their relative abondances (e.g., Peters & Bork, 1999;Matthews, 2000;Belmaker & Hovers, 2011;Comay & Dayan, 2018) or variations in the amount of stable isotopes in organic tissues (Anderson and Arthur, 1983;Tieszen, 1991;Royer et al., 2013). Those approaches allow one to 10 characterize the local environment in which organisms lived. However, drastic discrepancies can be observed between those proxies (e.g., Jeffrey, 2016;Davis & Pineda-Munoz, 2016). Relying on the chronological context can also help to understand global environmental conditions. However, depending on the dating method used, the age of the stratigraphic layers may vary and also lead to discrepancies in the chronoclimatic reconstructions (Ben Arous et al., 2020b).

15
Besides paleoenvironmental reconstructions, paleoclimate modeling is another field for which one of the objectives is to describe past environments. To do so, it simulates paleoclimates using physically based models that describe the dynamics and thermodynamics of the atmosphere, ocean, and continental and frozen surfaces coupled with water and carbon cycles. Those models efficiently describe the directions and large-scale patterns of past climate changes (Braconnot et al., 2012). They characterize 20 past climates through quantified variables and have long been used to investigate past environmental changes (Kutzbach and Otto-Bliesner, 1982;Braconnot et al., 2012;Duplessy and Ramstein, 2013;Schmidt et al., 2014;Harrison et al., 2015). Conversely to paleoenvironment proxies, paleoclimate simulations are much influenced by large scale climate processes because they rely on the global dynamic of the Earth's climate system.

25
We propose to combine those two fields and to use climate simulations to discuss and refine the environmental and chronological context of El Harhoura 2 (EH2) cave, an archeological site located on the north-Atlantic coast of Morocco, whose microvertebrate assemblages have been extensively studied. https://doi.org/10.5194/cp-2021-185 Preprint. Discussion started: 5 January 2022 c Author(s) 2022. CC BY 4.0 License.
What are the paleoclimate changes describe by climate simulations over the sequence at EH2? Are they consistent with paleoenvironmental inferences available in the literature made from faunal assemblages 30 and isotopes? To answer these questions, we produced a set of climate simulations corresponding to EH2 stratigraphic layers. In order to discuss the chronological framework, we assessed the consistency between the climatic sequences describe by simulations according to the varied dating methods and paleoenvironmental proxies from previous studies, such as isotopes and faunal assemblages. We thus expected that the congruence between paleoclimate simulations and paleoenvironmental proxies will 35 allow us to distinguish between those two dating hypotheses. We then discussed the previously raised discrepancies between the different paleoenvironmental proxies in the light of the results of climate simulations.

El Harhoura 2 cave
The archeological site EH2 is located on the Moroccan Atlantic coast in the Rabat-Témara region (33°57'08.9" N / 6°55'32.5" W). Dating, stratigraphic and paleoenvironmental information are summarized in Fig1. This site is dated from the Late Pleistocene to the mid-Holocene. Its stratigraphy is currently divided into 11 layers (each layer is abbreviated in "L" followed by the layer number), 45 among which the first eight are well studied and considered in this paper. Those layers revealed an impressive taxonomic richness and delivered an important amount of large and small vertebrate remains (Michel et al., 2009;Stoetzel et al., 2011Stoetzel et al., , 2012b. The site was used several time as a model to explore the relationship between the past diversity, the phenotype of organism and their environments (Stoetzel et al., 2010(Stoetzel et al., , 2011(Stoetzel et al., , 2012a(Stoetzel et al., , 2013(Stoetzel et al., , 2017Campmas et al., 2015;Cornette et al., 2015;Terray et al., 2021).

50
At EH2, paleoenvironments have mainly been inferred based on two different kind of proxies: species presence and isotopes. Regarding species presence, environments were reconstructed using palaeoecological indices such as the Taxonomic habitat index (THI), which is a method based on the habitat preference of small vertebrates (Stoetzel, 2009;Stoetzel et al., 2011Stoetzel et al., , 2014. The presence and/or abundance of particular taxa can also be a strong indicator of certain types of environment, such as 55 https://doi.org/10.5194/cp-2021-185 Preprint. Discussion started: 5 January 2022 c Author(s) 2022. CC BY 4.0 License. amphibians for humid context, or gerbils and gerboas for arid context. Isotope-based inferences were performed by Jeffrey (2016) on Meriones teeth from Layers 2 to 8 of EH2, providing varied indications about paleoenvironments such as aridity, relative humidity, seasonal variations and vegetal cover (Longinelli and Selmo, 2003;Navarro et al., 2004;Royer et al., 2013). Some inconsistencies were observed between species and isotopes based proxies on certain stratigraphic layers of the site (Jeffrey, 60 2016;Stoetzel et al., 2019). On two layers (5 and 7), while THI indicates arid conditions, isotopes suggest a more humid and temperate climate. Given the mosaic-like landscape described at EH2 and the fact that the two methods do not record the same information (more global for faunal communities, at the scale of a limited number of individuals of a single species for isotopes), it is not surprising to meet some inconsistencies, but it makes it difficult to reliably infer the global environmental conditions of 65 the site.
The current climatic context of EH2 is complex. Morocco, due to its geographical location, has a rich mosaic of climate types (Sobrino and Raissouni, 2000). In particular, the Rabat-Témara region where the cave is located is subject to climatic influences from the Atlantic Ocean in the west, the arid Sahara in the south and the Mediterranean region in the north. Its climate is marked by a strong annual and 70 interannual irregularity with a hot and dry summer where precipitations are almost absent and evaporation is particularly high, and a soft and short winter with concentrated precipitations (Sobrino and Raissouni, 2000;Lionello et al., 2006). The site also suffered from important climatic fluctuations over time (e.g. Hooghiemstra et al., 1992;deMenocal, 1995deMenocal, , 2004Le Houérou, 1997;Carto et al., 2009;Trauth, Larrasoaña, & Mudelsee, 2009;Drake et al., 2011;Drake, Breeze, & Parker, 2013;Blome et 75 al., 2012;Kageyama et al., 2013;Couvreur et al., 2020). The Late Pleistocene -Holocene period (Marine Isotopic Stages (MIS) 5 to 1) is marked by significant global climate changes, which resulted in a succession of relatively humid/arid and open/closed environments at EH2 (Stoetzel, 2009;Stoetzel et al., 2011Stoetzel et al., , 2012a. As a result, paleolandscapes of the Late Pleistocene are described as open steppe or savanna-like lands with patches of shrubs, woodlands and water bodies, the latter expanding during    Pierre-Simon Laplace) model (Marti et al., 2010;Dufresne et al., 2013;Boucher et al., 2020). All 100 versions consider the full coupling between ocean, atmosphere, land-surface and sea-ice components, as well as an interactive carbon cycle. The exact complexity and resolution vary depending on model version. The different simulations also differ by Earth's orbit parameters, atmospheric trace gazes composition and ice sheet configurations in order to represent the climate conditions of the different periods (Table1). We used a total of six simulations:

105
-For the climate in L1 we used the reference midHolocene IPSL simulation (mean age estimated at ~6ka) run as part of the PMIP4-CMIP6 ensemble of simulations (Kageyama et al., 2017;Braconnot et al., 2021) .It is called CM6mh in the following and has been performed with IPSL-CM6A-LR (Boucher et al., 2020). The atmosphere has a resolution of 144 points in longitude, 143 points in latitude and 79 vertical levels (144x143xL79). The ocean has a nominal 1°resolution with refined grid in the equatorial 110 region and 75 vertical levels.
Unfortunately, these simulations are not directly comparable because they were performed using different versions of the global coupled climate model. Indeed, those different versions, when used to simulate current climate, showed slightly different biases when compared to observations. This is due to the fact that the various versions of the model are characterized by different physical representations, 125 resolutions and tuning. The differences between the collected simulations could then result from the difference in bias between the various versions of the model instead of representing climate differences between periods.
In order to make reliable comparisons between periods, we ran an entirely new set of simulations for the key periods of EH2: midH (for the mid-Holocene period), earlyH (for the early Holocene period), 130 midMIS3 (for the mid MIS3 period), lateMIS4 (for the late MIS4 period), midMIS4 (for the mid MIS4 period) and MIS5d (for the MIS5d period). Configurations information are summarized in Table1.

Model
We ran those new simulations using a unique model: LMDZOR6A. This model is a sub-configuration 135 of the coupled model IPSL-CM6v1.11-LR, it considers the atmosphere-land surface component of the IPSL-CM6A-LR coupled model. It was chosen because it best simulates the climate on our area of interest (SM Fig 1 and 2). It also has the finest grid and improved representation of the atmospheric and land surface processes.
First, we produced a control simulation Ctrl on current days using the clim_pdControl experiment. This  However, no SST reconstructions are available for the desired periods. Consequently, we imposed the simulated SST from pre-existing coupled simulations. However, those SST are also subjected to biases 150 due to different model versions and they must be corrected. First, we estimated the differences between the SST biases of the three model versions. To do that we compared the SST simulated per each coupled model on current days with the SST of clim_pdControl. Those SST are the ones used for AMIP simulations (Atmospheric Model Intercomparison Project), which are a mean annual cycle of SST from observations. Secondly, we explored the consequences of those biases on the simulated climatology.

155
We ran three test simulations StanCM5ASST, StanCM5A2SST and StanCM6ASST. They are similar to Ctrl except that we replaced the default SST (the SST of AMIP simulations) with the SST simulated per each of the three coupled model versions.
Finally, we computed the anomaly of each model version, i.e. the difference between the simulated SST and the SST used for AMIP simulations (issued from observations). For each coupled simulation, we 160 corrected the simulated SST by the anomaly corresponding to the model version used. Corrected SST were obtained according to the formula: Where SSTamip are the SST used for AMIP simulations, SSTsim the SST from the coupled simulation, SSTmod the SST of the model version for current days and SSTcor the corrected SST.

165
In total, the simulations ran for 50 years. Three monitoring variables were used: bils (surface total heat flux), nettop (net dn radiative flux at top of the atmosphere) and mrso (total soil moisture). Bils and nettop allow to ensure that the energy balance of the system has reach equilibrium, and mrso allow to verify the stabilization of the hydrological reservoir. In global the model reach equilibrium after eight years of simulation and at local scale on EH2 after 17 years (SM Fig 3). Consequently, we worked the    w10m (wind at 10 meters) (Pellegrino et al., 2013;Chapman et al., 2011;Tanner et al., 1991;McNeil, 1991), sols (solar radiation at surface) (Monteith, 1972;Fyllas et al., 2017), drysoil_frac (fraction of 185 visibly dry soil) (Paz et al., 2015) and humtot (total soil moisture) (Paz et al., 2015). Two additional variables were computed. The daily temperature amplitude tsol_ampl_day (Alhajeri and Steppan, 2016) computed from tsol_max (day maximum temperature) and tsol_min (day minimum temperature) as: tsol_max -tsol_min; and the hydric stress hyd_stress (Martínez-Blancas and Martorell, 2020) computed from evapot (potential evaporation) and evap (evaporation) as: evapotevap. List of abbreviations and 190 units is provided in Table 2. For each variable, we worked on a seasonal cycle averaged on the last 30 years of simulation. We explored variations in monthly means (mean value of the variable) and monthly standard deviations (amplitude of seasonal variation). Climate variations over the EH2 sequence are presented on a common scale to ease the reading and better visualize covariation between climate variables (all variables were centered and reduced). In additions, to visualize climate proximity/differences between EH2 layers we performed principal component analyses (Jolliffe and Cadima, 2016). It allowed us to visualize EH2 layers into a climate-200 space by presenting the data along the principal components. These are new uncorrelated variables that successively maximize variance. They are computed from the eigenvectors and eigenvalues of the covariance/correlation matrix. The proximity of layers in the climate-space represents the similarity between the climate of the layers.
We tested the congruence between paleoenvironmental information available in the literature and our 205 results from climatic simulations. As climate variables we used the principal components obtained before, and as paleoenvironmental proxies we used isotopes (δ 18 O and δ 13 C) means (Jeffrey, 2016), reconstructed mean annual precipitations (Jeffrey, 2016) and percentages of represented habitats indicated by the THI (Stoetzel et al., 2014) presented in Fig1. First, to test the global covariance between climate simulations and paleoenvironmental proxies we carried out two-block partial least squares (two-210 block PLS), a method for exploring the patterns of covariation between two sets of variables (Sampson et al., 1989;Streissguth et al., 1993). Then, to refine our results, we performed pairwise correlation tests.
Regarding the THI, the oasis habitat was not tested because its percentage of presence does not vary over the EH2 sequence.  temperature (ts) and precipitations (pr) (Boucher et al., 2020). In each case a cold bias of ~5°C is observed over our area of interest. Some dry bias in precipitations are also noticed in winter. CM5 presents the strongest cold bias, up to 6 °C. CM5A2 display highly variable seasonal biases, with a cold one that reaches 6°C in winter and a hot one that reaches 2°C in summer. The climate simulated with 240 the SST from CM6 is the closest to the climate simulated with the SST of AMIP.  The last climatic transition, subtler, is between L1 and Act. The environment in Act seems close to the 295 one in L8 with more seasonal variability in temperature, solar radiations and water stress.
Overall, it seems there is an alternation of two main climate types. This partition is confirmed by the principal component analysis showed in Fig8. The first principal component explains 65.85% of the observed variance and split the layers of EH2 into two climate groups. The first regroups L3, L4a, L6 and L7 and is defined by humid and windy conditions with a high seasonal variability. The second 300 includes Act, L1, L2 and L8 and is characterized by hot and dry conditions. The second axis, explaining 30.28% of the variability, divides the latter group into two subgroups: L1 and L2 with high seasonal variability, and Act and L8 with a lower seasonal variability.
Regarding DH2, we observe three important climatic transitions (Fig7). The first happened between L5 and L4a. In L4a, the climate is much windier and temperature are colder with a substantial increase of 305 daily temperature variations. Precipitations and soil moisture increase importantly while hydric stress decreases. The climate also presents an overall higher seasonal variability. These tendencies are even more pronounced in L3. The second transition is between L3 and L2. The soil is drier and the hydric stress increases greatly as well as solar radiations and temperatures. Precipitations and soil moisture are less important. Conditions in L1 are close from those in L2, a bit colder with a relatively less important 310 water stress. The last climate transition, slighter than the previous ones, occurs between L1 and Act and is mainly marked by a global decrease of seasonal variation.
Results associated to DH2 suggest that three types of climates succeeded one another at EH2. The first group is composed of L8, L7, L6 and L5, the second by L3 and L4a and the third by L2, L1 and Act. As  (2016)); Forest, Bush, Steppe, Wetland, Rocky: relative % of representation according to the THI (from Stoetzel et al. (2014) and Jeffrey (2016)).

4 Discussion
Paleoclimate simulations allow us to discuss several important climate changes at EH2 over the Late Pleistocene to mid-Holocene period. Isotopes are correlated to particular climate variables of DH1. Both δ 18 O and δ 13 C are mainly related to seasonal variation in temperature and water stress. It is known that a high amount of δ 18 O is an indicator of aridity (Longinelli and Selmo, 2003;Levin et al., 2006 The large differences between the climate simulations and the fact that they provide a consistent view of the relationship between the different climatic variables allow us to discuss the inconsistencies existing between paleoenvironmental proxies at EH2. The two major ones concern L5 and L7. In both 385 cases, isotope surveys and mean annual precipitations reconstructions from Jeffrey (2016) indicate humid conditions with important precipitations. On the contrary, the THI as well as the presence of the steppic species Jaculus cf. orientalis (often used as an indicator of particularly arid conditions) and the scarcity of aquatic species support a dry climate (Stoetzel, 2009;Stoetzel et al., 2014). Large mammals would also support this hypothesis, with an increase in the representation of gazelles and alcelaphines, 390 and a decrease in the representation of bovines in both layers (Stoetzel et al., 2012a(Stoetzel et al., , 2014

395
Concerning L7, the climate describe by DH1 agrees with isotope surveys (Jeffrey, 2016). Those conclusions are supported by the increased abundance of Crocidura russula, a shrew species associated with Mediterranean climates (Cornette et al., 2015). In addition, Jaculus cf. orientalis can also be considered as an indicator of more continental conditions, such as the distance from the coastline, rather than a marker of arid environments.

400
An important difference is noticed between the climate describe by DH1 and the THI. DH1 describes an arid climate on L1, in contradiction with the composition of small and large mammal communities

410
The second issue relates to the compatibility of the spatial and temporal resolution of a climate simulation and a stratigraphic layer of an archeological site. The spatial resolution of the atmospheric grid of the IPSL-CM6A-LR model is 157 km in average (Boucher et al., 2020). Consequently, the four cell grids used in this study cover a quite wide area. Conversely, EH2 represents a precise locality, and most species whose presence was recorded have a lifetime dispersal range largely inferior to 157km 415 (e.g., the jird Meriones shawii has a home range estimated between 200-1000 m 2 (Ghawar et al., 2015)).
Then, the climate described by simulations might be a bit too coarse to faithfully describe microclimate variations at EH2. This may induce a mismatch between climate variables and other paleoenvironmental proxies.
The difference in temporal resolution between a stratigraphic layer and a paleoclimate simulation raises 420 even more concerns. The older the dated layer, the higher the dating uncertainty is. For EH2, this uncertainty goes until 11 ka on L8. If we compare the dates estimated by different methods, the difference is much more significant. For example, it is the case for L3 dated at ~40 ka by combined US-ESR and at ~60 by OSL (i.e. a difference of 20 ka). Radical climate changes can occur in 20 ka: ~60 ka places L3 in the last ice age while ~40 ka places it in the interglacial. As this paper demonstrates by 425 comparing DH1 and DH2, those differences in estimated dates can result in completely different climate sequence. EH2 was an ideal case because a lot of studies focusing on the dating of its stratigraphic layers are found in the literature (e.g., Jacobs et al., 2012;Janati-Idrissi et al., 2012;Ben Arous et al., 2020a,b).
However, for archeological and paleontological sites, whose chronological framework is less referenced, it can lead to a completely wrong climate description. Furthermore, layers are not necessarily snapshots representing the faunal assemblages, cultural presence, etc. of specific moments in the past.
A layer is a stratigraphic/sedimentary unit that can cover more or less long periods and undergo microclimatic variations that cannot be disentangled.

435
The use of climate simulations allowed us to provide a quantify description of the paleoclimate at El Harhoura 2. Our results enabled to discuss the varied dating methods used: in our case, the climate sequence corresponding to combined US-ESR ages was much more consistent with paleoenvironmental inferences than that corresponding to OSL ages. Regarding discrepancies between paleoenvironmental inferences based on faunal assemblages and isotopic studies, our results are more consistent with 440 conclusions drawn from the latter. But more than that, they highlighted the difference in scale between the information provided by each of these indicators. This study demonstrates that the combination of different sources of environemental data and climate simulations has a great potential for refining paleoenvironmental and chronological context of archeological and paleontological sites.

445
Code and data availability The code and data used in this article are available on reasonnable request to the corresponding author.