# The use of paleoclimate simulations to refine the environmental and chronological context of archaeological/paleontological sites

Léa Terray, Emmanuelle Stoetzel, Eslem Ben Arous, Masa Kageyama,
Raphaël Cornette, Pascale Braconnot

<figure><img src="image_1.png" /></figure>

## To cite this version:

Léa Terray, Emmanuelle Stoetzel, Eslem Ben Arous, Masa Kageyama, Raphaël Cornette, et al.. The use of paleoclimate simulations to refine the environmental and chronological context of archaeological/paleontological sites. 2023. hal-03992234

HAL Id: hal-03992234

https://hal.science/hal-03992234

Preprint submitted on 6 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

Léa Terray<fnref n="1" />, Emmanuelle Stoetzel<fnref n="2" />, Eslem Ben Arous<fnref n="3" />,<fnref n="4" />,<fnref n="2" />, Masa Kageyama<fnref n="5" />, Raphaël Cornette<fnref n="1" /> and
Pascale Braconnot<fnref n="5" />

<footnote n="1">¹ Institut de Systématique, Evolution, Biodiversité (ISYEB) – UMR 7205, Muséum National d'Histoire Naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, Paris, France.</footnote>

<footnote n="2">² Histoire Naturelle de l'Homme Préhistorique (HNHP) – UMR 7194, CNRS, Muséum National d'Histoire Naturelle, Sorbonne Université, UPVD, Musée de l'Homme, Paris, France.</footnote>

<footnote n="3">³ Pan-African Evolution Research Group (Pan-Ev), Max Planck Institute for the Science of Human History, Jena, Germany.</footnote>

<footnote n="4">⁴Geochronology lab, Centro Nacional de Investigación sobre la Evolución Humana, Burgos, Spain</footnote>

<footnote n="5">⁵Laboratoire des Sciences du Climat et de l'Environnement (LSCE) – UMR 8212 / Institut Pierre Simon Laplace (IPSL) – UMR 8112, CEA/CNRS/UVSQ, Centre CEA-Saclay, Gif-sur-Yvette, France</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

**Abstract.**

This study illustrates the strong potential of combining paleoenvironmental reconstructions and paleoclimate modeling to refine the paleoenvironmental and chronological context of archaeological and paleontological sites. We focus on the El Harhoura 2 cave (EH2), an archeological site located on the North-Atlantic coast of Morocco that covers a period from the Late Pleistocene to the mid-Holocene. On several stratigraphic layers, inconsistencies are observed between species- and isotope-based inferences used to reconstruct paleoenvironmental conditions. The stratigraphy of EH2 also shows chronological inconsistencies on older layers between age estimated by Optical Stimulated Luminescence (OSL) and Combination of Uranium Series and Electron Spin Resonance methods (combined US-ESR). We performed paleoclimate simulations to infer the global paleoclimate variations over the EH2 sequence in the area, and we conducted a consistency approach between paleoclimate reconstruction estimated from simulations and available from EH2 paleoenvironmental inferences. Our main conclusion show that the climate sequence based on combined US-ESR ages is more consistent with paleoenvironmental inferences than the climate sequence based on OSL ages. We also evidence that isotope-based inferences are more congruent with the paleoclimate sequence than species-based inferences. These results highlight the difference in scale between the information provided by each of these paleoenvironmental proxies. Our approach is transferable to other sites due to the increase number of available paleoclimate simulations.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

# 1 Introduction

The reconstruction of paleoenvironments has long been a subject of great interest, particularly to study the past biodiversity. Archeological sites provide unique opportunities to infer paleoenvironments from faunal and/or vegetal remains (e.g., Avery, 2007; Denys et al., 2018; Stoetzel et al., 2011; Comay and Dayan, 2018; Matthews, 2000; Marquer et al., 2022) or stable light isotope composition of sediments or biological remains (e.g. Tieszen, 1991; Royer et al., 2013). These approaches allow one to characterize the past landscapes and biotic environments on the localities of the sites. However, inconsistencies between isotope compositions, the different types of remains and available stratigraphy are often to be faced and prevent to properly assess the possible relationships between biodiversity and paleoenvironmental changes. This is the case for the El Harhoura 2 (EH2) cave, an archeological site located on the North-Atlantic coast of Morocco, that we use here as a case study.

At EH2, paleoenvironments have mainly been inferred based on two different kinds of proxies: species presence and stable light isotope composition. Because of species habitat preferences, the presence and/or abundance of particular taxa is a strong indicator of certain types of environment, such as amphibians for more humid contexts, or gerbils and jerboas for more arid contexts (e.g., Fernandez-Jalvo et al., 1998; Stoetzel et al., 2011). Stable light isotope composition of small mammals teeth provide varied indications about diet and paleoenvironments (Longinelli and Selmo, 2003; Navarro et al., 2004; Royer et al., 2013), and thus about aridity, seasonal variations and vegetal cover. At EH2, inconsistencies are observed between species- and isotope-based inferences (Jeffrey, 2016; Stoetzel et al., 2019). On two stratigraphic layers, while species presence suggests drier conditions than usual, isotopic composition of Meriones teeth suggest a more humid and temperate climate. Such discrepancies are not fundamentally surprising, because species presence and isotopic composition do not deliver the exact same information: species presence carries a signal at the scale of faunal communities, while isotopic composition reflects the diet preferences of a limited number of individuals of a single species. But these mixed messages make it difficult to extrapolate the global environmental conditions of the site.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

The stratigraphy of EH2 also shows some chronological inconsistencies. Radiocarbon dating (AMS-
¹⁴C) is the most reliable dating method, however it cannot date remains older than 50 ka. Beyond this
reach, other dating methods must be applied, such as Optical Stimulated Luminescence (OSL) and
Combination of Uranium Series and Electron Spin Resonance methods (combined US-ESR). It has been
shown that OSL and combined US-ESR methods display important differences in the estimated age of
the same stratigraphic layers (Ben Arous et al., 2020b). These differences are related to the fact that
these dating methods do not date the same objects. OSL estimates the last time quartz sediment was
exposed to light, while combined US-ESR estimates the age of fossil teeth. To date, the respective
reliability of these methods is difficult to establish.

These kind of paleoenvironmental and chronological discrepancies are widespread in archeology and
paleontology. Model simulations may inform us about the broader climate influences over the region,
and therefore might enrich our understanding of the large-scale climate changes. Climate models allow
to simulate paleoclimates using the physical laws that describe the dynamics and thermodynamics of
the Earth system, and model data comparisons have shown that the large scale patterns are consistent
with paleoclimate reconstructions, even though regional feature are underestimated or affected by model
biases (Kutzbach and Otto-Bliesner, 1982; Braconnot et al., 2012; Duplessy and Ramstein, 2013;
Schmidt et al., 2014; Harrison et al., 2015). They thus offer a consistent framework to test the
consistency between climate drivers and environmental changes recorded at EH2 depending for the
different stratigraphy options.

In this paper, we make use of simulations run over the last 10 years with the IPSL model to derive a set
of atmosphere alone simulations with the latest version of the model (Boucher et al., 2020; Hourdin et
al., 2020). This version is chosen because it of its skill in reproducing the climate in the region. However,
because of the differences in scale and resolution between the archaeological record and the
paleoclimate simulations, several caveats must be considered beforehand. Firstly, archaeological and
paleoclimate data present different temporal resolution. Paleoclimate simulations present snapshots of
the hypothetical climate states at precise ages. Conversely, an archeological layer is a
stratigraphic/sedimentary unit that can cover shorter or longer periods and undergoes microclimatic
<figure><img src="image_10.png" /></figure>

variations that cannot be disentangled. Nevertheless, because EH2 is dated from the Late Pleistocene to the mid-Holocene period (Marine Isotopic Stages (MIS) 5 to 1), which in the area is marked by important global climatic fluctuations over time (i.e. the last glacial-interglacial transition; e.g., Hooghiemstra et al., 1992; deMenocal, 1995, 2004; Le Houérou, 1997; Carto et al., 2009; Trauth et al., 2009; Drake et al., 2011, 2013; Blome et al., 2012; Kageyama et al., 2013; Couvreur et al., 2020), we do not expect these differences to prevent a global consistency approach. Secondly, archaeological and paleoclimate data present different spatial resolution. The spatial resolution of the atmospheric grid use here is ~150 km (Boucher et al., 2020). Conversely, EH2 represents a precise locality, and most species whose presence was recorded have a lifetime dispersal range largely inferior to 150 km (e.g., the jird Meriones shawii has a home range estimated between 200-1000 m² (Ghawar et al., 2015)). 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. We make the hypothesis that the results of the global climate simulations are sufficient to capture the large climate changes we are interested in. Regional dynamical or statistical approaches would be derived from these simulations and would add additional unknowns due to the fact that the tuning of these methods has to be done using present day observations at the site. The cave is now imbedded in an urban area that experiences rapid climate and environmental changes over the last century and it is unclear that tuning over this recent period would be valid for the past conditions we are considering in this study.

In order to discuss and refine the paleoenvironmental and chronological context of EH2, we conduct a consistency approach between paleoclimate reconstruction from simulations and paleoenvironmental inferences based on EH2 content from the literature. To overcome the issue of the differences between dates estimated with different dating methods, we considered two separate chronological sequences based on different dating method. Finally, for each of the chronological sequences, we test the consistency between paleoclimate variables extracted from simulations, and species presence and stable light isotope composition. We expect this consistency approach to discriminate paleoenvironmental inconsistencies between species- and isotope-based proxies, and also eventually to distinguish between the two chronological sequences, and therefore between the two dating methods.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

The remainder of the manuscript is organized as follow. We present the EH2 cave and the different choices made to run the climate set of paleoclimates in section 2. In section 3 we present the consistency analyses and the results. The discussion in section 4 highlights the major results and the proof of concept of the proposed approach, before the general conclusion (section 5).

## 2 Material and methods

### 2.1 El Harhoura 2 cave

#### 2.1.1 Presentation of the site

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). The climate of the area is marked by a strong seasonal variability (Fig 1) with a relatively warm and dry summer during which precipitations are almost absent and evaporation is particularly high, and a relatively cool and short winter, which is also the rainy season (Sobrino and Raissouni, 2000; Lionello et al., 2006). Inferences based on species presence and abundances suggest that, in the past, the site underwent important climatic fluctuations, which resulted in a succession of relatively humid/arid and open/closed environments at EH2 (Stoetzel, 2009; Stoetzel et al., 2011, 2012a, b). 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 wet periods, especially during the mid-Holocene (Stoetzel, 2009; Stoetzel et al., 2012a, 2014).
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

100

<figure><img src="image_15.png" /><figcaption>Fig. 1. Maps of mean temperature (unit: °C), wind speed and direction at 850hPa (represented by arrows, length is proportional to wind speed) and corrected pressure at sea level (unit: hPa) per season of the global atmospheric reanalysis ERA-interim (ERAi; Berrisford et al., 2011) and the *historical* simulations of IPSL-CM6-LR models on the region of EH2. Data were averaged on 30 seasonal cycles (1980-2009). EH2 cave location is represented by a star in the upper left panel. DJF: December, January, February (winter); JJA: June, July, August (summer).</figcaption></figure>

105
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 2.** Summary diagram displaying stratigraphy, paleoenvironmental proxies and the two dating hypotheses of El Harboura 2 (EH2). **A:** stratigraphy of EH2, unused layers are in dark grey; **B:** relative % of THI values (adapted from Stoetzel *et al.* 2014) and Jeffrey (2016); **C:** mean $\delta^{18}\text{O}$ values in *Meriones* teeth (from Jeffrey 2016); **D:** mean $\delta^{13}\text{C}$ values in *Meriones* teeth (from Jeffrey 2016); **E:** Mean Annual Precipitations (MAP; from Jeffrey (2016)); **F:** dating hypotheses DH1 and DH2 for the different layers of EH2 (Nespoulet & El Hajaoui, 2012; Jacobs *et al.*, 2012; Janati-Idrissi *et al.*, 2012; Ben Arous *et al.*, 2020b, a; Marquer *et al.*, 2022).</figcaption></figure>
<figure><img src="image_19.png" /></figure>

### 2.1.2 Chronostratigraphy and dating hypotheses

EH2 displays a high resolution stratigraphy and its layers have revealed an impressive taxonomic richness and delivered an important amount of large and small vertebrate remains (Michel et al., 2009; Stoetzel et al., 2011, 2012b). Its stratigraphy is currently divided into 11 layers (Fig2 A) (each layer is abbreviated as “L” followed by the layer number), among which eight are well studied and considered in this study. All of these eight layers have been dated (Ben Arous et al., 2020b). Three different methods were used: AMS-¹⁴C (Nespoulet and El Hajraoui, 2012; Marquer et al., 2022), OSL (Jacobs et al., 2012) and combined US-ESR (Janati-Idrissi et al., 2012; Ben Arous et al., 2020a). AMS-¹⁴C is the most reliable dating method of the three and is used for recent layers (L1 and L2). However, for older layers beyond the reach of AMS-¹⁴C (L3, L4a, L5, L6, L7 and L8), OSL and combined US-ESR method are used. However, when applied to a same layer, these two methods present discrepancies. For example, it is the case of L5 dated at ~60 ka by combined US-ESR and at ~100 ka by OSL (Fig2 F). In addition, when two consecutive layers are dated with these different methods, dates can be inconsistent with the relative position of the layers in the stratigraphy, as it is the case for the L8 dated at ~90 ka using combined US-ESR and L7 at ~110 ka using OSL (Fig2 F). To overcome this issue, we choose to set two separate dating hypothesis (DH) for the stratigraphic sequence: DH1 refers to the chronological sequence based on AMS-¹⁴C and combined US-ESR dates, and DH2 on AMS-¹⁴C and OSL dates. The DH are presented on Fig2 F.

110

115

120

125

### 2.1.3 Paleoenvironmental variables

Paleoenvironments at EH2 have been inferred based on two different kinds of proxies: species presence and stable light isotope composition. Regarding species presence, paleoenvironments are reconstructed using the Taxonomic habitat index (THI). The THI is a palaeoecological index which allows one to reconstruct the paleolandscape based on the species composition of the community (Stoetzel, 2009; Stoetzel et al., 2011, 2014). Each species is associated to its preferred habitat through an actualist approach. From that is deduced a qualitative composition of the paleolandscape, expressed as a

130
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

percentage. Data are from Stoetzel et al. (2014) and are presented in **Fig2 B**. Note that the oasis habitat
was not considered because its percentage of presence does not vary over the EH2 sequence.

135

Isotope-based inferences are from *Meriones* teeth. They provide information about the diet and the environment of studied individuals (Longinelli and Selmo, 2003; Navarro et al., 2004; Royer et al., 2013). Two isotopic fractions are considered: δ¹⁸O and δ¹³C. Plants consumed by small mammals are sensitive to environmental conditions which shows up in their δ¹³C and δ¹⁸O values. Because of that, they are good indirect indicators of aridity, seasonal variation and vegetal cover (Longinelli and Selmo, 2003; Blumenthal et al., 2017; Blumenthal, 2019). We used δ¹⁸O and δ¹³C means, minimums and maximums, and reconstructed mean annual precipitations (MAP) computed from δ¹⁸O from Jeffrey (2016). Data are presented in **Fig2 C, D** and **E**.

140

145

**2.2 Paleoclimate reconstruction**

**2.2.1 Climate model**

We ran a new set of paleoclimate simulations covering the Late Pleistocene to mid-Holocene period using the model LMDZOR6A (Hourdin et al., 2020). This model is the atmosphere-land surface component of the IPSL-CM6A-LR coupled model (Boucher et al., 2020) that has been used to run the CMIP6 ensemble of past, present and future climate simulations (Eyring et al., 2016), including the mid-Holocene (Braconnot et al., 2021), last interglacial (Sicard et al., 2022; Otto-Bliesner et al., 2021) and Pliocene (Haywood et al., 2016) periods. It has an atmospheric resolution of 144 points in longitude, 143 points in latitude and 79 vertical levels (144x143xL79). Compared to previous IPSL model version it has a finer spatial and vertical resolution and an improved representation of atmospheric and land surface processes. In addition, the simulated climatology of temperature and precipitation over the region is consistent with observations when the model is run for present day conditions with sea-surface temperature prescribed to the monthly climatology of the observed Atmospheric Model Intercomparison Project (AMIP) SST fields (Boucher et al., 2020, 2018) (Fig 1 and 3). Comparison of the simulated annual mean cycle of surface air temperature and precipitation match quite well the observed one,

150

155
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

160

despite an overestimation of summer temperature and slight 0-1 month shift in the seasonal cycle of precipitation (Fig3).

<figure><img src="image_24.png" /><figcaption>**Fig. 3.** Monthly variations in mean temperature (unit: C°) and precipitations (unit: mm/month) of the ERAi and GPCP reanalyzed data, the AMIP CM6A simulations (atmosphere only) and *historical* coupled simulations of IPSL-CM5A-LR, IPSL-CM5A2-LR and IPSL-CM6-LR models on the four grid cells containing EH2 cave. Data were averaged on 30 seasonal cycles (1980-2009). Error bars are estimated from interannual variation over the averaged 30 years and is visualized by quartiles. As a reference for current climate we used monthly data from the global atmospheric reanalysis ERA-interim for monthly mean temperature (Berrisford et al., 2011) and from the GPCP v2.3 (Global Precipitation Climatology Project) for monthly mean precipitations (Adler et al., 2018).</figcaption></figure>

165

170

### 2.2.2 Paleoclimate simulations

We produced a total of six paleoclimate simulations with the LMDZOR6A model for the key periods of the EH2 sequence: *midH* (for the mid-Holocene period), *earlyH* (for the early Holocene period), *midMIS3* (for the mid MIS3 period), *lateMIS4* (for the late MIS4 period), *midMIS4* (for the mid MIS4 period) and *MIS5d* (for the MIS5d period). The different simulations differ by prescribed Earth's orbital parameters, atmospheric trace gases composition, ice-sheet configuration and sea surface temperature (SST) in order to represent the climate conditions of the different periods (see **Table 2** for details). For

175
<figure><img src="image_25.png" /></figure>

this we make use of preexisting paleoclimate and control simulations that have been run in the last 8 years with different versions of the IPSL model (Marti et al., 2010; Dufresne et al., 2013; Boucher et al., 2020) (details are available in **Table 1**).

<table><caption>**Table 1.** Preexisting coupled simulations from the IPSL repository.</caption><thead><tr><th>Model</th><th>IPSL-CM6A-LR</th><th>IPSL-CM5A-LR</th><th>IPSL-CM5A-LR</th><th>IPSL-CM5A-LR</th><th>IPSL-CM5A-LR</th><th>IPSL-CM5A2-LR</th></tr></thead><tbody><tr><th>Date</th><td>~ 6 kyr BP</td><td>~ 9 kyr BP</td><td>~ 40 kyr BP</td><td>~ 60 kyr BP</td><td>~ 66 kyr BP</td><td>~ 115 kyr BP</td></tr><tr><th>References</th><td>Kageyama et al. (2017)<br>Braconnot et al. (2021)</td><td>Le Mézo et al. (2017)</td><td>Le Mézo et al. (2017)</td><td>Le Mézo et al. (2017)</td><td>Le Mézo et al. (2017)</td><td>Sicard personal communication</td></tr></tbody></table>

We also ran a control simulation *Ctrl* (see **Table 2** for details), representative of the present-day climate. The Sea Surface Temperature (SST) boundary conditions used in *Ctrl* correspond to the mean annual cycle of the SST estimated from current observations used for the AMIP (Boucher et al., 2018) and repeated in time. This *Ctrl* simulation will be considered as the reference for the current climate in our ensemble of new atmospheric simulations.

Unfortunately, simulated paleoclimate SST are not directly comparable because they were performed using different versions of the IPSL model (**Table 1**). This is due to the fact that the various versions of the model are characterized by different physical representations, resolutions and tuning. In particular, these different model versions show different present-day SST biases when compared to observations (**Fig 4**). They translate in different representation of the seasonal cycle of surface air temperature and precipitation at EH2 in simulations with LMDZOR6 when used as boundary condition instead of the AMIP SST field (**Fig 3**). The magnitude of the seasonal cycle is overestimated and precipitation are underestimated with almost no seasonality when IPSLCM5A or IPSLCM5A2 control SST are used as boundary condition. This differences that can be attributed exclusively to the differences in the SST biases are significant and may deteriorate the simulation of other variables of interest, thus complicating the intercomparison between our new LMDZOR6A simulations. Indeed, the differences between these simulations for the periods of interest at EH2 could then result from the difference in bias between the various versions of the model rather than representing significant climate differences between periods.
<figure><img src="image_26.png" /></figure>

205 In order to produce an homogeneous set of simulations we thus apply a correction to the SST fields to that the results better fit the observations for present day and the information of the SST changes for the different period is kept as part of the prescribed paleoclimate SST boundary conditions.

<figure><img src="image_27.png" /></figure>

210 **Fig. 4.** SST biases of IPSL-CM5A-LR, IPSL-CM5A2-LR and IPSL-CM6A-LR relative to AMIP's SST (issued from current observations) (unit: °C). The seasons are DJF: December, January, February (winter); MAM: March, April, May (spring); JJA: June, July, August (summer); SON: September, October, November (autumn).

215 **2.2.3 Sea-surface boundary conditions**

We correct the simulated SST of each coupled simulation for the systematic bias of the model. This is done by removing the SST bias corresponding to the model version used. In other words, corrected SSTs were obtained according to the formula:

$$SST_{cor} = SST_{sim} + (SST_{amip} - SST_{mod})$$

220 Where SSTsim are the SSTs from the coupled simulation, SSTamip are the AMIP SSTs, SSTmod the SSTs of the model version for current days and SSTcor the corrected SSTs, which will be used as boundary conditions in our new LMDZOR6A simulations. The underlying hypothesis in this correction scheme is that the mean annual SST bias cycle, e.g. SSTmod – SSTamip for each, is stationary in time.
<figure><img src="image_28.png" /></figure>

Note also that boundary conditions files of pre-existing (coupled) simulations were interpolated on a 143x144xL79 grid to be compatible with the grid of LMDZOR6A.

Configuration details for the new set of simulations are summarized in **Table 2**. The length of all these new simulations is 50 years, which is long enough given the fact that LMDZOR6A is an atmospheric model. In the following annual mean cycles are estimated from the last 30 years of each LMDZOR6A simulation.

<table><caption>**Table. 2.** Forcing and boundary conditions of the simulations produced in this study.</caption><thead><tr><th colspan="8">Simulations</th></tr><tr><th>Name</th><th>Ctrl</th><th>midH</th><th>earlyH</th><th>midMIS3</th><th>lateMIS4</th><th>midMIS4</th><th>MIS5d</th></tr></thead><tbody><tr><td>Model</td><td>LMDZOR6A</td><td>LMDZOR6A</td><td>LMDZOR6A</td><td>LMDZOR6A</td><td>LMDZOR6A</td><td>LMDZOR6A</td><td>LMDZOR6A</td></tr><tr><td>Date</td><td>Current days</td><td>~ 6 kyr BP</td><td>~ 9 kyr BP</td><td>~ 40 kyr BP</td><td>~ 60 kyr BP</td><td>~ 66 kyr BP</td><td>~ 115 kyr BP</td></tr><tr><td colspan="8">Orbital parameters*</td></tr><tr><td>Eccentricity</td><td>Same as<br>clim_pdControl</td><td>0.018682</td><td>0.01935</td><td>0.016715</td><td>0.024345</td><td>0.021311</td><td>0.041421</td></tr><tr><td>Obliquity<br>(degrees)</td><td>Same as<br>clim_pdControl</td><td>24.105</td><td>24.231</td><td>23.441</td><td>22.391</td><td>22.493</td><td>22.404</td></tr><tr><td>Perihelion – 180</td><td>Same as<br>clim_pdControl</td><td>0.87</td><td>303.03</td><td>102.7</td><td>80.09</td><td>174.82</td><td>110.88</td></tr><tr><td>Solar constant<br>(W/m2)</td><td>Same as<br>clim_pdControl</td><td>Same as<br>clim_pdControl</td><td>Same as<br>clim_pdControl</td><td>1365.6537</td><td>1365.6537</td><td>1365.6537</td><td>1361.20</td></tr><tr><td colspan="8">Gaz concentration</td></tr><tr><td>Carbon dioxide<br>(ppm)</td><td>Same as<br>clim_pdControl</td><td>264</td><td>287</td><td>205</td><td>230</td><td>195</td><td>274</td></tr><tr><td>Methane<br>(ppb)</td><td>Same as<br>clim_pdControl</td><td>597</td><td>791</td><td>500</td><td>450</td><td>450</td><td>505</td></tr><tr><td>Nitrous oxide<br>(ppb)</td><td>Same as<br>clim_pdControl</td><td>262</td><td>275</td><td>260</td><td>230</td><td>217</td><td>251</td></tr><tr><td>SST</td><td>Same as<br>clim_pdControl</td><td colspan="7">Simulated tsol_oce from pre-existing simulations (corrected for the models' systematic bias)</td></tr><tr><td>Geography</td><td>Same as<br>clim_pdControl</td><td colspan="7">Same as pre-existing simulations</td></tr><tr><td>Vegetation</td><td>Same as<br>clim_pdControl</td><td colspan="7">Same as pre-existing simulations</td></tr></tbody></table>

### 2.2.4 A subset of key paleoclimate variables

To characterize the large-scale climate over the area, we worked on the mean annual cycle of the four grid cells containing EH2. From simulations *Ctrl*, *midH*, *earlyH*, *midMIS3*, *lateMIS4*, *midMIS4* and *MIS5d* we extracted nine output variables. We choose to focus on variables that are likely to directly or

<footnote>* The term "orbital parameters" refers to variations in the eccentricity of the Earth and longitude of perihelion as well as changes in its axial inclination (obliquity).</footnote>
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

indirectly influence the landscape and/or biotic environment. Based on a review of the ecological literature, we selected: tsol (temperature at surface, in C°) (Gillooly et al., 2001; Yom-Tov and Geffen, 2006; Ebrahimi-Khusfi et al., 2020), precip (precipitations, in mm.month⁻¹) (Yom-Tov and Geffen, 2006; Alhajeri and Steppan, 2016; Ebrahimi-Khusfi et al., 2020), qsurf (specific humidity, in kg.kg⁻¹) (Hovenden et al., 2012; Alhajeri and Steppan, 2016), w10m (wind speed at 10 meters, in m.s⁻¹) (McNeil, 1991; Tanner et al., 1991; Chapman et al., 2011; Pellegrino et al., 2013), sols (solar radiation at surface, in W.m⁻²) (Monteith, 1972; Fyllas et al., 2017), drysoil Frac (fraction of visibly dry soil, in %) (Paz et al., 2015) and humtot (total soil moisture, in kg.m⁻²) (Paz et al., 2015). Two additional variables were computed. The diurnal temperature range tsol_ampl_day (in C°; 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 (in mm.day⁻¹; Martínez-Blancas and Martorell, 2020) computed from evapot (potential evaporation) and evap (evaporation) as: evapot - evap. We explore variations in monthly means (mean value of the variable) and monthly standard deviations (amplitude of seasonal variation). Correlation between these variables was tested and not found to be significant, suggesting that they all bring complementary information on environmental conditions over the regions

The association of paleoclimate simulations with the stratigraphic layers of EH2 is based on age
proximity. This step is performed for both DH1 and DH2 (the DH presented in section 2.1.2). As a
result, we obtain two hypothetic paleoclimate sequences corresponding to EH2 sequence.

## 2.3 Consistency analyses

To test the consistency between the paleoclimate simulations and paleoenvironmental inferences that have been made from the content of the EH2 site we use two kind of analyses: two-block Partial Least Squares (2B-pls) and pairwise correlation tests. For these analyses, the climate variables are the principal components provided by the principal component analysis described above. The different paleoenvironmental proxies considered are δ¹⁸O and δ¹³C means, minimums and maximums, and MAP, as well as percentages of represented habitats indicated by the THI.
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

265 The 2B-pls tests the global covariance between climate variables and paleoenvironmental variables. This multivariate method explores patterns of covariation between two sets of variables, i.e. two blocks (Sampson et al., 1989; Streissguth et al., 1993). Axes of maximum covariance between the two blocks are generated, thus reducing data dimensionality. A coefficient (r-PLS) is computed and represents the strength of covariation. The r-PLS is in the range of (0,1). The closer the r-PLS is to one, the stronger the covariation. P-values indicating the statistical significance of r-PLS were calculated based on 1000 permutations against the null hypothesis (absence of covariation between the two sets of variables). Then, to refine our results, we performed pairwise correlation tests.

270

275

In order to easily visualize the climate proximity/differences between EH2 layers we used a principal component analysis (Jolliffe and Cadima, 2016). This analysis allows one to find new uncorrelated variables that successively maximize variance (the principal components). They are computed from the eigenvectors and eigenvalues of the covariance/correlation matrix (correlation matrix in our case, as all the variables are standardized). Then, the data can be visualized along the leading principal components, which maximizes the amount of information displayed. It allows us to visualize EH2 layers into a climate-space, meaning that the proximity of stratigraphic layers in the plot simply represents the similarity between the climate of the layers. For each climate variable, we considered both the annual mean and seasonal variation, estimated by the annual standard deviation.

280

## 3 Results

### 3.1. Simulated climate changes

285 Changes in atmospheric circulation and pressure over the sequence are presented in Fig 5, and plots of monthly precipitations and temperatures for the region of EH2 are available in Fig 6. The similarity/dissimilarity between the climate of each period and the current climate is showed on similarity maps on Fig 7.

290

The dynamic over the region differ between winter and summer. Winter is marked by a vertical pressure gradient and dominated by a zonal circulation from west to east. In summer, the pressure gradient is
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

more horizontally-oriented, with the installation of a depression on the North African continent and the northward drift of high pressures on the Atlantic Ocean (Fig 5). All over the sequence, these processes shift slightly, especially between lateMIS4, midMIS4, midMIS3 (period from ~66ka to ~40ka) and Ctrl, midH, earlyH (period from ~9ka until today) (Fig 5 and 7), at the exception of MIS5d (~115ka) which conditions are relatively close from current ones on Fig 7. 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 (Fig 6).

From ~115ka (*MIS5d*) until ~40ka (*midMIS3*), the climate was colder than today. From ~66ka (*midMIS4*) to ~40ka (*midMIS3*), there is also more precipitations in winter (Fig 6). These conditions share similarities with what can be found in slightly higher latitudes today (Fig 7), it is concomitant with the pressure along the North Africa (and Morrocan) coast (Fig 5). In ~115ka (*MIS5d*), conditions are also cold, but as dry as in the Holocene (Fig 5). They are not related to changes in the atmosphere dynamic (Fig 5) but rather to changes in insolation, with current analogs limited to North Africa (Fig 7). Starting from ~9ka (*earlyH*) the seasonal temperature variation is enhanced. Conditions in winter and spring are similar to the current climate (Fig 7), but with a warmer autumn and a much warmer summer (Fig 6). At ~6ka (*midH*), climate conditions are close from today (Fig 7), but with a slightly more important seasonal temperature variation (Fig 6).

295

300

305
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. 5.** Maps of wind speed and direction at 850hPa (represented by arrows, length is proportional to wind speed) and corrected pressure at sea level (unit: hPa) per season simulated with *midH*, *earlyH*, *midMIS3*, *lateMIS4*, *midMIS4* and *MIS5d*. EH2 cave location is represented by a star in the upper left panel. DJF: December, January, February (winter); JJA: June, July, August (summer).</figcaption></figure>

310
<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>**Fig. 6.** Graph of monthly variations in mean temperature (unit: °C) and precipitation (unit: mm/month) averaged from the four grid cells containing EH2 cave in the different simulations. Interannual variation over the averaged 30 years is visualized by quartiles.</figcaption></figure>

315

## 3.2 Consistency between paleoclimate simulations and paleoenvironmental inferences

### 3.2.1 Association of paleoclimate simulations and stratigraphic layers

320

The two hypothetic paleoclimate sequences corresponding respectively to DH1 and DH2 are presented in Fig 8. The two principal components analyses are shown in Fig 9. To visualize how climate variables are associated to structure the climate-space of these two principal component analyses, biplots are available in SM Fig 2.
<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

325

**Fig. 7.** Maps presenting the similarity between past periods climate and current climate. Each map results from a distance matrix computed between climate variables on EH2 in the *Ctrl* simulation and the values of each cells of the past period simulation. Blue cells indicate localities where climate is at least 98% similar to current climate on EH2 locality. EH2 cave location is represented by a star in the upper left panel.
<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /><figcaption>**Fig. 8.** Climate variation over the EH2 sequence according to DH1 (Dating Hypothesis 1; upper), and DH2 (Dating Hypothesis 2; lower). Climate variables are centered and reduced and share a common scale to ease the reading and better visualize covariation between climate variables. Maximum/minimum refers to the maximum/minimum value for all simulations. “L” is the abbreviation for Layer, “Act” for Actual, “m” for mean and “sd” for standard deviation.</figcaption></figure>

330

335
<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /><figcaption>Fig. 9. Principal component analyses performed on climate variables according to DH1 (Dating Hypothesis 1; upper), and DH2 (Dating Hypothesis 2; lower). The contributions of the (standardized) climate variables to the two first principal components are presented (variables above the horizontal red line contribute significantly). “L” is the abbreviation for Layer and “Act” for Actual.</figcaption></figure>

340
<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

Based on DH1, our results indicate four major climate transitions (Fig 8). The first occurs between L8 and L5. In L5 the climate is wetter and colder with more precipitation, more soil humidity, increased wind speed, less hydric stress and less portions of dry soil. Humidity, precipitations and wind speed also show an important seasonal variability. The second transition, less marked, is between L5 and L4a. Climate in L4a is rather similar to the climate in L5, but precipitations and humidity have increased. Seasonal variations are globally more pronounced. The third transition is between L3 and L2. The climate changes drastically between these two periods, with hotter and drier conditions in L2 corresponding to the end of the last deglaciation. Temperature, solar radiation, water stress and soil dryness all increase, coupled with a decrease in precipitation, soil moisture and wind speed. All these changes are consistent with aridification or desertification from L2. Surprisingly, however, specific humidity is enhanced, which is partly consistent with the decrease in precipitation and wind speed at this coastal location, and is concomitant with large changes of the atmospheric circulation patterns over this region from L2 (Fig 5). The last climatic transition is more subtle and occurs between L1 and Act. The environment in Act seems closer to the one in L8 with more seasonal variability in temperature, solar radiation and water stress (Fig 8).

Overall, there is an alternation of two main climate types. This partition is confirmed by the principal component analysis shown in Fig 9. The first principal component explains 65.85% of the observed variance (of the standardized variables) and splits the layers of EH2 into two climate regimes. The first regroups L3, L4a and L5 and is defined by humid and windy conditions with a high seasonal variability. The second 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 and abrupt climatic transitions (Fig 8). The first happened between L5 and L4a. In L4a, the climate is much windier and temperatures are colder with a substantial increase of diurnal temperature range. Precipitations and soil moisture increase importantly while hydric stress decreases. The climate also presents an overall higher seasonal variability. These tendencies persist in L3. The second transition is between L3 and L2. The soil is drier and the hydric stress increases

345

350

355

360

365
<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

370 greatly as well as solar radiation and temperature. Precipitation and soil moisture are less important. Conditions in L1 are close from those in L2, a bit colder with a less marked water stress. The last climate transition, smoother than the previous ones, occurs between L1 and Act and is mainly marked by a global decrease of seasonal variations.

375 Results associated with DH2 suggest that three types of climate 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 for DH1, this partition is supported by the principal component analysis results presented in Fig 9. The first principal component, which explains 52.10% of the observed variance, separates the group containing L8, L7, L6 and L5 with the one composed of L3 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 present a hot environment associated with an important water stress.

### 3.2.2 Consistency analyses

385 Results of 2B-pls and pairwise correlations between climate and paleoenvironmental variables are presented in Fig 10. Regarding DH1, 2B-pls shows that there is a statistically significant covariation between THI values and all climate variables, and between isotope data and the second principal component (climate variables contributing to the second principal component are indicated in Fig 9). Specifically, this second principal component is positively correlated with δ¹⁸O mean, maximum and minimum values and δ¹³C maximum values. Conversely, DH2 presents no statistically significant result for 2B-pls nor correlations, meaning that for DH2 paleoclimate inferences and paleoenvironmental proxies are not consistent.
<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

DH1

Isotopes

Correlation tests

<figure><img src="image_55.png" /></figure>

2B-pls

■ Climate variables

P-value = 0.5
r-PLS = 0.5

■ PC1

P-value = 0.9
r-PLS = 0.08

■ PC2

<figure><img src="image_56.png" /></figure>

THI

Correlation tests

<figure><img src="image_57.png" /></figure>

2B-pls

■ Climate variables

<figure><img src="image_58.png" /></figure>

■ PC1

P-value = 0.2
r-PLS = 0.8

■ PC2

P-value = 0.7
r-PLS = 0.5

DH2

Isotopes

Correlation tests

<figure><img src="image_59.png" /></figure>

2B-pls

■ Climate variables

P-value = 0.3
r-PLS = 0.5

■ PC1

P-value = 0.2
r-PLS = 0.5

■ PC2

P-value = 0.3
r-PLS = 0.5

THI

Correlation tests

<figure><img src="image_60.png" /></figure>

2B-pls

■ Climate variables

P-value = 0.4
r-PLS = 0.5

■ PC1

P-value = 0.5
r-PLS = 0.5

■ PC2

P-value = 0.3
r-PLS = 0.6

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /><figcaption>Fig. 10. Graphical representation of the 2B-pls and correlations between isotopes (Jeffrey, 2016) and THI values (Stoetzel et al., 2014) and climate variables according to DH1 (Dating Hypothesis 1; upper panel), and DH2 (Dating Hypothesis 2; lower panel). Crosses indicate cases where correlation is not significant (p-value > 0.05) and colors represent the strength of the correlations. “L” is the abbreviation for Layer, “m” for mean and “sd” for standard deviation. δ¹⁸O: mean δ¹⁸O values in *Meriones* teeth (from Jeffrey (2016)); δ¹³C: mean δ¹³C values in *Meriones* teeth (from Jeffrey (2016)); MAP: Mean Annual Precipitations (from Jeffrey (2016)); Forest, Bush, Steppe, Wetland, Rocky: relative % of representation according to the THI (from Stoetzel *et al.* (2014) and Jeffrey (2016)).</figcaption></figure>
<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

## 4 Discussion

Despite the differences in scale and resolution between climate and paleoenvironmental data, we find statistically significant and meaningful results. The climate changes considered seem to be large enough for a consistency to be detected between the climate and environmental data. First, we discuss the paleoclimate changes described by simulations over the period, and the underlying dynamical processes. Then, we address the contribution of climate simulations to the discussion of chronological and paleoenvironmental discrepancies of EH2.

405

410

### 4.1 Paleoclimate variation and underlying forcings

Paleoclimate simulations allow us to discuss several important climate changes at EH2 over the Late Pleistocene to mid-Holocene period. These paleoclimate variations may result from mixed influences from global and regional dynamical processes.

The variation of radiative-related variables (as sols, tsol or tsol_ampl_day) seems to depend more of global processes related to large-scale climate changes. They clearly separate interglacial climate (Ctrl, midH, earlyH) from glacial climate (midMIS3, lateMIS4, midMIS4, MIS5d). The warm/cold differences between these two periods is due to the size of the ice-sheet and variation in gaz concentration. For example, mean insolation and temperature significantly increase from early Holocene, along with higher gaz concentration and the retreat of the ice sheet. The amplitude of the seasonal variability of insolation and temperature also increases at this period, along with variations in obliquity and precession of the Earth orbit. Note that MIS5d is a peak of glacial sub-stage where conditions where quite mild, which explains the proximity of its climate with current ones.

415

420

The variation of the humidity-related variables (as precip, hydric_stress or drysoil-frac) seems to be mainly explained by a translation of the regional atmospheric dynamic processes over the sequence. For example, mean precipitations and soil moisture are higher from mid-MIS4 to mid-MIS3. At this time, the wet and cold westerly winds descend further south and blow across the area of EH2. This is explained by the position of the Açores High: the proximity of the high pressures creates a stronger pressure

425
<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

gradient, thus favoring a stronger zonal circulation, which brings cooler and more humid air on the
North African coast. Therefore, for this region, the dynamic effect and changes in moisture advection
over the region have a significant impact on whether a warmer climate would lead to increased moisture
content and precipitation.

430

## 4.2 Paleoclimate simulations and chronostratigraphy

The climate sequence varies significantly depending on the DH it is based on. Based on DH1, there is an alternation of semi-arid and temperate climates. The climate succession based on DH2 is quite different, with the presence of three main climate types and rapid transitions between them. The climate sequences based on DH1 and DH2 are not equally congruent with paleoenvironmental proxies from the literature. Indeed, several statistically significant correlations and covariation are found between the climate sequence based on DH1 and the paleoenvironmental inferences, while none are found for the one based on DH2. Thus, our results suggest that, with respect to EH2, combined US-ESR dating may be more reliable than OSL dating. As this dating process relies on quartz grains and that their chronology and origin is difficult to establish in the context of karstic coastal caves (as discussed in Ben Arous et al. (2020a)), OSL ages might have been overestimated (i.e. the age of the layers may be younger than estimated by OSL). Moreover, OSL dates can in some cases be internally inconsistent, meaning that they can have a number of reversals. For instance, it is known that water content influences the determinations. Considering the location of the cave on the coastline, inundations related to rising in sea level or very high tides may have occurred.

435

440

445

## 4.3 Paleoclimate simulations and paleoenvironmental inferences

Since the climate sequence based on DH1 is the only one displaying significant results, we focus on it
in the following discussion. Isotope fractions are mainly correlated to seasonal variation in water stress
(hyd_stress), insolation (sols) and temperature (tsol, tsol_ampl_day). This result is expected as δ¹⁸O is
an indicator of aridity (Longinelli and Selmo, 2003; Blumenthal et al., 2017; Blumenthal, 2019) and as

450
<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

455 $\delta^{13}$C is related to seasonal variation of temperature and water stress (O'Leary, 1988; Lin, 2013; Smiley et al., 2016). On the contrary to isotopes, the THI is only globally related to climate variations (only the 2B-pls performed on all climate variables is significant). It is not surprising either, as isotopes reveal fluctuations in particular variables that are temperature- and precipitation-related provided by a limited number of individuals of a single species, while THI is an estimate of the global type of the environment provided by the whole microvertebrate communities. Thus, they do not deliver information at the same resolution.

460

Another explanation of the improved correlations between climate and isotope fractions may be that biologically-derived proxies (as the THI) are more complex functions of physical drivers than the isotope signal. Isotopes directly reflect to the magnitude of seasonal variation in insolation, water stress, temperature and diurnal temperature range, because these variables condition the presence of essential elements for plants survival (e.g. sunlight, water in the soil). Thus, they are directly related to the type of vegetation. On the other hand, the THI relies on ecological preferences of species. Altogether, the variables of the THI give information about the proportion of biomes (e.g. forest, bush, steppe), and thus the spatial distribution and density of the vegetation. Consequently, its relationship to climate is more indirect than for isotopes.

465

470

The large differences between the climate simulations and the fact that they provide a physical consistent view of the relationships 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 cases, isotope surveys and MAP from Jeffrey (2016) indicate more humid conditions with more important precipitation than on other layers. 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 drier climate than on other layers (Stoetzel, 2009; Stoetzel et al., 2014). Large mammals would also support this last hypothesis, with an increase in the representation of gazelles and alcelaphines, and a decrease in the representation of bovines in both layers (Stoetzel et al., 2012a, 2014). Unfortunately, no combined US-ESR ages are available for L7 to date. Considering that climate conditions associate to L8 display less precipitation than currently, we could hypothesize a

475

480
<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

similar climate for L7. In that case, this would support inferences from species presence. Nevertheless, we cannot exclude that a microclimatic event could have induced particular climatic conditions on L7. Concerning L5, the climate described by the sequence based on DH1 agrees with isotope surveys (Jeffrey, 2016). These conclusions are supported by the 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.

An important difference is noticed between the climate sequence and the THI on L1. Paleoclimate simulations indicate on a quite dry climate relatively to other layers, in contradiction with the composition of small and large mammal communities that suggest an expansion of forests and wetlands (Stoetzel *et al.*, 2014). This difference could be explained by the location of EH2: the cave is located at the interface of varied climate influences, as described previously. Because global climate models describe general climate characteristics, a local climate phenomenon could have generated a wetter environment in the surroundings of EH2.

485

490

495

## 5 Conclusions

Considering together paleoclimate simulations and paleoenvironmental inferences allows us to provide new insights into the chronostratigraphy and paleoenvironmental reconstruction of El Harhoura 2 cave. We find that the climate sequence based on combined US-ESR ages is more consistent with paleoenvironmental inferences than the climate sequence based on OSL ages. We also evidence that, overall, isotope-based paleoenvironmental inferences are more congruent with the paleoclimate sequence than species-based inferences. But, above all, we highlight the difference in scale between the information provided by each of these paleoenvironmental proxies. This study demonstrates that the combination of different sources of environmental data and climate simulations has a great potential for refining the paleoenvironmental and chronological context of archeological and paleontological sites. Even thought, its applicability on periods marked by less important climate changes remains to be tested. Although, our approach may concern a limited number of well-studied sites. However, with more

500

505
<figure><img src="image_71.png" /></figure>

510

powerful statistical tools, it could be extended to other sites whose context is less referenced, and
perhaps even help establish it.

Author contribution

Conceptualization: Léa Terray, Pascale Braconnot, Raphaël Cornette, Emmanuelle Stoetzel.

Formal analyses: Léa Terray.

515

Methodology: Léa Terray, Pascale Braconnot, Masa Kageyama.

Funding acquisition: Léa Terray.

Writing – original draft preparation: Léa Terray.

Writing – review & editing: Pascale Braconnot, Raphaël Cornette, Emmanuelle Stoetzel, Eslem Ben Arous.

520

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgments

525

The microvertebrate remains from El Harhoura 2 cave were recovered and studied within the framework
of the El Harhoura-Témara Archaeological Mission (dir. R. Nespoulet and M.A. El Hajraoui), under the
administrative supervision of the Institut National des Sciences de l'Archéologie et du Patrimoine
(Rabat, Morocco) and financial support from the Ministère des Affaires Etrangères et du Développement
International (France) and the Ministère de la Culture (Maroc). The authors would like to thank Marie

530

Sicard for sharing the *lig115* climate simulation, part of her PhD project. Climate simulations were
carried out on the Jean Zay supercomputer at IDRIS and the computing time was provided by the project
gen2212. The authors would like to thank Julia Lee-Thorp (School of Archeology, University of Oxford)
for sharing her expertise on stable light isotopes and for her suggestions on the manuscript. We are
<figure><img src="image_72.png" /></figure>

535

grateful to Pascal Terray (LOCEAN, IRD) for his comments on the manuscript and advice. This research
was funded by the Université 463 Paris Descartes, the Ecole Doctorale FIRE - Programme Bettencourt.

All operations on NetCDF files (the standard file format for the outputs of the IPSL models) were
performed using CDO (Climate Data Operators) (Schulzweida, 2019). Maps, plots and analyses were
produced using the R free software (R Development Core Team, 2018) and the libraries *ncdf4*
(http://dwpierce.com/software), *raster* (Hijmans and van Etten, 2012), *FactoMineR* (Lê et al., 2008),

540

corrplot (Wei and Simko, 2021) and ggplot2 (Wickham, 2015).
<figure><img src="image_73.png" /></figure>

References

Adler, R., Sapiano, M., Huffman, G., Wang, J.-J., Gu, G., Bolvin, D., Chiu, L., Schneider, U., Becker, A., Nelkin, E., Xie, P., Ferraro, R., and Shin, D.-B.: The Global Precipitation Climatology Project (GPCP) Monthly Analysis (New Version 2.3) and a Review of 2017 Global Precipitation, Atmosphere, 9, 138, https://doi.org/10.3390/atmos9040138, 2018.

545

Alhajeri, B. H. and Steppan, S. J.: Association between climate and body size in rodents: A phylogenetic test of Bergmann’s rule, Mamm. Biol., 81, 219–225, https://doi.org/10.1016/j.mambio.2015.12.001, 2016.

Avery, D. M.: Pleistocene micromammals from Wonderwerk Cave, South Africa: practical issues, J. Archaeol. Sci., 34, 613–625, https://doi.org/10.1016/j.jas.2006.07.001, 2007.

550

Ben Arous, E., Falguères, C., Tombret, O., El Hajraoui, M. A., and Nespoulet, R.: Combined US-ESR dating of fossil teeth from El Harhoura 2 cave (Morocco): New data about the end of the MSA in Temara region, Quat. Int., 556, 88–95, https://doi.org/10.1016/j.quaint.2019.02.029, 2020a.

555

Ben Arous, E., Falguères, C., Nespoulet, R., and El Hajraoui, M. A.: Review of chronological data from the Rabat-Temara caves (Morocco): Implications for understanding human occupation in north-west Africa during the Late Pleistocene, in: Not just a corridor. Human occupation of the Nile Valley and neighbouring regions between 75,000 and 15,000 years ago, edited by: Leplongeon, A., Goder-Goldberger, M., and Pleurdeau, D., Paris, 177–201, 2020b.

560

Berrisford, P., Dee, D., Poli, P., Brugge, R., Fielding, M., Fuentes, M., Kållberg, P., Kobayashi, S., Uppala, S., and Simmons, A.: The ERA-Interim archive Version 2.0, ERA Report Series, 2011.

Blome, M. W., Cohen, A. S., Tryon, C. A., Brooks, A. S., and Russell, J.: The environmental context for the origins of modern human diversity: A synthesis of regional variability in African climate 150,000–30,000 years ago, J. Hum. Evol., 62, 563–592, https://doi.org/10.1016/j.jhevol.2012.01.011, 2012.

565

Blumenthal, S. A.: Isotopic records of climate seasonality in equid teeth, Geochim. Cosmochim. Acta, 20, 2019.

570

Blumenthal, S. A., Levin, N. E., Brown, F. H., Brugal, J.-P., Chritz, K. L., Harris, J. M., Jehle, G. E., and Cerling, T. E.: Aridity and hominin environments, PNAS, 6, 2017.

Boucher, O., Denvil, S., Levavasseur, G., Cozic, A., Caubel, A., Foujols, M.-A., Meurdesoif, Y., Cadule, P., Devilliers, M., Ghattas, J., Lebas, N., Lurton, T., Mellul, L., Musat, I., Mignot, J., and Cheruy, F.: IPSL IPSL-CM6A-LR model output prepared for CMIP6 CMIP amip, https://doi.org/10.22033/ESGF/CMIP6.5113, 2018.

575

Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., Bekki, S., Bonnet, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Caubel, A., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., D'Andrea, F., Davini, P., Lavergne, C., Denvil, S., Deshayes, J., Devilliers, M., Ducharne, A., Dufresne, J., Dupont, E., Éthé, C., Fairhead, L., Falletti, L., Flavoni, S., Foujols, M., Gardoll, S., Gastineau, G., Ghattas, J., Grandpeix, J., Guenet, B., Guez, L., E., Guilyardi, E., Guimberteau, M., Hauglustaine, D., Hourdin, F., Idelkadi, A., Joussaume, S., Kageyama, M., Khodri, M., Krinner, G., Lebas, N., Levavasseur, G., Lévy, C., Li, L., Lott, F., Lurton, T., Luyssaert, S., Madec, G., Madeleine, J., Maignan, F., Marchand, M., Marti, O., Mellul, L., Meurdesoif, Y., Mignot, J., Musat, I., Ottlé, C., Peylin, P., Planton, Y., Polcher, J., Rio, C., Rochetin, N., Rousset, C., Sepulchre, P., Sima, A.,

580

585
<figure><img src="image_74.png" /></figure>

Swingedouw, D., Thiéblemont, R., Traore, A. K., Vancoppenolle, M., Vial, J., Vialard, J., Viovy, N., and Vuichard, N.: Presentation and Evaluation of the IPSL-CM6A-LR Climate Model, J. Adv. Model. Earth Syst., 12, https://doi.org/10.1029/2019MS002010, 2020.

590 Braconnot, P., Harrison, S. P., Kageyama, M., Bartlein, P. J., Masson-Delmotte, V., Abe-Ouchi, A., Otto-Bliesner, B., and Zhao, Y.: Evaluation of climate models using palaeoclimatic data, Nat. Clim. Change, 417–424, 8, 2012.

595 Braconnot, P., Albani, S., Balkanski, Y., Cozic, A., Kageyama, M., Sima, A., Marti, O., and Peterschmitt, J.-Y.: Impact of dust in PMIP-CMIP6 mid-Holocene simulations with the IPSL model, Clim. Past, 17, 1091–1117, https://doi.org/10.5194/cp-17-1091-2021, 2021.

Carto, S. L., Weaver, A. J., Hetherington, R., Lam, Y., and Wiebe, E. C.: Out of Africa and into an ice age: on the role of global climate change in the late Pleistocene migration of early modern humans out of Africa, J. Hum. Evol., 56, 139–151, https://doi.org/10.1016/j.jhevol.2008.09.004, 2009.

600 Chapman, J. W., Klaassen, R. H. G., Drake, V. A., Fossette, S., Hays, G. C., Metcalfe, J. D., Reynolds, A. M., Reynolds, D. R., and Alerstam, T.: Animal Orientation Strategies for Movement in Flows, Curr. Biol., 21, R861–R870, https://doi.org/10.1016/j.cub.2011.08.014, 2011.

605 Comay, O. and Dayan, T.: From micromammals to paleoenvironments, Archaeol. Anthropol. Sci., 10, 2159–2171, https://doi.org/10.1007/s12520-018-0608-8, 2018.

Cornette, R., Stoetzel, E., Baylac, M., Moulin, S., Hutterer, R., Nespoulet, R., El Hajraoui, M. A., Denys, C., and Herrel, A.: Shrews of the genus Crocidura from El Harhoura 2 (Témara, Morocco): The contribution of broken specimens to the understanding of Late Pleistocene–Holocene palaeoenvironments in North Africa, Palaeogeogr. Palaeoclimatol. Palaeoecol., 436, 1–8, https://doi.org/10.1016/j.palaeo.2015.06.020, 2015.

610
Couvreur, T. L. P., Dauby, G., Blach-Overgaard, A., Deblauwe, V., Dessein, S., Droissart, V., Hardy, O. J., Harris, D. J., Janssens, S. B., Ley, A. C., Mackinder, B. A., Sonké, B., Sosef, M. S. M., Stévart, T., Svenning, J., Wieringa, J. J., Faye, A., Missoup, A. D., Tolley, K. A., Nicolas, V., Ntie, S., Fluteau, F., Robin, C., Guillocheau, F., Barboni, D., and Sepulchre, P.: Tectonics, climate and the diversification of the tropical African terrestrial flora and fauna, Biol. Rev., 96, 16–51, https://doi.org/10.1111/brv.12644, 2020.

615

deMenocal, P. B.: Plio-Pleistocene African Climate, Science, 270, 53–59, https://doi.org/10.1126/science.270.5233.53, 1995.

620

deMenocal, P. B.: African climate change and faunal evolution during the Pliocene– Pleistocene, Earth Planet. Sci. Lett., 220, 3–24, https://doi.org/10.1016/S0012- 821X(04)00003-2, 2004.

Denys, C., Stoetzel, E., Andrews, P., Bailon, S., Rihane, A., Huchet, J. B., Fernandez-Jalvo, Y., and Laroulandie, V.: Taphonomy of Small Predators multi-taxa accumulations: palaeoecological implications, Hist. Biol., 30, 868–881, https://doi.org/10.1080/08912963.2017.1347647, 2018.

625

Drake, N. A., Blench, R. M., Armitage, S. J., Bristow, C. S., and White, K. H.: Ancient watercourses and biogeography of the Sahara explain the peopling of the desert, Proc. Natl. Acad. Sci., 108, 458–462, https://doi.org/10.1073/pnas.1012231108, 2011.

Drake, N. A., Breeze, P., and Parker, A.: Palaeoclimate in the Saharan and Arabian Deserts during the Middle Palaeolithic and the potential for hominin dispersals, Quat. Int., 300, 48–61, https://doi.org/10.1016/j.quaint.2012.12.018, 2013.

630
<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>

Dufresne, J.-L., Foujols, M.-A., Denvil, S., Caubel, A., Marti, O., Aumont, O., Balkanski, Y., Bekki, S., Bellenger, H., Benshila, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., de Noblet, N., Duvel, J.-P., Ethé, C., Fairhead, L., Fichefet, T., Flavoni, S., Friedlingstein, P., Grandpeix, J.-Y., Guez, L., Guilyardi, E., Hauglustaine, D., Hourdin, F., Idelkadi, A., Ghattas, J., Joussaume, S., Kageyama, M., Krinner, G., Labetoulle, S., Lahellec, A., Lefebvre, M.-P., Lefevre, F., Levy, C., Li, Z. X., Lloyd, J., Lott, F., Madec, G., Mancip, M., Marchand, M., Masson, S., Meurdesoif, Y., Mignot, J., Musat, I., Parouty, S., Polcher, J., Rio, C., Schulz, M., Swingedouw, D., Szopa, S., Talandier, C., Terray, P., Viovy, N., and Vuichard, N.: Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5, Clim. Dyn., 40, 2123–2165, https://doi.org/10.1007/s00382-012-1636-1, 2013.

635

640

Duplessy, J.-C. and Ramstein, G.: Paléonclimatologie, tome II : enquêter sur les climats anciens., CNRS Editions., 2013.

645

Ebrahimi-Khusfi, Z., Mirakbari, M., and Khosroshahi, M.: Vegetation response to changes in temperature, rainfall, and dust in arid environments, Environ. Monit. Assess., 192, 691, https://doi.org/10.1007/s10661-020-08644-0, 2020.

650

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958, https://doi.org/10.5194/gmd-9-1937-2016, 2016.

Fernandez-Jalvo, Y., Denys, C., Andrews, P., Williams, T., Dauphin, Y., and Humphrey, L.: Taphonomy and palaeoecology of Olduvai Bed-I (Pleistocene, Tanzania), J. Hum. Evol., 34, 137–172, 1998.

655

Fyllas, N. M., Bentley, L. P., Shenkin, A., Asner, G. P., Atkin, O. K., Díaz, S., Enquist, B. J., Farfan-Rios, W., Gloor, E., Guerrieri, R., Huasco, W. H., Ishida, Y., Martin, R. E., Meir, P., Phillips, O., Salinas, N., Silman, M., Weerasinghe, L. K., Zaragoza-Castells, J., and Malhi, Y.: Solar radiation and functional traits explain the decline of forest primary productivity along a tropical elevation gradient, Ecol. Lett., 20, 730–740, https://doi.org/10.1111/ele.12771, 2017.

660

Ghawar, W., Zaätour, W., Chlif, S., Bettaieb, J., Chelghaf, B., Snoussi, M.-A., and Salah, A. B.: Spatiotemporal dispersal of Meriones shawi estimated by radio-telemetry, Int. J. Multidiscip. Res. Dev., 2, 211–216, 2015.

Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M., and Charnov, E. L.: Effects of Size and Temperature on Metabolic Rate, Science, 293, 2248–2251, 2001.

665

Harrison, S. P., Bartlein, P. J., Izumi, K., Li, G., Annan, J., Hargreaves, J., Braconnot, P., and Kageyama, M.: Evaluation of CMIP5 palaeo-simulations to improve climate projections, Nat. Clim. Change, 5, 735–743, https://doi.org/10.1038/nclimate2649, 2015.

Haywood, A. M., Dowsett, H. J., Dolan, A. M., Rowley, D., Abe-Ouchi, A., Otto-Bliesner, B., Chandler, M. A., Hunter, S. J., Lunt, D. J., Pound, M., and Salzmann, U.: The Pliocene Model Intercomparison Project (PlioMIP) Phase 2: scientific objectives and experimental design, Clim. Past, 12, 663–675, https://doi.org/10.5194/cp-12-663-2016, 2016.

Hijmans, R. J. and van Etten, J.: raster: Geographic analysis and modeling with raster data. R package version 2.0-12., 2012.

670

Hooghiemstra, H., Stalling, H., Agwu, C. O. C., and Dupont, L. M.: Vegetational and climatic changes at the northern fringe of the sahara 250,000–5000 years BP: evidence from 4 marine

675
<figure><img src="image_77.png" /></figure>

pollen records located between Portugal and the Canary Islands, Rev. Palaeobot. Palynol., 74,
1–53, https://doi.org/10.1016/0034-6667(92)90137-6, 1992.

680

Hourdin, F., Rio, C., Grandpeix, J., Madeleine, J., Cheruy, F., Rochetin, N., Jam, A., Musat, I., Idelkadi, A., Fairhead, L., Foujols, M., Mellul, L., Traore, A., Dufresne, J., Boucher, O., Lefebvre, M., Millour, E., Vignon, E., Jouhaud, J., Diallo, F. B., Lott, F., Gastineau, G., Caubel, A., Meurdesoif, Y., and Ghattas, J.: LMDZ6A: The Atmospheric Component of the IPSL Climate Model With Improved and Better Tuned Physics, J. Adv. Model. Earth Syst., 12, https://doi.org/10.1029/2019MS001892, 2020.

685

Hovenden, M. J., Vander Schoor, J. K., and Osanai, Y.: Relative humidity has dramatic impacts on leaf morphology but little effect on stomatal index or density in Nothofagus cunninghamii (Nothofagaceae), Aust. J. Bot., 60, 700, https://doi.org/10.1071/BT12110, 2012.

690

Jacobs, Z., Roberts, R. G., Nespoulet, R., El Hajraoui, M. A., and Debénath, A.: Single-grain OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2, Morocco: Implications for Late Pleistocene human-environment interactions along the Atlantic coast of northwest Africa, J. Hum. Evol., 62, 377–394, https://doi.org/10.1016/j.jhevol.2011.12.001, 2012.

695

Janati-Idrissi, N., Falgueres, C., Nespoulet, R., El Hajraoui, M. A., Debénath, A., Bejjit, L., Bahain, J.-J., Michel, P., Garcia, T., Boudad, L., El Hammouti, K., and Oujaa, A.: Datation par ESR-U/th combinées de dents fossiles des grottes d'El Mnasra et d'El Harhoura 2, région de Rabat-Temara. Implications chronologiques sur le peuplement du Maroc atlantique au Pléistocène supérieur et son, Quaternaire, 23, 25–35, https://doi.org/10.4000/quaternaire.6127, 2012.

700

Jeffrey, A.: Exploring palaeoaridity using stable oxygen and carbon isotopes in small mammal teeth: a case study from two Late Pleistocene archaeological cave sites in Morocco, North Africa, 2016.

705

Jolliffe, I. T. and Cadima, J.: Principal component analysis: a review and recent developments, Philos. Trans. R. Soc. Math. Phys. Eng. Sci., 374, 20150202, https://doi.org/10.1098/rsta.2015.0202, 2016.

710

Kageyama, M., Braconnot, P., Bopp, L., Caubel, A., Foujols, M.-A., Guilyardi, E., Khodri, M., Lloyd, J., Lombard, F., Mariotti, V., Marti, O., Roy, T., and Woillez, M.-N.: Mid-Holocene and Last Glacial Maximum climate simulations with the IPSL model—part I: comparing IPSL_CM5A to IPSL_CM4, Clim. Dyn., 40, 2447–2468, https://doi.org/10.1007/s00382-012-1488-8, 2013.

715

Kutzbach, J. E. and Otto-Bliesner, B. L.: The Sensitivity of the African-Asian Monsoonal Climate to Orbital Parameter Changes for 9000 Years B.P. in a Low-Resolution General Circulation Model, J. Atmospheric Sci., 39, 1982.

Le Houérou, H. N.: Climate, flora and fauna changes in the Sahara over the past 500 million years, J. Arid Environ., 37, 619–647, https://doi.org/10.1006/jare.1997.0315, 1997.

Lê, S., Josse, J., and Husson, F.: FactoMineR: an R Package for multivariate analysis, J. Stat. Softw., 25, https://doi.org/10.18637/jss.v025.i01, 2008.

Lin, G.: Chapter 4: Research on stable isotope and carbon cycle (1st ed.), in: Stable Isotope Ecology, Beijing, 89–123, 2013.

720

Lionello, P., Malanotte, P., and Boscolo, R. (Eds.): Mediterranean Climate Variability, Elsevier., Elsevier, Amsterdam, 2006.
<figure><img src="image_78.png" /></figure>

Longinelli, A. and Selmo, E.: Isotopic composition of precipitation in Italy: a first overall map, J. Hydrol., 270, 75–88, https://doi.org/10.1016/S0022-1694(02)00281-0, 2003.

Marquer, L., Otto, T., Ben Arous, E., Stoetzel, E., Campmas, E., Zazzo, A., Tombret, O., Falgueres, C., El Hajraoui, M. A., and Nespoulet, R.: The first use of olives in Africa around 100,000 years ago, Nat. Plants, 8, 204–208, 2022.

725

Marti, O., Braconnot, P., Dufresne, J.-L., Bellier, J., Benshila, R., Bony, S., Brockmann, P., Cadule, P., Caubel, A., Codron, F., de Noblet, N., Denvil, S., Fairhead, L., Fichefet, T., Foujols, M.-A., Friedlingstein, P., Goosse, H., Grandpeix, J.-Y., Guilyardi, E., Hourdin, F., Idelkadi, A., Kageyama, M., Krinner, G., Lévy, C., Madec, G., Mignot, J., Musat, I., Swingedouw, D., and Talandier, C.: Key features of the IPSL ocean atmosphere model and its sensitivity to atmospheric resolution, Clim. Dyn., 34, 1–26, https://doi.org/10.1007/s00382-009-0640-6, 2010.

730

Martínez-Blancas, A. and Martorell, C.: Changes in niche differentiation and environmental filtering over a hydric stress gradient, J. Plant Ecol., 13, 185–194, https://doi.org/10.1093/jpe/rtz061, 2020.

735

Matthews, T.: Predators, prey and the palaeoenvironment, South Afr. J. Sci., 95, 22–24, 2000.

McNeil, J. N.: Behavioral Ecology of Pheromone-Mediated Communication in Moths and Its Importance in the Use of Pheromone Traps, Annu. Rev. Enthomology, 36, 407–30, 1991.

740

Michel, P., Campmas, É., Stoetzel, E., Nespoulet, R., Abdeljalil El Hajraoui, M., and Amani, F.: La macrofaune du Pléistocène supérieur d’El Harhoura 2 (Témara, Maroc) : données préliminaires, L’Anthropologie, 113, 283–312, https://doi.org/10.1016/j.anthro.2009.04.003, 2009.

Monteith, J. L.: Solar Radiation and Productivity in Tropical Ecosystems, J. Appl. Ecol., 9, 747, https://doi.org/10.2307/2401901, 1972.

745

Navarro, N., Lécuyer, C., Montuire, S., Langlois, C., and Martineau, F.: Oxygen isotope compositions of phosphate from arvicoline teeth and Quaternary climatic changes, Gigny, French Jura, Quat. Res., 62, 172–182, https://doi.org/10.1016/j.yqres.2004.06.001, 2004.

Nespoulet, R. and El Hajraoui, M. A.: Excavation report, 2012.

750

O’Leary, M. H.: Carbon Isotopes in Photosynthesis, BioScience, 38, 328–336, https://doi.org/10.2307/1310735, 1988.

Otto-Bliesner, B. L., Brady, E. C., Zhao, A., Brierley, C. M., Axford, Y., Capron, E., Govin, A., Hoffman, J. S., Isaacs, E., Kageyama, M., Scussolini, P., Tzedakis, P. C., Williams, C. J. R., Wolff, E., Abe-Ouchi, A., Braconnot, P., Ramos Buarque, S., Cao, J., de Vernal, A., Guarino, M. V., Guo, C., LeGrande, A. N., Lohmann, G., Meissner, K. J., Menviel, L., Morozova, P. A., Nisancioglu, K. H., O’ishi, R., Salas y Mélia, D., Shi, X., Sicard, M., Sime, L., Stepanek, C., Tomas, R., Volodin, E., Yeung, N. K. H., Zhang, Q., Zhang, Z., and Zheng, W.: Large-scale features of Last Interglacial climate: results from evaluating the <i>lig127k<i/> simulations for the Coupled Model Intercomparison Project (CMIP6)– Paleoclimate Modeling Intercomparison Project (PMIP4), Clim. Past, 17, 63–94, https://doi.org/10.5194/cp-17-63-2021, 2021.

755

Paz, H., Pineda-García, F., and Pinzón-Pérez, L. F.: Root depth and morphology in response to soil drought: comparing ecological groups along the secondary succession in a tropical dry forest, Oecologia, 179, 551–561, https://doi.org/10.1007/s00442-015-3359-6, 2015.

760

Pellegrino, A. C., Peñaflor, M. F. G. V., Nardi, C., Bezner-Kerr, W., Guglielmo, C. G., Bento,

765
<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

J. M. S., and McNeil, J. N.: Weather Forecasting by Insects: Modified Sexual Behaviour in Response to Atmospheric Pressure Changes, PLoS ONE, 8, e75004, https://doi.org/10.1371/journal.pone.0075004, 2013.

770 R Development Core Team: R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL: http://www.R-project.org., 2018.

775 Royer, A., Lécuyer, C., Montuire, S., Amiot, R., Legendre, S., Cuenca-Bescós, G., Jeannet, M., and Martineau, F.: What does the oxygen isotope composition of rodent teeth record?, Earth Planet. Sci. Lett., 361, 258–271, https://doi.org/10.1016/j.epsl.2012.09.058, 2013.

Sampson, P. D., Streissguth, A. P., and Bookstein, F. L.: Neurobehavioral Effects of Prenatal Alcohol: Part II. Partial Least Squares Analysis I, Neurotoxicol. Teratol., 11, 477–491, 1989.

780 Schmidt, G. A., Annan, J. D., Bartlein, P. J., Cook, B. I., Guilyardi, E., Hargreaves, J. C., Harrison, S. P., Kageyama, M., LeGrande, A. N., Konecky, B., Lovejoy, S., Mann, M. E., Masson-Delmotte, V., Risi, C., Thompson, D., Timmermann, A., Tremblay, L.-B., and Yiou, P.: Using palaeo-climate comparisons to constrain future projections in CMIP5, Clim. Past, 10, 221–250, https://doi.org/10.5194/cp-10-221-2014, 2014.

Schulzweida, U.: CDO User Guide (Version 1.9.8), 2019.

785 Sicard, M., Kageyama, M., Charbit, S., Braconnot, P., and Madeleine, J.-B.: An energy budget approach to understand the Arctic warming during the Last Interglacial, Clim. Past, 18, 607–629, https://doi.org/10.5194/cp-18-607-2022, 2022.

Smiley, T. M., Cotton, J. M., Badgley, C., and Cerling, T. E.: Small-mammal isotope ecology tracks climate and vegetation gradients across western North America, Oikos, 125, 1100–1109, https://doi.org/10.1111/oik.02722, 2016.

790 Sobrino, J. A. and Raissouni, N.: Toward remote sensing methods for land cover dynamic monitoring: Application to Morocco, Int. J. Remote Sens., 21, 353–366, https://doi.org/10.1080/014311600210876, 2000.

Stoetzel, E.: Les microvertébrés du site d'occupation humaine d'El Harhoura 2 (Pleistocene supérieur - Holocène, Maroc) : systématique, évolution, taphonomie et paléoécologie., 437, 2009.

795 Stoetzel, E., Marion, L., Nespoulet, R., El Hajraoui, M. A., and Denys, C.: Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Témara, Morocco), J. Hum. Evol., 60, 1–33, https://doi.org/10.1016/j.jhevol.2010.07.016, 2011.

800 Stoetzel, E., Bougariane, B., Campmas, E., Ouchaou, B., and Michel, P.: Chapitre V. Faunes et paléoenvironnements., in: Préhistoire de la région de Rabat-Témara., vol. 3, Rabat, 35–51, 2012a.

805 Stoetzel, E., Denys, C., Bailon, S., El Hajraoui, M. A., and Nespoulet, R.: Taphonomic Analysis of Amphibian and Squamate Remains from El Harhoura 2 (Rabat-Témara, Morocco): Contributions to Palaeoecological and Archaeological Interpretations: Taphonomic Study of Amphibian and Squamate Fossil Remains, Int. J. Osteoarchaeol., 22, 616–635, https://doi.org/10.1002/oa.1275, 2012b.

Stoetzel, E., Campmas, E., Michel, P., Bougariane, B., Ouchaou, B., Amani, F., El Hajraoui, M. A., and Nespoulet, R.: Context of modern human occupations in North Africa:

810 Contribution of the Témara caves data, Quat. Int., 320, 143–161, https://doi.org/10.1016/j.quaint.2013.05.017, 2014.
<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /></figure>

815

Stoetzel, E., Lalis, A., Nicolas, V., Aulagnier, S., Benazzou, T., Dauphin, Y., El Hajraoui, M. A., El Hassani, A., Fahd, S., Fekhaoui, M., Geigl, E.-M., Lapointe, F.-J., Leblois, R., Ohler, A., Nespoulet, R., and Denys, C.: Quaternary terrestrial microvertebrates from mediterranean northwestern Africa: State-of-the-art focused on recent multidisciplinary studies, Quat. Sci. Rev., 224, 105966, https://doi.org/10.1016/j.quascirev.2019.105966, 2019.

Streissguth, A. P., Bookstein, F. L., Sampson, P. D., and Barr, H. M.: The enduring effects of prenatal alcohol exposure on child development: Birth through seven years, a partial least squares solution., The University of Michigan Press., 1993.

820

Tanner, E. V. J., Kapos, V., and Healey, J. R.: Hurricane Effects on Forest Ecosystems in the Caribbean, Biotropica, 23, 513, https://doi.org/10.2307/2388274, 1991.

Tieszen, L. L.: Natural variations in the carbon isotope values of plants: Implications for archaeology, ecology, and paleoecology, J. Archaeol. Sci., 18, 227–248, https://doi.org/10.1016/0305-4403(91)90063-U, 1991.

825

Trauth, M. H., Larrasoaña, J. C., and Mudelsee, M.: Trends, rhythms and events in Plio-Pleistocene African climate, Quat. Sci. Rev., 28, 399–411, https://doi.org/10.1016/j.quascirev.2008.11.003, 2009.

Wei, T. and Simko, V.: R package “corrplot”: Visualization of a Correlation Matrix. (Version 0.90), 2021.

830

Wickham, H.: ggplot2: elegant graphics for data analysis., Springer, 2015.

Yom-Tov, Y. and Geffen, E.: Geographic variation in body size: the effects of ambient temperature and precipitation, Oecologia, 148, 213–218, https://doi.org/10.1007/s00442-006-0364-9, 2006.