Analysis of stratigraphical sequences at Cocina Cave (Spain) using rare earth elements geochemistry

This study investigates the stratigraphical sequence of Cocina Cave (Spain) employing and testing for the first time the capability of rare earth elements as markers of human activities in caves. Located in Dos Aguas (Valencian Community, Spain), Cocina Cave is characterized by the presence of several Holocene archaeological deposits from the final Mesolithic to the present day and is a pivotal site for understanding the socio‐ecological dynamics of the last hunter‐gatherer inhabitants of the Iberian Peninsula and the transition to pastoral and agricultural economies in the Western Mediterranean. However, the identification of strata from particular time‐periods in the cave is often difficult due to the homogeneity of layers, the poor archaeological record in some strata and the presence of severe disturbance phenomena. The methodological approach of this study consisted of cross‐referencing rare earth elements and other chemical markers with the archaeological stratigraphical interpretation, in an attempt to not only support the identification of the anthropic contribution to the formation of Cocina Cave strata, but also to characterize and confirm different natural and occupational episodes, particularly those associated with hunter‐gatherer, early agriculturalist, and shepherd activities. Sediments were collected from different excavation areas and analysed for major elements, trace elements, rare earth elements (REE), soil organic matter (SOM) amounts and pH. Multivariate statistics were employed to group samples according to their elemental profile, and these were then compared to the archaeological temporal interpretation. The obtained results showed that REE amount and fractionation geochemical processes were regulated by carbonates, phosphates and pH. The use of REE as markers was particularly useful as their concentrations and their calculated ratios and anomaly distributions were demonstrated to be highly consistent with the archaeological stratigraphical interpretation.

The sediment sources in caves result from physico-chemomechanical processes including dissolution, granular disintegration, external and internal inputs driven by wind, water transportation and mass movement. These are further modified by biological drivers including animals, plants and human actions, introducing and processing both organic and inorganic materials. The actions of these different factors result in a variety of sediments that are often prone to postdepositional movements and erosive processes that can truncate complete sequences (Goldberg & Sherwood 2006). In karst environments, the understanding of the cave stratigraphy is further complicated by the solubility of ash, carbonates and other minerals (Karkanas et al. 2000;Rasbury & Cole 2009;Canti & Huisman 2015;Ward et al. 2017;Berna et al. 2020).
Multi-element analysis of archaeological sediments has proven to be a valuable tool for understanding postdepositional processes and the interpretation of archaeological stratigraphies by supporting traditional fieldwork techniques (Lima da Costa & Kern 1999;Oonk et al. 2009;Pastor et al. 2016). Metal (i.e. calcium, strontium, copper and zinc) and nonmetal (i.e. carbon and phosphorous) elemental concentrations have often been employed as stratigraphical markers in open-air human settlements occupied during different periods (Parnell et al. 2002;Oonk et al. 2009;Fleisher & Sulas 2015;Dreibrodt et al. 2017). However, multi-elemental analyses in cave sediments to better understand stratigraphical processes and particularly to observe the human contribution to the development of the strata have been carried out in just a few cases. Finlayson et al. (2006) and Theden-Ringl & Gadd (2017) identified human occupation in Gorham's Cave (Gibraltar) and Wee Jasper (Australia), respectively, employing phosphorous, copper and zinc as anthropic markers. Homsey & Capo (2006) and Monge et al. (2015Monge et al. ( , 2016 employed phosphorous, strontium, zinc, copper, nickel and lead as markers to correlate sediment elemental concentration with the materials found in the excavations and identify anthropic activities. McAdams et al. (2020) developed an integrated geoarchaeological approach, including sediment multi-elemental analysis, to reconstruct the natural and human contribution in the sediment formation at Moong Cave (North Vietnam). Sheltered sites are subjected to different weathering conditions compared to open-air field sites (Goldberg & Macphail 2006) and consequently show different enrichment and leaching of elemental markers (Friesem et al. 2016).
Cocina Cave is located in Dos Aguas (Valencian Community, Spain; Fig. 1) and characterized by the presence of several Holocene archaeological deposits from the final Mesolithic to the present day. The presence of Mesolithic, Neolithic and Bronze Age preserved strata secures Cocina Cave as a pivotal site for understanding the socio-ecological dynamics of the last hunter-gatherer inhabitants of the Iberian Peninsula and the transition to pastoral and agricultural economies in the Western Mediterranean (García-Puchol et al. 2018a). However, the identification of anthropogenic strata in Cocina Cave, as is common in caves, is often difficult, due to the postdepositional processes that mask the original sediment properties, such as colour, texture and composition, as well as the depositional structures (Goldberg & Sherwood 2006). Furthermore, in the studied cave the homogeneity of the strata, the unequal archaeological record in some of these and the presence of severe bioturbation make the work of identifying anthropic layers very difficult for field archaeologists.
Rare earth elements (REE) have traditionally been employed with other trace elements to distinguish sediments associated with human occupation belonging to different contexts (Entwistle et al. 2000;Cook et al. 2006;Elmaleh et al. 2012;Nielsen & Kristiansen 2014;Butler & Dawson 2018;Dergacheva et al. 2018;Sulas et al. 2019). However, REE signatures have been recently employed by Gallello and colleagues (Gallello et al. 2013;2019a, b) in archaeological strata to determine differences between natural and human-induced In the Serpis Valley region (Spain), a marl-based geological area, REE soil analyses have been demonstrated to provide significant details about anthropogenic layers, suggesting that the archaeological interpretation of some dark brown deposits is evidence of a Neolithic region-wide agricultural system (Gallello et al. 2019b). REE results obtained in soils/sediments from the African subtropical region of Konso (Ethiopia) showed stratigraphical differences, already defined through field observations, and highlighted changes within the deposits developed by both natural and anthropogenic activities that were also affected by postdepositional processes (Gallello et al. 2019a). The combination of soil components affected by anthropogenic activities leaves 'fingerprints' in the REE chemistry and makes it possible to use these elements to distinguish between natural, disturbed and anthropogenic deposits. Therefore, REE are demonstrated to be excellent proxies for human activities in certain contexts, especially when compared to the more soluble/mobile calcium, sodium, strontium and potassium markers.
The aim of the current study was to develop a methodological approach based on multi-elemental analyses and focused especially on cross-referencing REE and other elemental marker concentrations with the archaeological stratigraphical interpretation in order to identify the anthropic contribution in the formation of Cocina Cave strata and to characterize different natural and occupational episodes. This analysis encompassed the rich archaeological sequences of Cocina Cave across hunter-gatherers (Mesolithic period), early agriculturalists (Neolithic period) and Bronze Age shepherd activities, whilst also taking natural and artificial sediment disturbance phenomena into account.

Geological and archaeological setting
Cocina Cave (from the Spanish name: cueva de 'La Cocina') is characterized by a wide cavity (30×22 m) carved into the La Ventana ravine, an irregular stream tributary to the Barranco Falón that joins the Xúquer river. The main course of the river at this point leaves the Iberic Range and travels 30 km at lower elevations to the Mediterranean Sea (Fig. 1).
The site is located in the massive calcareous Cretaceous ranges surrounding the Canal Valley, limited to the north by the Caballón Mountains; by the Dos Aguas Depression and the Caroig massif to the south, and the Quaternary deposits that form the alluvial plain of Valencia in the east. Granulometry and carbonate analyses were carried out by Fumanal (1979Fumanal ( , 1986 to study the sedimentary processes related to the cave, providing useful data to explain the formation process of the deposits. The obtained sedimentological data confirmed that the formation of the cave was related to karstic processes and linked to its drain function, while the posterior visible accumulated sediments from the Pleistocene to the Holocene mainly correspond to the fluvial action of the La Ventana ravine (Fumanal 1986). This can be observed in the trench area excavated by Fortea and also in different excavation profiles that reached Pleistocene levels, exposed in the northeast corner, close to the current entrance. Level X (10th) corresponds to the top of the Pleistocene deposits (without archaeological evidence) and it was radiocarbon dated in context C468 to 26 733AE2181 a BP and in C469 to 27 466AE2398 a BP. Holocene occupation levels are developed over a poorly interpretable depositional hiatus in the cave. Sporadic flows of water through the cave have been documented, as observed in 1978 and described by Fumanal (1986).
The remarkable archaeological discoveries of the 1940s (Pericot 1945) and the research subsequently carried out by Fortea (1973) have contributed to these archaeological deposits being considered one of the most important late Mesolithic sequences in the Western Mediterranean region. Since then, Cocina Cave has been a key site to understand the evolutionary history of the last hunter-gatherers in Mediterranean Iberia, and the appearance and development of farming and pastoral practices. Archaeological studies consisted of Pericot's excavation (1941Pericot's excavation ( -1945 (Pericot 1945), Fortea's investigation (1974Fortea's investigation ( -1981 (Fortea et al. 1987), and the recent campaigns (2015-2018) by García-Puchol, McClure and Juan Cabanilles (García-Puchol et al. 2015;2018b;Diez-Castillo et al. 2017;Pardo-Gordó et al. 2018;Cortell-Nicolau et al. 2020).
During the recent fieldwork, funded by the Museum of Prehistory of Valencia a systematic sampling strategy of vertical profiles selected from excavated trenches was undertaken. The chronological study carried out in the excavated deposits showed a wide development of the final Mesolithic occupations from the beginning of the 9th to the middle of the 8th millennium (García-Puchol et al. 2018a, b;Pardo-Gordó et al. 2018). Recently, radiocarbon dating has revealed the presence of human activities since the last centuries of the 8th millennium (Neolithic) and also during the 5th (Chalcolithic) and the 4th millennium (Bronze Age). However, a low level of preservation was observed in the Neolithic and Bronze Age layers, only identified in some residual parts of the cave, probably due to natural and modern anthropic induced taphonomic processes (Pardo-Gordó et al. 2018).

Material and methods
Cocina Cave soil samples were collected from crosssections of different excavated areas, including cleaned Fortea's trench profiles and newly excavated pits, which were subsequently analysed to determine major, minor and trace element (including REE) contents. Additionally, pH and soil organic matter (SOM) content of each sample were measured. Conventionally employed soil elemental anthropic markers (Ca, P, Ba, Cd, Zn, Cu, Ni, Mn and pH) results have been cross-referenced with REE data and the chemical processes influencing the REE behaviour were assessed. The data were processed through multivariate statistics techniques such as cluster analysis (CA) and factor analysis (FA) to understand the differences and similarities among the samples.

Sampling
Seventy-one samples (Table 1) were collected from the cross-sections of five (FT, S3A, S3B, S3C and S4) excavation areas (Fig. 2). The studied sections (all of them between 0.7 and 1 m deep), were scraped clean and sampled. Then, a representative amount of sediment (approximately 20 g) was removed with a lab spoon from each layer following the identified archaeological stratigraphy, working from the bottom to the top. Each sample was given a unique number, with the lowest layer having the lowest number. The colour of each dry soil sample was classified using the Munsell Color System. Each sample taken was recorded indicating the possible origin (natural and disturbed) and strata-related period (Modern, Bronze, Neolithic and Mesolithic) based on the archaeologists' fieldwork interpretation. These attributions were based on the excavated material's cultural and biological charred remains, radiocarbon dating, and stratigraphical and archaeological information derived from fieldwork carried out across the site. When the attributions to a period or possible origin of a sample is uncertain a '?' was added (Table 1).
Twenty-three samples were collected from a trench originally excavated during Fortea's fieldwork (FT) that was reopened and extended in the 2015 excavation season. The samples were collected from three crosssections, characterized by brown and orange strata, along four soil columns (FTa1-9, FTb1-4, FTc1-5, FTd1-5) and corresponding to the archaeological units CO-0, CO-1, CO-4 and CO-5. According to the archaeological interpretation of the stratigraphy, most of the samples were collected from natural sediments, from modern stable strata or from disturbed earth. Samples FTc1-2 from unit C0-1 and FTb1-2 from unit CO-0 were interpreted as potentially belonging to the Mesolithic period after the obtained radiocarbon dates (8200 to 7700 cal. a BP; García-Puchol et al. 2018a).
The area S3 is stratigraphically more complex. It was divided into three different subareas (S3A, S3B, S3C) excavated during 2015 fieldwork next to another of Fortea's trenches (sondeo 3, Fig. 2) close to FT. This area (S3) was covered with more than 1 m of modern upper layers excavated during previous fieldwork campaigns. The identified archaeological layers, in the three well-preserved sections (S3A, S3B, S3C), were 15-20 cm thick. The archaeological record has been attributed to span the Mesolithic, the Neolithic, and to the Bronze Age and radiocarbon dated to the early Neolithic (7300 cal. a BP) (García-Puchol et al. 2018a). Six samples were collected from S3A (1-6). The stratigraphy is characterized by the presence of natural orange sediments from archaeological unit CO-0 at the bottom (S3A1), and by subsequent orange or brown strata from unit CO-1 (S3A2) interpreted as Mesolithic, unit CO-2 as early Neolithic (S3A4) and CO-3 (S3A5 and S3A6) interpreted as late Neolithic or Bronze Age. There is uncertainty about the archaeological interpretation of the sample S3A3 taken in CO-2, but it may belong to an early Neolithic deposit. Another six samples were collected from orange and dull brown strata from S3B, stratigraphically related to units CO-0 (S3B1), interpreted as natural, CO-1 (S3B2) Mesolithic, CO-2 (S3B4-5) early Neolithic, and CO-3 (S3B6) late Neolithic or Bronze Age; S3B3 from CO-1 was interpreted with some uncertainty as belonging to the Mesolithic period. S3C is also characterized by orange to brown strata, here six samples were taken from units CO-0 (S3C1) interpreted as natural, CO-1 (S3C2) as Mesolithic, CO-2 (S3C4-5) as early Neolithic and CO-3 (S3C6) as Bronze Age. Sample S3C3 from CO-2 was archaeologically interpreted as belonging to the early Neolithic.
Thirty samples were collected from area S4, excavated during the 2018 fieldwork, close to the current entrance of the cave and sampled from four columns (S4a1-8, S4b1-8, S4c1-8, S4d1-6). While the samples from S4a and S4b, employed as control samples, were taken from the archaeological unit CO-5 and are probably made up of disturbed earth, related to the built modern wall, most of the samples from the other two columns (S4c1-8, S4d1-6) belonging to unit CO-1 are attributed to Mesolithic phases of occupation. Two samples (S4c1-2) from the bottom of the column S4c in unit CO-1 were interpreted as natural sediments. S4c8 and S4d6 from unit CO-1 are of uncertain attribution and it is not clear from the archaeological record if they are from Mesolithic or disturbed deposits. Most of the samples presented a brownish colour except S4c2-3, which were orange. It should be noted that samples from 'disturbed earth' may be a composite of modern and ancient material and therefore have a mixed geochemical signature from the various components.
Portable X-ray fluorescence analysis (pXRF) All analyses were carried out on sediment samples that had been pulverized and homogenized with an agate mortar and pestle. The elemental concentrations of the BOREAS  Table 1. List of the samples in stratigraphical order for each sampled section, area of provenance, archaeological units (Units), archaeological interpretation (A.I.), sediment description (S.D.), and Munsell colour code (Colour). Colour (according to Munsell Color System): 7.5YR5/8 and 5/6 = bright brown; 7.5YR6/8 and 6/6 = orange; 7.5YR6/4, 7/3 = dull orange; 7.5YR6/3, 5/4 and 5/3 = dull brown; 7.5YR6/2 = greyish brown; 7.5YR7/2 = light brownish grey. Archaeological interpretation: NA = natural; MO = modern stable (livestock soils); ST = disturbed earth; ME = Mesolithic; NE = Neolithic; BR = Bronze Age; ? = uncertain archaeological interpretation. powdered samples were obtained by using a S1 Titan portable energy dispersive pXRF from Bruker (Kennewick, Washington DC, USA) equipped with an Rh X-ray tube and X-Flash ® SDD detector. The Geochem-trace calibration, as supplied by the manufacturer, was used to perform the quantitative analyses and determine Al, Si, P, K, Ca, Ti, Fe and Zr concentrations. The accuracy and precision of the analysis were tested using NIM-GBW07408 Soil certified reference material (Table S1). All the obtained elemental results have less than 7% error compared to the certified results (Table S1), except that of P, which was less than the detection limit.
Inductively coupled plasma mass spectrometry (ICP-MS) Inductively coupled optical emission spectrometry (ICP-OES) and inductively coupled plasma mass spectrometry (ICP-MS) are the analytical techniques commonly employed to measure chemical elements with the advantage of lower instrumental limits of detection compared to the other techniques such as Xray fluorescence spectroscopy (XRF). Several sediment preparation methods such as weak acid extraction (Oonk et al. 2009;Vyncke et al. 2011) and both partial (Dirix et al. 2013(Dirix et al. , 2016Linderholm et al. 2019) and  (Sulas et al. 2019), while others avoid more aggressive acid digestion that may cause an overwhelming geogenic signature masking the anthropic signature (Vyncke et al. 2011;Gallello et al. 2019a).
In this work a preliminary study was undertaken to identify the preferred digestion method, to take into account the geological environment of the cave. The homogenized and powdered sediment samples were digested both using an aqua regia acid partial attack and a HF based total attack for trace elemental analysis by ICP-MS. This test comparing the results of the two digestions showed that no significant improvements were obtained through the employment of the total digestion method, which was longer and more hazardous due to multiple steps and to the use of hydrofluoric acid (Table S2A). On the contrary, elemental differences, based on the total digestion results (Table S2B), within the anthropogenic stratawere diffuse and less easy to interpret. As suggested by others (e.g. Middleton 2004), we observed that a total digestion provides an elemental profile that is dominated by the geogenic mineral fraction that overwhelmed the anthropogenic contribution, particularly due to the presence of clay minerals found in the Cocina Cave sedimentological environment (Fortea et al. 1987).
Finally, the digestion method employed to obtain the results presented in this study consisted of the addition of 1.35 mL HCl (37%, Scharlab) and 0.45 mL HNO 3 (69%, Scharlab) to 0.15 g sample in a glass tube; subsequently heated in a boiling water bath for approximately 40 min.
Once cooled to room temperature, the solutions were poured into plastic tubes and made up to a volume of 25 mL with purified water. The concentrations of Ba, Bi, Cd, Cr, Co, Cu, Pb, Li, Mn, Mo, Ni, Sr, Tl, V, Zn, REE (La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb and Lu), Sc and Y were determined in the solution after filtration. Two multi-element stock solutions (5% HNO 3 , Sharlab) for the ICP analysis containing 100 mg L −1 of the analysed elements were used after dilution as standards for calibration. The concentration of trace elements in individual standards ranged between 1 and 600 μg L −1 except for REE, Y and Sc, which ranged between 1 and 100 μg L −1 . As an internal standard, to correct for matrix suppression, 10 μg L −1 of Rh was added to each sample and calibration point; internal standard normalization was carried out automatically by instrument software.
The analyses were performed using an Elan DRCII inductively coupled plasma mass spectrometer by Perkin Elmer (Concord, Ontario, Canada). NIM-GBW07408 Soil certified reference material was measured to test the accuracy and precision of the analyses. Almost all the analysed elements (aqua regia) are statistically in the range of the certified results (Table S2A). Although the total digestion method (HF) also provided results very close to the certified ones, they were not employed in this study.

Soil organic matter (SOM) and pH analyses
Soil organic matter (SOM) was estimated by oxidation through potassium dichromate as suggested by Radojevic & Bashkin (2006). Sediment pH was measured in a 1:2 (soil/distilled water) extract, after shaking for 5 min using a MicropH2000 pH-meter by Crison.

Statistical data analysis
Raw data are reported in Table S3 (major elements, SOM and pH), Table S4 (trace elements) and Table S5 (REE data and REE ratios).
The first exploratory analysis consisted of identifying the relationships between total REE (TREE) and the different major elements (Figs 3, S1, S2). Pearson correlation indexes (r) and Wilcoxon rank sum test with a standard α = 0.05 as the level of rejection were used (Table S6).
In order to evaluate the presence of different groups related to the phases of occupation of the cave, a hierarchical cluster analysis (CA) was run. CA was employed separately on the FT, S3A, S3B, S3C and S4 areas. Data were autoscaled prior to analysis, Euclidean distances between each pair of samples were calculated and complete linkage was employed for the clustering method. CAs (Figs 4-6) were carried out employing as variables REE (Figs 4A, 5A, 5C, 5E, 6A) and Ca, P, Ba, Cd, Zn, Cu, Ni, Mn and pH (Figs 4B, 5B, 5D, 5F, 6B) due to their previously reported use as anthropic markers (Table 2). Data analysis was performed in R (version 3.6.2; R Core Team 2019) and dendrograms were made by the R package 'factoextra' (version 1.0.7; Kassambara & Mundt 2020).
To reinforce the reliability of data interpretation, factor analysis (FA) was also employed to observe differences between the CO-0, CO-1, CO-2, CO-3, CO-4, CO-5 strata without separation by areas (Fig. 7A, B). After computing the eigenvalues for each factor, n = 3 factors were used, which accounted for 82.46% of the variance. The data were rotated using the Varimax method to aid interpretation of the role of variables. In order to deal with abnormalities within the data, each variable was winsorized, trimming its edges at the 0.05 and 0.95 quantiles, thus adjusting possible outliers. One observation (S4b8) remained far off the group mean according to Mahalanobis distance even after winsorization andwas removed (for more details see Data S1). After this, data were scaled, as is a standard procedure for carrying out FA with variables of different magnitudes. A correlation matrix was used, where seven variables (K, Fe, TREE, SOM, pH, Ce/Ce* and Sm/Yb) were removed due either to their low contribution to the model or their high degree of multicollinearity. By doing this, an observations-variables (70-13) ratio of 5 was obtained, within the recommended range (Habing 2003). Sample adequacy was further tested by calculating the overall Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy (MSA) statistic, which was found to be 0.71, supporting the suitability of the data for FA. Thus, the final FA suite of variables were: Al, P, Ca, Pb, Ba, Cd, Zn, Cu, Ni, Mn, Eu/Eu*, La/Yb, La/Sm. The main R packages used for FA were 'psych' (version 2.0.12; Revelle 2020) and 'PerformanceAnalytics' (2.0.4; Peterson & Carl 2020).

Rare earth elements (REE) parameters
In order to observe REE fractionation, lanthanide concentrations were normalized by using post Archaean Australian shale (PAAS) REE values, following the values reported by McLennan (1989) (e.g. La n : normalized lanthanum).
Ce and Eu anomalies, resulting from differential reactions associated with reduction-oxidation (Ce 3+ or 4+ and Eu 2+ or 3+ ), were calculated as suggested by Hinz & Kohn (2010): (1) Anomalies are considered positive or negative depending on whether the values are above or below 1, respectively.

Chemical analysis
Higher concentrations of rare earth elements (REE , Table  S5) were detected in FT (total REE concentration, TREE: 128AE26 mg kg −1 ), S3C (TREE: 117AE17 mg kg −1 ), S3B (TREE: 93AE13 mg kg −1 ) than in S3A (TREE: 72AE6 mg kg −1 ) and S4 (TREE: 68AE13 mg kg −1 ). Similar results were observed in LREE, MREE and HREE, all of them These results broadly coincide with expectations from the relationship between TREE and different major elements, SOM and pH (Fig. 3 for S4; Figs S1, S2 for FT, S3A, S3B and S3C). Significant positive correlations between TREE and Al, Si, K, Ti and Fe are consistent throughout the data analysis and confirmed by Wilcoxon analysis (Table S6). However, TREE is negatively correlated with Ca. Although similar Pearson's coefficients can be observed for most of the elements in S3A, S3B and S3C (Fig. S2), due to the small sample size in these sectors (n = 6) at the boundary of acceptability for nonparametric tests and usually insufficient for linear correlations, we focus the analysis on FT (Fig. S1) and S4 (Fig. 3) samples. Even in these cases, some observations should be made. For example, in the correlation TREE~P on the FT sector we must take into account the high number of <LOD values in the original data. We have substituted these values with min(P)/2 values, which then weights the data close to zero. Furthermore, the outlier value of FTb3 in SOM clearly influences the linear correlation, but in this case we have preferred to maintain the outliers as part of the data rather than reduce the sample size or moderate them. Most of these problems are not apparent in S4 (Fig. 3).
The obtained results show elemental differences among the sampled areas, probably related to the micro-environmental characteristics in the different areas of the cave.

Cluster analysis (CA) of the studied areas
To observe the capability of REE to distinguish strata and evaluate differences, CA was carried out analysing each area separately (Figs 4-6).
The results of CA for FT using only the REE or the anthropic markers as variables are shown in Figs 4A-B. Four main clusters can be observed in the REE dendrogram (Fig. 4A), but the CA could not separate units CO-0, CO-1, CO-4 and CO-5. However, while in three clusters the different sediments are grouped independently of their origin, in the cluster at the right of the dendrogram units CO-0 (natural) and CO-1 (Mesolithic) form one subgroup and CO-4 (modern) and CO-5 (disturbed) form another. If we examine the dendrogram created by using anthropic markers (Ca, P, Ba, Cd, Zn, Cu, Ni, Mn and pH) ( Table 2) as variables (Fig. 4B) no clear grouping between the different units can be found. The difficulty in separating the different archaeological units can be related to postdepositional phenomena. This area of the cave was affected by later activities, such as a modern animal pen, excavating sediment for manuring fields outside of the cave and earth moving, that may have increased elemental migration between the horizons, making the interpretation of the stratigraphy more difficult.
Concerning the area S3, according to the field observations, S3A is defined by four units (CO-0, CO-1, CO-2 and CO-3). The REE CA results follow the stratigraphical sequence (Fig. 5A). S3A6 (CO-3), recognized as a Bronze Age layer, has the greatest distance from the rest of the samples (S3A1-5) due to its very low REE levels and S3A1 (CO-0) groups together with sample S3A2 (CO-1). Another subgroup is composed by S3A3 (CO-2) separated from S3A4 (CO-2) and S3A5 (CO-3). Samples S3A3-5 have lower concentrations in almost all REE compared to natural samples S3A1 and S3A2 (Fig. 5A). Differences between the S3A3 and S3A4 samples are due to a slightly higher concentration for almost all REE in the latter sample. S3A3-4 cluster together in the dendrogram of the anthropic markers (Fig. 5B), separated from S3A5, as they are enriched in Ca, P, Ba, Cd, Zn the abovementioned samples (Tables S3, S4). The sample S3A6 also clusters separately from the rest of the samples and shows the highest concentrations for all the selected elements (except Ca) as well as the highest pH (Table S3).
The area S3B is also defined by the same sequence as S3A and REE CA results follow the stratigraphy (Fig. 5C). S3B1 (CO-0) has the greatest distance from the rest of the samples (S3B2-6) due to its very high levels of REE. S3B2-3 (CO-1) and S3B4-5 (CO-2) are clearly clustered separately and are distinct from the S3B6 sample (CO-3), which is identified as aBronze Age deposit. S3B2-3 samples from CO-1 are characterized by higher REE levels than those of S3B4-5 from unit CO-2. The dendrogram of the anthropic markers (Fig. 5D) clusters S3B6 separately from the other samples, and this sample is more enriched in P, Ba, Cd and Zn (Tables S3,  S4). However, the dendrogram did not show any differences amongst units CO-0, CO-1 and CO-2.
The area S3C is defined by CO-0, CO-1, CO-2 and CO-3. The REE CA results show a certain coherence with the identified layers (Fig. 5E). The sample S3C1 (CO-0) clusters separately but is somewhat close to S3C2 (CO-1) and it presents the highest levels of REE. Another group is composed by samples S3C3-5 (CO-2) and S3C6 (CO-3), and is characterized by lower concentrations in REE compared to the CO-0 and CO-1 units. S3C4-5 are also separated from S3C3 and S3C6. The dendrogram of the anthropic markers (Fig. 5F) is slightly more consistent with the stratigraphy than REE cluster. S3C4-5 (CO-2) are classified close to S3C3 (CO-2) and separated from S3C6 (CO-3), the latter with higher concentrations of P, Ba, Cd, Zn and Cu compared to the other samples (Tables S3, S4).
The archaeologists suspected that these areas (S3A, S3B and S3C) of the cave were employed as a pen during the Bronze Age and this would give rise to inadvertent manuring by animals. The higher concentrations of anthropic markers in CO-3 layers compared to CO-1 and CO-2 units suggest a consistent change in the impact of human activity carried out in this part of the cave, which BOREAS is compatible with the above-mentioned hypothesis. (Wilson et al. 2007;Oonk et al. 2009;Rondelli et al. 2014;Lancelotti et al. 2017).
The results of the REE CA for area S4 (divided into soil columns a, b, c and d) are shown in the dendrogram of Fig. 6A. Units CO-0, CO-1 and CO-5 were identified in this area. In particular, the sampled columns S4a and S4b from CO-5 (disturbed earth) related to the construction of the enclosing fence in the cave entrance. In order to avoid misleading results, the samples S4a6-8 (CO-5) were excluded from the CA due to the anomalous concentrations of the anthropic markers, probably related to modern contamination. The REE CA shows that S4c1-2 are grouped together close to S4b8 and S4c3. Except for S4c3, all the samples from the unit CO-1 (S4c4-7, S4d1-6) were identified as belonging to Mesolithic deposits, clustering also together with samples S4c8 (CO-1). Furthermore, the samples S4a1-5, S4b1 and S4b5-7 (CO-5) are separated from the CO-1 unit, clustering together and separated from S4b2-4 (also from CO-5). Samples S4c1-2 (CO-0) show the highest REE concentrations, while sediments from CO-5 have the lowest levels and samples from CO-1 are mid-way between the other two units (Table S3).
In the anthropic markers CA (Fig. 6B), S4c1-2 (CO-0) group with S4b8. S4c and S4d samples (CO-1) cluster closely, including sample S4c8 (interpreted as Mesolithic with uncertainty). In general, with some exceptions, samples from unit CO-1 showed higher concentrations of Ca, P, Ba, Cd, Zn and Cu compared to unit CO-0, while lower concentrations were measured for most markers compared to samples from CO-5 (Table S4).

REE parameters and anthropic markers defining differences between strata
Plotting factor analysis scores for individual samples proved informative (Fig. 7A, B). In Fig. 7A, it may be observed that FA-1 clearly separates CO-0 (natural) and CO-4 (modern) samples from archaeological samples CO-1 (Mesolithic); CO-2 (Neolithic); CO-3 (Bronze Age) and . The mixed array of CO-0 and CO-4 suggested modern samples were likely natural derived material. FA-2 scores clearly separate CO-2 (Neolithic) and CO-3 (Bronze Age) from CO-1 (Mesolithic), which tends to group with CO-5 (disturbed). CO-5 presents the highest score values in FA-2.
In Fig. 7B FA-3 indicates a degree of separation of CO-3 (Bronze Age) from CO-1 (Mesolithic) and CO-2 (Neolithic); the highest FA-3 values are also associated with CO-5 (disturbed) material. The combination of FA-1 and FA-3 scores confirms the relationship between CO-0 (natural) and CO-4 (modern) material, while also providing a neat separation between CO-0 and CO-4 on the one hand, and archaeological samples on the other. Specific samples of note included: FTc1 and FTc2, which were described as either CO-0 (natural) or CO-1 (Mesolithic), can be confirmed as CO-0 (natural), while samples S3A1 and S3C1 are more likely to be CO-1 (Mesolithic) than CO-0 (natural). Sample S3A3 was classified as CO-2 (early Neolithic) with some reservations, whereas its values seem to point to a CO-1 (Mesolithic) attribution; contrary to sample S3B3, initially classified as CO-1 (Mesolithic) but which seems to group better with CO-2 (Neolithic) samples. The initial CO-1 (Mesolithic) attribution for samples S4c8 and S4d6 is reinforced by the FA results.
Inspection of the FAvariable loadings in Table 3 gives insight into the controlling geochemistry. FA-Ld-1 was dominated by Al and Ca in an inverse relationship with P and Cu also significant. This suggested clay (Al) and carbonates (Ca) and possibly bone (Ca, P) proportions were a mineralogical factor. FA-Ld-2 was dominated by Ni, La(LREE)/Sm(MREE); Eu/Eu*; with La(LREE)/ Yb(HREE) also significant, demonstrating the role of REE geochemistry in separating the different archaeological strata. FA-Ld-3 was dominated by Cd, Ba and Mn, all classic anthropic elements.

Geochemical synthesis
We observed that in the FT area the unit CO-4 is affected by recent activities such as a modern pen and unit CO-5 by sediment removal that has almost  certainly increased the elemental migration through the section and made the chemical interpretation of the stratigraphy more difficult. Nevertheless, except for FT, in each area the samples collected from the anthopic units CO-1, CO-2 and CO-3 are depleted in REE compared to the samples from corresponding natural levels (CO-0) (Fig. 8). This may be because REE have been mobilized then sorbed onto the clay mineral phase (Gallello et al. 2019a), resulting in higher levels of REE in natural sediments. We can observe that CO-0 is often characterized by slightly higher Ce n /Ce* ratios, which are negatively correlated to P (correlation coefficient 'r' for P vs. Ce n / Ce* is −0.90 for S3A, −0.93 for S3B and −0.83 for S3C, p < 0.05). Positive Ce anomalies are typical of clay with high Fe content ('r' for Fe vs. Al is 0.93 in S3A, 0.94 in S3B, 0.90 in S3C and 0.95 in the selected samples of S4, p < 0.05) (Pattan et al. 2005). In S3A, S3B and S3C, natural layers are different from the anthropic layers especially because natural strata have higher LREE levels and higher Ce positive anomaly values.
The slight depletion in Ce compared to its neighbouring elements in anthropic units CO-1, CO-2, CO-3 could be explained by the presence of authigenic phosphate minerals linked to the diagenesis of human activity products such as bone or ash (Schiegl et al. 1996;Karkanas 2010). Authigenic phosphate minerals are known to retain moderate Ce negative anomalies in sediments (Pattan et al. 2005).
S3A, S3B, S3C samples are also characterized by slightly negative to neutral Eu anomalies and most samples from anthropic strata have a slightly higher enrichment in Eu than the natural ones (Fig. 8).
Although all S4 samples have slightly positive Eu anomalies, S4c1 has the lowest Eu n /Eu*. the Eu anomalies are likely driven by phosphate minerals, which can complex Eu during diagenesis (Hu et al. 2006). Weathered products and secondary phosphate minerals can contribute to the development of the strong REE migration during the weathering process of carbonate rocks (Pattan et al. 2005;Jiyan & Ruidong 2010). LREE and MREE are enriched over HREE in S3A, S3B, S3C and S4 (Fig. 8). The depletion in HREE is higher in CO-0 (natural) for S3A, S3B and S3C areas where there is a clear negative correlation of La n /Yb n and Sm n /Yb n with P ('r' is −0.51 and −0.42 for S3A, −0.88 and −0.90 for S3B, and −0.70 and −0.80 for S3C). Furthermore, Ca concentrations are lower in CO-0 than the anthropic units CO-1, CO-2 and CO-3 and this may explain the lowest values of La n /Yb n and Sm n /Yb n ratios in some anthropic units as it may be controlled by the precipitation of carbonate minerals enriched in HREE (Compton et al. 2003;Zhaozhou et al. 2012). Carbonate and phosphate minerals can also form complexes with REE under alkaline conditions, especially with LREE enrichment over HREE. Since the carbonates and phosphates in Cocina anthropic layers are probably precipitates resulting from prolonged human activities, the REE ratios are indicative of an anthropic impact.
In the S3A, S3B and S3C areas, we can also observe that the main differences between CO-1 (identified as Mesolithic) and CO-2 (early Neolithic) strata are related to the lower La n /Yb n fractionation ratio of CO-2 (Table 3). By contrast, in unit CO-3 (identified as a Bronze Age deposit), the lowest LREE levels were measured, which is reflected in a lower La n /Yb n ratio compared to CO-1 and CO-2.  CO-1 ME/ST? ME ME ME S4c7

Synthesis of the archaeological information
CO-1 ME ME ME ME S4c6 CO-1 ME ME ME ME S4c5 CO-1 ME ME ME ME S4c4 CO-1 ME ME ME ME S4c3 CO-1 ME NA/ST? ME ME? S4c2 CO-1 ME/ST? ME ME ME S4d5 CO-1 ME ME ME ME S4d4 CO-1 ME ME ME ME S4d3 CO-1 ME ME ME ME S4d2 CO-1 ME ME ME ME S4d1 CO-1 ME ME ME ME made (A.I.+REE+Markers) taking into account the elemental results and the archaeological interpretations; inconsistences in the interpretation are also highlighted (black border rectangle). Differences between CO-1 (Mesolithic) and CO-2 (early Neolithic) confirmed by REE (Table 4) have allowed the development of hypotheses regarding the human occupation impact in Cocina Cave by crossreferencing these data with the published archaeological literature. As suggested by the remains found in the cave (Fortea et al. 1987) in early Neolithic deposits, hunting was still part of the daily activity of the inhabitants. Therefore, it is possible that the activities in Cocina Cave were very similar during the Mesolithic and Neolithic and that the identified small differences are just due to the intensity of the occupation in this area of the cave rather than a reflection of different activities as evidenced by the minor fractionation of LREE over HREE observed in the S3B and S3C early Neolithic strata. Furthermore, the low REE levels, low fractionation ratio of LREE and MREE over HREE together with low Ca and high levels of P and also Zn and Cu (Oonk et al. 2009) in unit CO-3 (interpreted as late Neolithic-Bronze strata) confirm anthropic activity changes in the cave. These data, crossreferenced with the archaeological records, support the hypothesis that Cocina Cave was employed as a stable for animals during the late Neolithic and Bronze Age. Sectors S4c and S4d follow a similar pattern to S3A, S3B, and S3C; CO-1 (Mesolithic layers) are higher in P and Ca and lower in REE concentrations, and there is a slightly negative Ce anomaly compared with unit CO-0 (natural), perhaps confirming activities based on hunting, bone manufacturing and meat processing ( Fig. 8 and Table 4), already archaeologically recorded by Fortea et al. (1987).

Conclusions
For the first time, rare earth elements (REE) analysis of sediments has been successfully carried out in a prehistoric cave, and specifically at the key site of Cocina Cave, where important traces of the last hunter-gatherers in Mediterranean Iberia have been found. The obtained data provide pivotal results about the impact of anthropic activities at the site and showed that the analysed cave areas present different sediment elemental compositions, likely caused by several factors such as environmental conditions inside the cave, human activities, and recent soil movement. However, a deeper understanding of the chemical mechanisms generating differences within strata is elucidated by cross-referencing REE with Ca, P, Ba, Cd, Zn, Cu, Ni and Mn data. REE amount and fractionation geochemical processes were regulated by carbonates, phosphates and pH. REE as geochemical markers proved to be particularly useful, with results that in some sections are more consistent with the archaeological stratigraphical interpretation than typical traditional anthropic elemental markers (Table 4). In both units CO-1 and CO-2, the intensity of Mesolithic and early Neolithic hunting activities in Cocina Cave is underlined by fractionation of LREE over HREE and seems to be reduced during the Neolithic period. This suggests that during the early Neolithic, hunting, although less intense, was still part of the daily activity of the inhabitants as confirmed by the remains found in the cave related to the manufacturing of bone artefacts or meat processing activities. Furthermore, REE concentrations and fractionation ratios, together with traditional markers Ca, P and also Zn and Cu, in unit CO-3 confirm the change of use of the space in the cave during the Bronze Age. These data when crossreferenced with the archaeological record confirm that Cocina Cave was used as a stable for animals.
This study shows the important role of REE in archaeological sediments as a complementary tool to better understand the human contribution to soil/sediment strata formation. Although the use of REE geochemistry in archaeology is growing, it is still limited but may prove a valuable tool to overcome the limits of traditional archaeological or geochemical methods to separate anthropic from natural processes, and to distinguish deposits formed by different human activities. The authors are grateful to the reviewers and the editor for their suggestions, comments and feedback that have greatly enhanced this manuscript. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Author contributions. -GG developed the project idea, the sampling strategy and wrote the paper. MR sampled the sediments, and, with MLC and AP, carried out the analysis and processed the analytical data. OGP, SBMcC and ADC wrote the archaeological setting and the archaeological interpretations. SC participated in the geochemical interpretation. ACN carried out data analysis and statistics.

Supporting Information
Additional Supporting Information may be found in the online version of this article at http://www.boreas.dk.   Table S1. The accuracy of pXRF analysis was evaluated from the analysisi of the soil CRM NIM GBW07408. The obtained and the certified concentrations of the analysed elements are shown in the table. Values expressed as percentages (%). Values expressed as mg kg −1 are marker with a star (*). <LOD: below the detection limit. Note: P levels measured in the standard soil are below the limit of detection (<0.01*).  Table S3. Major elements by pXRF and soil organic matter (SOM) concentrations, expressed as wt%, and pH. <LOD = less than the limit of detection. Table S4. Trace elements concentrations, after Aqua Regia digestion by ICP-MS, expressed as mg kg -1 . Zr concentrations were measured by pXRF. Table S5. REE concentrations, after Aqua Regia digestion by ICP-MS, expressed as mg kg -1 and derived parameters (dimensionless as ratios). TREE is the total (sum) of all REE. Table S6. Pearson correlation indexes (r) and Wilcoxon rank sum test with a standard α = 0.05 as the level of rejection.