# Towards a formation model of the Neanderthal symbolic accumulation of herbivore crania: Spatial patterns shaped by rockfall dynamics in Level 3 of Des-Cubierta Cave (Lozoya valley, Madrid, Spain)

Lucía Villaescusa, Enrique Baquedano, David M Martín-Perea, Belén Márquez, M Ángeles Galindo-Pellicena, Lucía Cobo-Sánchez, Ana Isabel Ortega, Rosa Huguet, César Laplana, M Cruz Ortega, et al.

► To cite this version:

Lucía Villaescusa, Enrique Baquedano, David M Martín-Perea, Belén Márquez, M Ángeles Galindo-Pellicena, et al. Towards a formation model of the Neanderthal symbolic accumulation of herbivore crania: Spatial patterns shaped by rockfall dynamics in Level 3 of Des-Cubierta Cave (Lozoya valley, Madrid, Spain). Archaeological and Anthropological Sciences, 2026, 18, pp.16. 10.1007/s12520-025-02382-5. hal-05458037

HAL Id: hal-05458037

https://hal.science/hal-05458037v1

Submitted on 14 Jan 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Towards a formation model of the Neanderthal symbolic accumulation of herbivore crania: Spatial patterns shaped by rockfall dynamics in Level 3 of Des-Cubierta Cave (Lozoya valley, Madrid, Spain)

Lucía Villaescusa¹,² · Enrique Baquedano²,³ · David M. Martín-Perea⁴ · Belén Márquez² ·

M. Ángeles Galindo-Pellicena¹,² · Lucía Cobo-Sánchez⁵ · Ana Isabel Ortega⁶ · Rosa Huguet⁷,⁸,⁹ · César Laplana¹,² ·

M. Cruz Ortega¹⁰ · Sandra Gómez-Soler²,¹¹,¹² · Abel Moclán³,¹³ · Nuria García¹⁰,¹⁴ · Diego J. Álvarez-Lao¹⁵ ·

Rebeca García-González¹⁶ · Laura Rodríguez¹⁶,¹⁷ · Alfredo Pérez-González³ · Juan Luis Arsuaga¹⁰,¹⁴

Received: 18 July 2025 / Accepted: 27 November 2025
© The Author(s) 2026

## Abstract

Understanding formation processes is crucial for interpreting sites with complex sedimentary histories and exceptional archaeological records, such as Level 3 of Des-Cubierta Cave. This Middle Paleolithic unit contains an assemblage of anthropically modified ungulate horned crania, Mousterian lithics, and evidence of fire use, all preserved in a clast-supported gravel deposit shaped by successive rockfalls. This study integrates geostatistical analyses with traditional spatial and taphonomic methods to examine the cone-shaped sedimentary structure that dominates the level and its influence on the spatial distribution and preservation of archaeological materials. The results reveal distinct spatial patterns for geological and archaeological materials, indicating separate formation dynamics. Size-based spatial analyses of boulders characterize the morphology of the conical structure, highlighting size sorting and intensity variations that may indicate sedimentary hiatuses. Bone refitting analysis suggests limited post-depositional movement, with material distribution shaped by the conical sedimentary structure and karst gallery morphology. Variability in crania preservation correlates with proximity to the cone's apex and elevation, with better-preserved specimens in central areas and more fragmented remains in zones affected by erosive and edaphic processes. These findings underscore the value of integrating geostatistical and traditional archaeological approaches to advance interpretations of spatial and temporal patterns in karst environments and provide a methodological approach for examining sites with similar sedimentary histories. Moreover, the results support the interpretation of a recurrent, culturally motivated Neanderthal behaviour centred on the deliberate accumulation of large ungulate crania, pointing to a symbolic dimension in these practices.

**Keywords** Formation processes · Depositional processes · Spatial analysis · Bone refit · Neanderthals

## Introduction

Understanding formation processes is essential for any archaeological site (Stein 2001), but it becomes particularly important in sites with complex sedimentary histories and unique cultural records, such as Unit 3 of Des-Cubierta Cave (Lozoya Valley, Madrid, Spain). The archaeological record here is dominated by at least 35 anthropically modified ungulate crania, accompanied by a Mousterian assemblage, and evidence of fire use (Baquedano et al. 2023). These crania,

along with associated materials, were found embedded in a 2-meter-thick lithostratigraphic unit, where the main sedimentary process involved successive rockfalls, resulting in a clast-supported gravel deposit that appears poorly sorted to the naked eye, with cobbles and boulders and minimal interstitial matrix (Baquedano et al. 2023; Martín-Perea et al. 2025). According to Bertran et al. (2019), this sedimentary context falls under the accumulation-transit type, characterized by reworking processes in which the final archaeological distribution reflects both the initial anthropic pattern and
subsequent sedimentary dynamics. Thus, rockfall processes not only deposited sediments but also likely influenced the redistribution of archaeological remains. In this context, where the archaeological record suggests non-subsistence behaviors linked to the symbolic world of Neanderthals, preserved within a multi-episodic depositional context indicating the persistence of this behavior over time, and where the spatial organization of remains may have been affected by sedimentary processes, a thorough study of formation processes is therefore crucial.

The study of formation processes is typically addressed through various approaches and methodologies, beginning with geoarchaeology to reconstruct the evolutionary histories of sites (Woodward and Goldberg 2001; Goldberg and Macphail 2005; Rapp and Hill 2006). Archaeological taphonomy, particularly in its assessment of site disturbances (Zilhão et al. 2015; Bertran et al. 2019), is applied to both bone (Isaac 1983; Lee Lyman 1985; Binford et al. 1988; Villa 2004; Domínguez-Rodrigo et al. 2007; Lloveras et al. 2011; Arriaza et al. 2017; Arilla et al. 2020; Luzón et al. 2021) and lithic materials (Villa 2004; Schoville 2014; de la Torre et al. 2018; Bel 2022). Refitting analyses of lithic and bone remains are valuable for evaluating the spatio-temporal integrity of sites (Villa 1982; Hofman 1986; Cziesla 1990; Todd and Standford 1992; Petraglia et al. 1994; Morrow 1996; Morin et al. 2005; Pollarolo et al. 2010; Bargalló et al. 2016; Vaquero et al. 2017; Deschamps and Zilhão 2018; García-Moreno et al. 2023; Falcucci et al. 2024; Sossa et al. 2025; Goder-Goldberger et al. 2025). Other significant contributions include the archaeo-stratigraphic method (Canals et al. 2003; Sañudo et al. 2016; Sossa-Ríos et al. 2022; Fraile-Márquez et al. 2022), fabric analysis of sedimentary and archaeological materials (Bertran et al. 1997; Lenoble and Bertran 2004; McPherron 2005, 2018; Benito-Calvo and de la Torre 2011; Domínguez-Rodrigo et al. 2012; de la Torre and Benito-Calvo 2013; Dibble et al. 2015; García-Moreno et al. 2016) and particle size analysis (Petraglia and Potts 1994; Bertran et al. 2010, 2012; de la Peña et al. 2022).

Recently, geostatistics has become increasingly prominent, particularly point pattern analysis at the intra-site level, which evaluates variables such as material type (Panera et al. 2019; Mendez-Quintas et al. 2019; Cobo-Sánchez 2020; Zilio et al. 2021; Moclán et al. 2023a; Arteaga-Brieba 2024; Merino-Pelaz et al. 2024; Goder-Goldberger et al. 2025), particle size (Zilio et al. 2021; Moclán et al. 2023a; Arteaga-Brieba 2024; Spagnolo et al. 2024), lithic typology (Domínguez-Rodrigo et al. 2014; Alperson-Afil 2017; Panera et al. 2019; Mendez-Quintas et al. 2019; Rabuñal Gayo 2021; Diez-Martín et al. 2021; Moclán et al. 2023a; Rabuñal et al. 2023; Arteaga-Brieba 2024), raw material (Spagnolo et al. 2016; Diez-Martín et al. 2021; Moclán et al. 2023a, b; Arteaga-Brieba 2024), taxa (Peters and van

Kolfschoten 2020; Zilio et al. 2021), taphonomic variables (Giusti and Arzarello 2016; Alperson-Afil 2017; Saladié et al. 2021; Moclán et al. 2023b; Mielgo et al. 2024; Spagnolo et al. 2024), or anthropogenic alterations (Peters and van Kolfschoten 2020; Zilio et al. 2021; Moclán et al. 2023b; Spagnolo et al. 2024; Goder-Goldberger et al. 2025). Spatial statistics have proven crucial in addressing key issues, including formation processes (Mendez-Quintas et al. 2019; Peters and van Kolfschoten 2020; Sánchez-Romero et al. 2020; Zilio et al. 2021; Méndez-Quintas et al. 2022; Rabuñal et al. 2023), palimpsest dissection (Spagnolo et al. 2016; Martín-Perea et al. 2020; Luzón et al. 2021; Martín-Perea 2021; Arteaga-Brieba et al. 2023; Rabuñal et al. 2023; Mielgo et al. 2024), and the reconstruction of anthropogenic behavior (Bargalló et al. 2016; Alperson-Afil 2017; Domínguez-Rodrigo and Cobo-Sánchez 2017; Spagnolo et al. 2019; Cobo-Sánchez 2020; Diez-Martín et al. 2021; Sossa et al. 2025). Ultimately, these analyses reveal the interplay between cultural and natural agents in site formation, and help uncover spatial relationships and underlying structures that could otherwise be imperceptible.

This study investigates the formation processes of Level 3 of Des-Cubierta Cave using spatial statistics to analyze the distribution of geological -clasts- and archaeological materials. Complementary cranial bone refit analyses and spatial assessments of cranial completeness aim to evaluate sedimentary dynamics and their impact on material preservation. The analysis aims to test the following hypotheses:

* Geological and archaeological materials, if originating from different sources, should exhibit distinct spatial patterns.

* Successive rockfalls may have created specific sedimentary structures, which could be detectable in the spatial distribution of the dolomitic clasts.

* The multi-episodic deposition of Unit 3 may have resulted in sedimentary hiatuses, potentially identifiable in the spatial distribution of the remains.

* The accumulation-transit sedimentary context might have influenced the fragmentation and transport of faunal remains, as well as the preservation of the crania, and had a detectable effect on their spatial distribution.

## The site: Des-Cubierta Cave

Des-Cubierta Cave is situated in the Guadarrama Range, a central Iberian mountain range within the Central System. Located on the Calvero de la Higuera hill at 1118 m a.s.l., the cave lies in the Upper Lozoya River Valley, near Piniña del Valle in Madrid region. This tectonic valley, oriented NE-SW, is flanked by mountain ranges reaching 1700–2400
m a.s.l., with Peñalara peak (2428 m a.s.l.) as the highest point. Des-Cubierta Cave is part of the archaeopaleontological sites of Pinilla del Valle, alongside Camino Cave (Alférez et al. 1982; Arsuaga et al. 2010, 2012; Álvarez-Lao et al. 2013), Navalmaíllo Rock Shelter (Márquez et al. 2013; Arriaza et al. 2017; Abrunhosa et al. 2019; Abrunhosa 2020; Moclán et al. 2020, 2021, 2023b, 2024; Blain et al. 2022), and Buena Pinta Cave (Arsuaga et al. 2010; Baquedano et al. 2012, 2016a, b; Laplana et al. 2015, 2016; Blain et al. 2025) (Fig. 1).

Des-Cubierta Cave is located in Upper Cretaceous limestone and dolostone on the right bank of the permanent Lozoya River, within a multilevel karst system characterized by subhorizontal conduits shaped by structural features and the Lozoya's base level. The carbonates hosting the karst reflect a depositional history of shoaling-upward sequences comprising three distinct units: a lower marl unit, less resistant to weathering; a middle unit, consisting of carbonate bars more porous and soluble, where endokarst development is most pronounced; and an upper unit formed by carbonates with heavy cementation, which makes it more resistant to weathering and results in the formation of cliffs and scarps (Pérez-González et al. 2010; Martín-Perea et al. 2025). Although the exact age of the Pinilla del Valle karst is unknown, it is mainly linked to Middle Pleistocene terraces, indicating a minimum age corresponding to that period (Baquedano et al. 2023; Supp. Inf.; Martín-Perea et al. 2025).

Des-Cubierta Cave is a narrow, sub-horizontal, zigzag-
shaped conduit, 80 m long and ranging from 1.5 to 4.5 m
wide. Its recharge zone is located in the southern part of
the Calvero, while the discharge zone is a few meters to the
north in the valley of Arroyo del Hontanar, near the entrance
to Buena Pinta Cave (Fig. 1b). Along the gallery, various

lithostratigraphic units were deposited containing archaeo-
logical and paleontological remains (Fig. 2d) (for a com-
plete description of the units, (see Baquedano et al. 2023;
Supp. Inf.; Martín-Perea et al. 2025).

Unit 3 formed above a speleothem (S1 speleothem) dated to MIS6a (Fig. 2b-d, Supplementary Information 1), likely during MIS4 or early MIS3, a cold and dry period indicated by pollen and microfaunal records (Baquedano et al. 2023; Martín-Perea et al. 2025). This 2-meters thick clast-supported unit primarily comprises angular and subangular cobbles and boulders, with minimal carbonate matrix, deposited during successive rockfall events. It is concentrated in the so-called Monumental sector, the gallery's widest part, and the southern narrowing zone, where gallery narrows from 4.5 to 2 m (Fig. 2a). The underlying paleosurface varies, lying over cryoclastic dolostone cobbles (Unit 4), sandy-clay deposits (Unit 5), or speleothems S1 and S2, dated to 185.35±4.7 ka and 135.7±1.9 ka, respectively (Baquedano et al. 2023; Supp. Inf. p. 86). In the southern narrowing zone, part of speleothem S1 remains in situ, while other fragments are imbricated and markedly verticalized, forming a couvette-shaped structure that acts as a southern boundary within the conduit (Supplementary Information 1). Later, during the latter part of the Upper Pleistocene, an intense erosive process dismantled the gallery roof, impacting the entire Calvero, including the karst system conduits and their fill deposits. This process resulted in the current roofless state of Des-Cubierta Cave. During the Holocene, Unit H, an Ap horizon, was formed, lying on different levels of gallery fill, including Unit 3, which was horizontally affected in the upper parts of the unit and laterally in favor of the sinkhole present in the east wall and lower levels of the conduit (Baquedano et al. 2023; Martín-Perea et al. 2025).

<figure><img src="image_2.png" /><figcaption>**Fig. 1 (a)** Map of the Iberian Peninsula showing the location of the Upper Lozoya Valley (b) Calvero de la Higuera hill with 1 m contour lines and the location of the archaeological sites (modified from Martín-Perea et al. 2025)</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** (a) Photogrammetric model of the main gallery of Des-Cubierta Cave, by Alfonso Dávila (modified from Baquedano et al. 2023) (b) Detail of the study area where Level 3 was deposited, with the spatial analysis window in plan view outlined in white. The S1 speleothem is indicated by a dashed pink background pattern, marking the segment that remains in situ in the narrowing zone. The excavation grid is shown in solid black lines, and UTM coordinate axes are displayed</figcaption></figure>

along the edges. A–A' indicates the limits of the section (c) Section view of the study area, with the spatial analysis window in section view outlined in black, the A-A' markers, and the in situ segment of the S1 speleothem highlighted. The excavation grid is shown as a dashed black line, with UTM coordinate axes along the edges (d) Detailed stratigraphic sections of Monumental Sector in Des-Cubierta Cave, in squares I'43, H'46. Modified from Martín-Perea et al. (2025)

Unit 3 contains an archaeological record consisting mainly of lithic industry and faunal remains, some of them affected by fire (for a complete description see Baquedano et al. 2023). The lithic industry is clearly Mousterian, with an assemblage formed by anvils, hammerstones, cores, flakes and shaped-tools, on raw materials available locally, highlighting the use of quartz, and gneiss for large tools (mainly anvils and hammerstones). The faunal record shows a clear over-representation of cranial elements, with 88% of the NISP belonging to this type of elements (Baquedano et al. 2023; Supp. Inf.). The cranial remains represent at least 35 individuals from species with defensive appendages, including 28 bovines, five cervids, and two rhinoceroses (Baquedano et al. 2023). These crania were deliberately brought into the cave, while the majority of the remaining carcass elements -including the postcranial skeleton, mandibles, and maxillae- was largely absent. The morphology of the anthropized crania in Level 3 aligns with the expected

outcome of extracting nutritional resources from the cranial region, excluding the brain (e.g., jaw meat, tongue, and eyes). This pattern is consistent with results from an experimental butchering of cattle heads (Baquedano et al. 2023; Supp. Inf., Álvarez-Fernández et al. 2025), which left only the corneal appendages, neurocranium, and splanchnocranium, with mandibles and maxillae removed.

The bone remains show a high fragmentation, being compatible with the mechanical action produced by falling boulders from the gallery roof and walls (Oliver 1989; Fisher 1995). Taphonomic analysis shows a bone assemblage with hardly any modifications by carnivores and with postdepositional modifications consistent with damp conditions with low-energy runoff typical of a karst environment, such as the presence of concretion on the bone surface, manganese oxide, and low degrees of rolling and polishing (Baquedano et al. 2023; Supp. Inf.). The discarding of other accumulating agents of these archaeological remains, such
as gravitational processes, water transport, biotic agents such as carnivores or anthropic activities with subsistence purposes, leads Baquedano et al. (2023) to propose an interpretation related to the symbolic world of the Neanderthals. The deliberate selection of species with cranial appendages for carcass processing outside the cave, followed by the intentional introduction of modified crania into this narrow gallery, and the repeated execution of these non-subsistence practices—evidenced by the accumulation of at least 35 skulls across multiple depositional episodes—makes this site a key location for investigating the behavioral complexity and symbolic practices of *Homo neanderthalensis*.

## Materials

The materials analyzed in this study consist of the point pattern materials from Level 3 of Des-Cubierta Cave consisting of dolomitic clasts (n=3047) and archaeological materials (n=6144), including bones (n=4780) and lithic artifacts (n=1719). All materials were recovered during field seasons conducted between 2009 and 2022. The materials from the Unit 3 area laterally affected by Unit H were excluded from the analysis.

During the fieldwork, the archaeological surface was divided into a 1×1 m grid. Archaeological and geological materials were mapped using a total station (TS Leica TCRP 1205 R400 model), capturing their spatial coordinates (x, y, z) and drawing them in map at a 1:10 scale. The grid followed the orientation previously established at the nearby site of Buena Pinta Cave (Fig. 2b), with a local north slightly offset from the geographical north. However, for clarity and consistency in the presentation of spatial analyses, the figures in this article use the axes of the UTM coordinate system, allowing direct reference to the original measurement units. The orientation and slope of the materials were categorized without the use of a compass or clinometer. Clasts larger than 200 mm and bones larger than 20 mm were recorded, while lithic artifacts were coordinated regardless of size. Archaeological materials recovered through sieving or washing, lacking spatial coordinates, were excluded from this study. Initially, bone remains and lithic artifacts were studied together as a single set of archaeological materials, distinct from the clasts or geological materials. Subsequently, the analysis focused on bone remains, particularly cranial remains, which proved to be as the predominant material in Level 3 (Baquedano et al. 2023).

The excavation surface with nearly the entire paleosurface on which Level 3 of Des-Cubierta Cave was deposited, except for control areas in squares I'44–I'45 and K'39, was documented through a photogrammetric model. Photographs were taken with a Nikon D500 camera. The photogrammetric model was carried out using Agisoft Photoscan 1.4.3 software.

## Methods

### Geostatistical analysis: point pattern analysis

Point pattern analysis was applied to examine the spatial organization of geological and archaeological materials from Level 3 of Des-Cubierta Cave, aiming to identify the formation processes underlying the observed distributions. Following Baddeley et al. (2016), we classified patterns as either unmarked (single-category samples, such as geological or archaeological materials) or marked (multi-category patterns, such as material type or clast granulometry). Geological clasts, ranging from 200 to 905 mm, were classified according to Blott and Pye's (Blott 2012) scale into small (200–256 mm, n=1936), medium (256–512 mm, n=1077), and large (512–905 mm, n=34) boulders, forming the basis for subsequent spatial analyses. Spatial analysis in archaeological contexts often encounters challenges such as high density, spatial inhomogeneity, and visualization complexities, particularly in three-dimensional environments. As current methods for correcting inhomogeneity in 3D remain limited, analyses were performed on two-dimensional projections—longitudinal (yz) and plan (xy)—given the ~2 m thickness of Level 3.

All analyses were performed in R (R Core Team 2024) using the spatstat package (Baddeley and Turner 2005; Baddeley et al. 2016). Data, code, and visual outputs are openly available in our research compendium (https://osf.io/frzts/, DOI https://doi.org/10.17605/OSF.IO/FRZTS), following the reproducible research framework outlined by Marwick (2017).

Spatial analysis focused on two complementary aspects: intensity and spatial correlation (Bevan 2020). Intensity describes variations in point density, while correlation evaluates spatial interactions between points. Departures from Complete Spatial Randomness (CSR) were tested through quadrat-based ($\chi^2$, Freeman-Tukey, modified Neyman) and nearest-neighbor-based approaches (Clark-Evans and Hopkins-Skellam). Edge effects were corrected using the cumulative distribution function method. The Clark-Evans test calculates the ratio (R) between observed and expected nearest-neighbor distances ($R<1$ = clustering; $R>1$ = dispersion), while the Hopkins-Skellam calculates the A index (A=0 for CSR, A<1 clustering, A>1 dispersion) and is less sensitive to spatial inhomogeneity. Robustness was evaluated through Monte Carlo simulations (39 simulations for Clark-Evans, 999 for Hopkins-Skellam). The rejection of CSR prompted the use of additional methods addressing spatial inhomogeneity and correlation.

Pattern intensity was estimated non-parametrically using Kernel Density Estimation (KDE) (Baddeley et al. 2016). Bandwidth ($\sigma$) was selected by likelihood cross-validation, with smoothing adjustments applied where necessary; detailed $\sigma$ values are listed in Supplementary Information
2. Edge effect corrections were included for all maps. The scan test was applied to identify high-intensity areas showing statistically significant density compared to the surrounding regions, effectively dividing the spatial window into zones with (TRUE) and without (FALSE) significant high-density clusters (p=0.05). Outputs from the KDE and scan test, together with density lines (z-coordinate for vertical and x/y-coordinates for horizontal axes), were combined to provide an integrated visualization of spatial structure. To extend the intensity analysis to multitype point patterns, we applied the relative risk function—optimized by cross-validation to select the smoothing bandwidth—which estimates the spatial variation in the probability of occurrence of each point type relative to others (Bevan 2012; Smith et al. 2015; Cobo-Sánchez 2020; Diez-Martín et al. 2021; Moclán et al. 2023b; Arteaga-Brieba 2024). The resulting heat maps display the relative probability of each material type, with tolerance contours (Hazelton and Davies 2009; Cobo-Sánchez 2020) highlighting areas where probabilities deviate significantly from the mean (p=0.05), as determined by Monte Carlo simulations (n=19). Dominant risk maps were used to identify zones where one material type predominates, revealing regions dominated by geological clasts in the central sector. These results were imported into QGIS 3.22.14 to intersect Level 3 archaeological materials with low-intensity zones for subsequent spatial analysis. As part of the multitype analysis, we also applied a segregation test to evaluate whether the spatial distributions of different types significantly deviated from a random labeling scenario. The test uses a t-statistic to quantify segregation strength, with significance assessed at p=0.05 through 39 Monte Carlo simulations based on random relabeling while keeping point locations fixed.

To assess the spatial correlation of the point pattern, we applied both unmarked and marked versions of the inhomogeneous *L* function, a transformed version of Ripley's *K* function (Ripley 1979), with edge-effect corrections. These functions evaluate departures from spatial randomness—CSR in unmarked patterns or its inhomogeneous equivalent—by measuring clustering or dispersion across multiple spatial scales. For unmarked patterns, we used the *Linhom* function, and for marked (multitype) patterns, the cross-type *Lij(r)* function, which tests for spatial interaction between types under the null hypothesis of random labeling. In both cases, statistical significance was assessed through 39 Monte Carlo simulations, with intensity estimated via likelihood cross-validation. Because of the high intensity and spatial inhomogeneity of the patterns, *L*-based methods may not fully capture local relationships. We therefore applied nearest-neighbor-based methods (*nncount* and *nnequal*), which focus on point proximity rather than intensity, calculating the proportion of neighboring points of each type and comparing these to expectations under random mixing. Statistical significance was again evaluated through 39

Monte Carlo simulations based on random relabeling (Cobo-Sánchez 2020; Diez-Martín et al. 2021).

## Cranial bone refitting

Bone refitting was conducted to evaluate whether rockfall events in Level 3 of Des-Cubierta Cave acted as agents of fragmentation, transport, or both, influenced by the morphology of the gallery and sedimentary structures. Bone refitting was conducted within the context of the restoration work and the archaeofaunal study of the Unit 3 assemblage. The restoration team's efforts were focused on the crania, with the goal of reconstructing them to restore their legibility, thereby enabling taxonomic and morphometric analysis, as well as preserving each specimen (Galindo-Pellicena et al. 2019; Baquedano et al. 2023). Postcranial elements were not considered. In addition, smaller fragments, trabecular bones, those with fracture morphologies unsuitable for refitting, and most burned fragments were excluded due to their low potential for refitting. Given the high fragmentation of the crania, extraction sometimes required facing techniques, which resulted in some elements lacking individual spatial coordinates. In these cases, a single XYZ coordinate was recorded for the entire group of associated fragments.

After refitting, spatial coordinates (x, y, z) for each fragment were recorded in a database, along with a unique refit ID. Using QGIS 3.22's Points to Path tool, refitting lines were generated to connect fragments of the same refit, with the z-coordinate determining the sequence of union, assuming higher elevation fragments were in their original position prior to fragmentation. For multi-fragment refits, the restoration order served as a supplementary guide. Metrics such as line length, orientation, and vertical displacement were calculated in QGIS, with further analysis conducted in Rstudio (RStudio 2024).

Spatial analysis of refitting lines incorporated intensity analysis and circular statistics, visualized through rose diagrams and kernel density plots for circular data. Analyses were conducted in R (R Core Team 2024) using the `circular` (Agostinelli and Lund Agostinelli 2023), `CircStats` (Lund and Agostinelli 2025), `ggplot` (Wickham 2016), and `spatstat` (Baddeley et al. 2016). Both the total set of refitting lines (n=124) and those ≥ 20 cm in length (n=28) were analyzed. Key parameters—including mean vector, concentration (kappa), circular variance, circular standard deviation, and standard error—were calculated to evaluate potential anisotropy. Confidence intervals for the mean direction were estimated, and directionality and anisotropy were tested using Rayleigh, Kuiper, and Watson's U² tests. All circular plots and calculations were based on degrees, with 0° (or 360°) aligned to geographic North. These methods follow previous applications in archaeological contexts (Cobo-Sánchez et al. 2014; García-Moreno et al. 2023).
## Preservation and completeness of crania

The 35 crania deposited in Level 3 exhibit varying degrees of integrity or completeness. This study aimed to determine whether these differences correlated with spatial location and structural features within the gallery. Based on the archaeological assemblage and on experimental butchering results (Baquedano et al. 2023), all crania were assumed to have been deposited in a standardized morphological state— without mandibles and maxillae. Crania were categorized into four completeness grades following restoration:

1. Completeness Grade I: Isolated skull fragments preserving only the cornual protuberances (horns or antlers).

2. Completeness Grade II: Skull fragments with horns/ antlers and some neurocranial bones, without splanchnocranial elements.

3. Completeness Grade III: Includes elements from Grade II, with a largely complete neurocranium and, in some cases, partial splanchnocranial bones.

4. Completeness Grade IV: Almost complete skulls, preserving both the neurocranium and splanchnocranium, with only the mandible and maxilla missing.

The spatial distribution of the crania was analyzed by plotting the digitized field drawings of the crania, categorized by completeness grade, in plan view. Similarly, the points corresponding to the primary fragments were plotted in yz section view, also classified according to their completeness. In both cases, these distributions were overlaid onto the cone-shaped sedimentary structure identified through spatial analysis.

## Results

### Geostatistical analysis

### Comparison between geological and archaeological pattern

The geological and archaeological point patterns showed distinct spatial characteristics. Archaeological materials exhibited higher intensity (180.87 points/m² in plan;

289.42 points/m² in section) compared to geological clasts (88.15 points/m² in plan; 144.32 points/m² in section). Statistical tests—including the χ² test, which indicated spatial inhomogeneity—suggest deviations from complete spatial randomness. Clark-Evans and Hopkins-Skellam tests also returned values consistent with aggregation; however, these results should be interpreted with caution, as such methods are sensitive to intensity variation and may reflect inhomogeneity rather than true spatial interaction, with these tendencies more pronounced in the archaeological materials (Table 1).

KDE revealed differences in spatial distribution, both in section and plan. Geological materials displayed a more evenly distribution, with higher intensity in the central-north zone (4530742<y<4530744) (Fig. 3e and 4e), corresponding to the widest section of the gallery. In contrast, archaeological materials were more spatially concentrated in the narrowing southern zone (4530738<y<4530742) (Fig. 3d and 4d). The scan test confirmed spatial segregation, identifying distinct high-intensity areas for each material type (Fig. 3f and 4f).

Further analyses using the inhomogeneous L function revealed contrasting trends between geological and archaeological materials, as well as differences between plan and section views (Supplementary Information 3a, b, f, g). However, these patterns lacked strong statistical significance, and the observed segregation by zones was better validated through other methods. The results from the relative risk function with tolerance contours, the dominant risk function, and the segregation test (T=909.70, *p*=0.05 for section; T=809.40, *p*=0.05 for plan) confirmed that the spatial segregation of geological and archaeological materials is statistically significant (Figs. 3g-i and 4g-i).

Beyond the north-south segregation identified in the intensity analyses of the unmarked patterns, a distinct area emerges at the base of the central zone (4530742<y<4530745, 1108<z<1108.50), where the probability of the geological type is clearly predominant (Figs. 3h-i and 4h-i). Aggregation and segregation were further validated by nearest-neighbor methods (nncount, nnequal), showing clustering of like materials and segregation between types (Supplementary Information 3d, e, i, j).

<table><caption>**Table 1** Results of Spatial tests applied to geological and archaeological materials in plan and section views. The tests evaluate Spatial patterns: χ² tests identify deviations from Spatial randomness (i.e., inhomogeneity in intensity), the Clark-Evans R index assesses point arrangement, and the Hopkins-Skellam A index quantifies clustering tendencies. The results suggest significant departures from Spatial randomness for both geological and archaeological materials, primarily due to inhomogeneous intensity. These effects are more pronounced in the archaeological datasets, particularly in the plan view</caption><thead><tr><th rowspan="2"></th><th>chi0</th><th>chi1</th><th>chi2</th><th>chi3</th><th colspan="2">Clark-Evans</th><th colspan="2">Hopkins-Skellam</th></tr><tr><th><i>p</i>-value</th><th><i>p</i>-value</th><th><i>p</i>-value</th><th><i>p</i>-value</th><th><i>R</i> Index</th><th><i>p</i>-value</th><th>A Index</th><th><i>p</i>-value</th></tr></thead><tbody><tr><td>Geological plan</td><td>1.2126E-174</td><td>6.9804E-249</td><td>0</td><td>0.05</td><td>0.90479</td><td>0.025</td><td>0.21274</td><td>0.001</td></tr><tr><td>Geological section</td><td>1.2126E-174</td><td>6.9804E-249</td><td>0</td><td>0.05</td><td>0.90499</td><td>0.025</td><td>0.21274</td><td>0.001</td></tr><tr><td>Archaeological plan</td><td>0</td><td>0</td><td>0</td><td>0.05</td><td>0.81009</td><td>0.025</td><td>0.11577</td><td>0.001</td></tr><tr><td>Archaeological section</td><td>0</td><td>0</td><td>0</td><td>0.05</td><td>0.81009</td><td>0.025</td><td>0.11779</td><td>0.001</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Section view. Geological materials: (a) Distribution map (d) KDE map (g) Relative risk map. Archaeological materials: (b) Distribution map (e) KDE map (h) Relative risk map. The relative risk maps include a tolerance contour (white line) indicating areas where probabilities significantly deviate from the mean. Combined geological and archaeological materials: (c) Distribution map (f) Scan test results shown as a superimposed map of areas with significant</figcaption></figure>

(p≤0.05) high density (TRUE, colored) and non-significant (FALSE, white) zones for both materials, highlighting the absence of overlap in high-density areas (i) Dominant risk map. High-density areas of geological and archaeological materials are spatially segregated, with no overlap between zones of maximum concentration, emphasizing the distinct depositional domains of clastic and archaeological components

## First rockfall and cave opening

Intensity analysis identified a high-density geological zone overlapping a low-density archaeological zone. This area, spanning coordinates 4530742.53°<y<4530744.85 and 1108.02°<z<1108.51 m a.s.l., contained 141 records: 140 geological clasts and a single bone (Fig. 5a-c). The only bone exhibited characteristics consistent with percolation or gravitational movement -namely, its small size (35×13×8 mm) and trabecular structure. Geological materials comprise 1 large boulder, 60 medium boulders, and 79 small boulders.

The spatial distribution of these clasts indicated a cone-shaped structure in section and a fan-shaped structure in plan, with higher density near the east wall, gradually decreasing westward. KDE, scan tests (Fig. 5.d-g) and the density graph on the x-axis (Supplementary Information 4.a) revealed spatial organization by size, with medium-sized boulders forming a central high-intensity focus and small boulders showing a surrounding ring-shaped distribution. These findings suggest the formation of a rockfall cone in the Monumental sector, representing the earliest rockfall event and potentially marking the creation of a skylight, located in the area 432,078°x<43,207,080, and 4,530,742°<y<4,530,744.

## Sedimentary structures: cone-shaped rockfall

Regarding the analysis of intensity of clasts divided by size, the KDE of medium-sized boulders revealed several underlying structures not previously detected in fieldwork or distribution maps. In section, a high-intensity cone-shaped area was observed between 4,530,742°<y<4,530,744 and 1108°<z<1109 m a.s.l., corresponding to the area where the first phase of cone formation had been identified. Above this, a decrease in intensity was observed, followed by another conical shape at 4,530,741°<y<4,530,745, with elevations up to 1109.75 m a.s.l. Finally, a high-intensity line appeared at 4,530,738°<y<4,530,746, sloping to the south (Fig. 6e and h). These conical morphologies were interpreted as successive accretion episodes, reflecting multiple phases of cone growth due to repeated rockfalls.

Plan views showed a lateral cone expanding westward from the apex near the east wall, in the widest area of the gallery (Fig. 7e and h.). A smaller, circular high-intensity focus was observed in the southern zone, the narrowest part of the gallery. Overlaying the photogrammetric model and the Mesozoic sedimentary unit map with intensity maps revealed that this hotspot coincided with a fracture in the S1 speleothem and the junction between the middle and lower marl units (Fig. 7j.). The S1 speleothem remained in situ until approximately
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Plan view. **Geological materials:** (a) Distribution map (d) KDE map (g) Relative risk map. **Archaeological materials:** (b) Distribution map, (e) KDE map (h) Relative risk map. The relative risk maps include a tolerance contour (white line) indicating areas where probabilities significantly deviate from the mean. **Combined geological and archaeological materials:** (c) Distribution map (f) Scan test results shown as a superimposed map of areas with significant</figcaption></figure>

(p≤0.05) high density (TRUE, colored) and non-significant (FALSE, white) zones for both materials, highlighting the absence of overlap in high-density areas. High-density areas of geological materials are concentrated in the widest part of the gallery, while those of archaeological materials are focused in the narrowing section highlighting the spatial segregation of depositional areas
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Area devoid of archaeological materials at the base of Level 3. (a) Relative risk map of archaeological materials in section view, highlighting the low-density area at the base of the level (b) Distribution map of geological and archaeological materials with the same area marked. Boulders within this area are divided by size in section</figcaption></figure>

view (c) and plan view (d). (e) KDE map of medium boulders (f) KDE map of small boulders in this area (g) Superimposed map of significant ($p \le 0.05$) high-density zones (TRUE, colored) and non-significant zones (FALSE, white) for both boulder sizes, illustrating size-based sorting

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Section view. Distribution map of large (a), medium (b), and small boulders (c). KDE of large (d), medium (e), and small boulders (f). Scan test $p<0.05$ of large (g), medium (h), and small boulders (i). The medium-boulder density patterns reveal a distinct conical morphology,</figcaption></figure>

in contrast to the more dispersed arrangements of large and small boulders, defining the main sedimentary structure of the deposit. Red frames highlight the key interpretative patterns

4530738.50 y and 1108.70 z, where it fractured, creating a 40 cm drop and a cuvette-shaped depression surrounded by imbricated fragments (Fig. 2b, c. Supplementary Information 1). This localized structural instability likely contributed to the higher density of medium-sized boulders in this area.

Small boulders showed a more evenly spread pattern, with a high-intensity area at $4,530,742 < y < 4,530,744$ and elevations between 1108 and 1110 m a.s.l. (Figure 6f and i). In the southern zone, three high-intensity lines, separated by small hiatuses, may indicate multiple depositional events.
In plan view, their distribution aligned with that of medium-sized boulders, both displaying circular patterns near the S1 fracture, though the smaller boulders exhibited a more linear trend in the central gallery (Fig. 7c, f and i).

Large boulders exhibited two primary high-intensity zones: one at medium elevations (1109–1110 m a.s.l.) in the southern gallery and another around 1110 m a.s.l. in the northern half (Fig. 6d and g). In plan, they formed an annular pattern along the cone's periphery, consistent with medium-sized boulders (Fig. 7d and g).

Vertical distribution analysis confirmed size-based sorting, with density peaks for large boulders at medium elevations (z≈1109.50 m a.s.l.) and multiple peaks for small boulders corresponding to decreases in medium boulders (Supplementary Information 3b), reflecting the sedimentary dynamics of repeated rockfalls.

## Bone refitting

### Number and characteristics

A total of 60 refits involving 194 bone fragments have been identified. The majority (98.34%) correspond to mechanical refits, where two or more fragments from the same fractured bone joined, likely as a result of collapse processes within the gallery leading to dry/diagenetic fractures. Only one anatomical refit (Refit19) was observed, involving bilateral maxilla fragments. Most refits (63.33%) involved two fragments, while one exceptional case (Refit31) included 19 fragments (Table 2), restoring a cranium with significant portions of the frontal and occipital bones and nearly complete horns (Fig. 8d.3).

<figure><img src="image_8.png" /><figcaption>**Fig. 7 Plan view.** Distribution maps of large (a), medium (b), and small boulders (c). KDE maps of large (d), medium (e), and small boulders (f). Scan test results ($p \le 0.05$) for large (g), medium (h), and small boulders (i). Raster and contour map of medium boulders' KDE overlaid on the gallery's photogrammetry, along with the Mesozoic sedimentary units and speleothems (j). (modified from Martín-Perea et al. 2025). The medium-boulder density pattern reveals in plan view</figcaption></figure>

the conical morphology of the deposit, with a prominent high-intensity area in the widest zone of the gallery, aligned with a scarp of the hill at the boundary between the middle (more porous) and upper (more resistant) carbonate units. A smaller hot spot is also visible in the narrowing zone, corresponding to the fracture of the S1 speleothem and the boundary between the middle and lower marl unit (less resistant to weathering). Red frames highlight the key interpretative patterns
<table><caption>**Table 2** Summary of bone refitting results from level 3 of Des-Cubierta Cave. The table shows the number of elements involved in each refit, the total number of refits for each category, and their corresponding percentages</caption><thead><tr><td>Number of refitting elements</td><td>Number of refits</td><td>%</td></tr></thead><tbody><tr><td>2</td><td>38</td><td>63.33</td></tr><tr><td>3</td><td>8</td><td>13.33</td></tr><tr><td>4</td><td>6</td><td>10</td></tr><tr><td>5</td><td>1</td><td>1.66</td></tr><tr><td>6</td><td>2</td><td>3.33</td></tr><tr><td>7</td><td>2</td><td>3.33</td></tr><tr><td>9</td><td>1</td><td>1.66</td></tr><tr><td>11</td><td>1</td><td>1.66</td></tr><tr><td>19</td><td>1</td><td>1.66</td></tr></tbody></table>

## Refitting lines

The spatial distribution of refitting fragments was analyzed in plan, generating 124 refitting lines (Fig. 9a., 9c. and Supplementary Information 5), with the descriptive statistics detailed in Table 3. More than half (54.84%) of these lines span a horizontal distance of up to 10 cm, and 94.35% are ≤50 cm. Only 7 refits exceed 50 cm, with 4 surpassing 1 m; one extends as far as 2.17 m. Refit3 consists of four fragments from the skull's basal area: the right and left petrous bones, the basioccipital with the right occipital condyle, and the left occipital condyle. Of these, two fragments share spatial information, thus resulting in only two refitting lines: one with a distance of 1.55 m (right petrous bone to basioccipital) and another of 2.17 m (basioccipital to left occipital condyle). These refits were located at the upper part of the level (z-values between 1109.65 and 1110.01 m a.s.l.), which was heavily affected by diagenetic processes such as edaphization and erosion (Martín-Perea et al. 2025). Similarly, Refit55 involves the right nasal and the frontal bone of Crania 6, with a 1.20-meter distance between the two fragments. This refit, located in the high zone of the level, also shows the impact of these post-depositional processes. The last long-refit, Refit 24 (Fig. 8b) joins the right horn fragments of Crania 14, spanning 1.22 m, and is located in the gallery's narrowing zone, near the fractured S1 speleothem.

The vertical distances of the refits show a mean of 5 cm, with a maximum distance of 48 cm. Only three refits exceed 40 cm in vertical distance: Refit24, Refit28 (two fragments of deer antler), and Refit58 (two antler fragments) (for more details see Table 3).

KDE of horizontal refitting lines revealed hotspots near the gallery's narrowing zone and along the west wall (Fig. 9b), while vertical KDE maps showed high intensity between 1109<z<1110 m a.s.l. (Fig. 9d), with the S1 speleothem fracture acting as a boundary to the south.

## Orientation of refits

The rose diagram and kernel density plot of the refit orientations (Fig. 9e and f) suggest a fairly consistent directional trend, with a slight preference for a southeast direction. When considering only the refits with distances greater than 20 cm (n=28), the trend is less pronounced, slightly favoring the E-SE direction (Fig. 9g and h). However, the dispersion measures (Table 4) and statistical tests (Rayleigh, Kuiper, and Watson's U2) (results in Table 5) indicate that the data are not tightly clustered, and the p-values suggest that uniformity cannot be ruled out, meaning the trend towards a southeast orientation may not be statistically significant. This orientation may be linked to the general growth direction of the accretion cone but does not appear to be strongly influenced by it.

## Preservation and integrity of crania

The analysis of crania completeness in Level 3 revealed that 20 out of 35 identified individuals (57.14%) exhibited the lowest degree of preservation—Grade 1. This category included 18 bovines and 2 red deers. Grade 2 included 8 crania, 7 of which were bovine and 1 cervid. Grade 3 included 6 crania: 3 bovids, 2 red deers, and 1 rhinoceros. Lastly, Grade 4 included only one rhinoceros' crania, notable for its high degree of completeness and anatomical integrity (Fig. 10 and Supplementary Information 7c).

The spatial distribution of crania across the different completeness grades was analyzed both in plan and section (Fig. 11 and Supplementary Information 6). Interestingly, in the apical zone of the cone-shaped deposit, no crania are present; they generally tend to be located on the periphery of the cone, toward the west part of the gallery. The most poorly preserved specimens (Grade 1) were primarily located in the southern half of the study area, particularly near the west wall. These specimens were often found at higher elevations within Level 3, where weathering and edaphization processes were intensified (Martín-Perea et al. 2025). Five of these crania (Cr. 16, 28, 29, 32, and 34) were located in the southern periphery of the cone, near or within the narrowing zone of the gallery, an area characterized by fragmented and vertically aligned S1 speleothems. Crania 31, also with grade 1 of completeness, was found to the north of the cone at lower elevations. The poor preservation of this specimen is likely related to the presence of runoff or water passage in that area.

In contrast, crania with Grade 2 integrity were consistently located at medium elevations. Six of these specimens were found in the southern half of the gallery, while one (Cr. 8) was situated in the central zone and another (Cr. 20) in the northern area. Crania with Grade 3 completeness were
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Selected examples of refits, with coordinated fragments outlined in white. (a) Refit 4 (b) Refit 24 (c) Refit 20 (d) Refit 31 (d.1) Lateral view of a large boulder above the right horn of Crania 24 -pointed out with a red arrow- with a porphyry beside it, pointed out</figcaption></figure>

with the green arrow (d.2) Top view of the right horn and the porphyry after removal of the boulder, revealing numerous cranial fragments (d.3) Crania 24 after restoration

predominantly located in the central and northern parts of the gallery, aligning with the cone's morphology: distributed along the cone's periphery in plan view and following a cone-shaped arrangement in section view.

The only individual categorized with Grade 4 completeness (Cr. 22) was found at the base of Level 3, in an area where the initial rockfalls were identified, directly above the cone's earliest phase. This *Stephanorhinus hemitoechus* crania, exhibited exceptional integrity, with all cranial bones present, except for one maxilla (Fig. 10 and Supplementary Information 7c). The preservation of this specimen allowed the identification of numerous cut marks on the zygomatic bone and basisphenoid, as well as clear evidence

of anthropic bone breakage (Baquedano et al. 2023). The anatomical elements missing from the specimen are attributed to anthropic activities rather than post-depositional processes. The second *Stephanorhinus hemitoechus* crania from Level 3 (Cr. 27) is also associated with the initial paleosurface formed after the first rockfalls (Supplementary Information 7a, 7b). This specimen showed poorer integrity compared to Cr. 22, lacking the zygomatic arches and part of the nasal bone. This crania was found in an inverted position (Supplementary Information 7d) south of the cone. This positioning, along with its poor preservation, suggests it may have been subject to post-depositional repositioning due to rockfall processes.
<figure><img src="image_10.png" /><figcaption>**Fig. 9 Plan view.** (a) Line refitting pattern (b) KDE of the line pattern, showing a pronounced high-intensity area along the west wall in the narrowing zone. **Section view.** (c) Line refitting pattern (d) KDE of the line pattern. **Orientation analysis.** (e) Rose diagrams displaying the directions of all refitting lines and (g) those longer than 20 cm.</figcaption></figure>

Kernel density estimates for circular data are represented in pie charts, where points are distributed according to vector angles, and an outer line indicates the peaks in directions with the highest density, for all refitting lines (f) and longer lines (h)

<table><caption>**Table 3** Descriptive statistics for the length (m) of bone refits from level 3 of Des-Cubierta Cave, presented for plan and vertical distances, showing greater plan than vertical distances, indicating limited vertical displacement</caption><thead><tr><th></th><th>Mean</th><th>Std. Error</th><th>Median</th><th>Mode</th><th>Standar Dev.</th><th>Variance</th><th>Min</th><th>Max</th><th>Range</th></tr></thead><tbody><tr><td>Plan distance</td><td>0.179</td><td>0.026</td><td>0.086</td><td>0.032</td><td>0.294</td><td>0.087</td><td>0.010</td><td>2.174</td><td>2.164</td></tr><tr><td>Vertical distance</td><td>0.057</td><td>0.008</td><td>0.020</td><td>0.000</td><td>0.096</td><td>0.009</td><td>0.000</td><td>0.480</td><td>0.480</td></tr></tbody></table>

<table><caption>**Table 4** Descriptive statistics for circular data of refitting lines from level 3, showing weak directional trends with mean vectors near 101.7° (all lines) and 82.6° (long lines, ≥ 0.2 m), high dispersion, and slightly stronger clustering for longer lines</caption><thead><tr><th></th><th>Total lines</th><th>Long lines</th></tr></thead><tbody><tr><td>Number of Observations</td><td>124</td><td>28</td></tr><tr><td>Mean Vector (μ)</td><td>101.717°</td><td>82.566°</td></tr><tr><td>Length of Mean Vector</td><td>0.155</td><td>0.111</td></tr><tr><td>Median</td><td>103.676°</td><td>98.013</td></tr><tr><td>Concentration</td><td>0.184</td><td>0.263</td></tr><tr><td>Circular Variance</td><td>0.845</td><td>0.889</td></tr><tr><td>Circular Standard Deviation</td><td>110.687°</td><td>120.226°</td></tr><tr><td>Standard Error of Mean</td><td>9.94°</td><td>22.720°</td></tr><tr><td>95% Confidence Interval (-/+)</td><td>101.701°<br>101.732°</td><td>82.497°<br>82.635°</td></tr><tr><td>99% Confidence Interval (-/+)</td><td>101.690°<br>101.744°</td><td>82.447°<br>82.684°</td></tr></tbody></table>

<table><caption>**Table 5** Results of statistical tests for uniformity applied to circular data of refitting lines from level 3 of Des-Cubierta Cave. Marginal directional trends are suggested for all lines by the Rayleigh test (p=0.0514), while kuiper's and watson's tests show no significant deviations</caption><thead><tr><th></th><th></th><th>Rayleigh Test</th><th>Kuiper's Test</th><th>Watson's (U²)</th></tr></thead><tbody><tr><td rowspan="2">Total lines</td><td>Statistic</td><td>0.1547</td><td>1.5285</td><td>0.175</td></tr><tr><td>p-value</td><td>0.0514</td><td>>0.15</td><td>0.05&lt;p&lt;0.10</td></tr><tr><td rowspan="2">Long lines</td><td>Statistic</td><td>0.1106</td><td>0.9775</td><td>0.0415</td></tr><tr><td>p-value</td><td>0.7135</td><td>>0.15</td><td>>0.10</td></tr></tbody></table>

## Discussion

The application of geostatistics to the geological and archaeological materials from Level 3 of Des-Cubierta Cave has offered valuable insights into its formation processes,

<figure><img src="image_11.png" /></figure>
particularly the dominant depositional mechanism—a recurring rockfall process—and its implications for the cave's paleosurface and the preservation of cranial remains.

## Conic sedimentary structure: causes and morphology

Geological studies had framed the formation of Unit 3 as primarily a result of successive rockfalls from the cave roof, occurring during a cold period, likely during MIS4 or early MIS3, in a mountainous environment at altitudes above 1100 m a.s.l. (Martín-Perea et al. 2025). Rockfall processes are typical in karst systems, particularly during vadose phases, where they modify existing morphologies and generate breakdown deposits. These deposits consist of poorly sorted angular clasts, ranging from sands to boulders (Ford and Williams 1991; Sasowsky 2007). Causes of rockfall include loss of underlying support, removal of overlying rock layers, ice wedging (cryoclastic processes), mineral wedging (e.g., gypsum), and seismic activity (White and White 2003; White 2007; Sasowsky 2007). At Des-Cubierta Cave, based on available environmental data, cryoclasticity is likely the dominant cause of rockfalls. The faunal and pollen records support this conclusion, indicating a cold, arid environment dominated by open landscapes with scattered forest patches. Micromammal remains are primarily composed of voles (Microtus arvalis, Microtus agrestis, Microtus gr. lusitanicus-duodecimcostatus), with field mice (Apodemus gr. A. sylvaticus-flavicollis), hedgehogs (Erinaceus europaeus), and toads (Bufo gr. bufo spinosus) also present, indicating an absence of thermophilic species and the presence of cold-adapted taxa such as steppe pika (Ochotona cf. pusilla) and Iberian frog (Rana iberica) (Baquedano et al. 2023). The macrofaunal record reflects a similar preference for open grasslands, with steppe bison (Bison priscus) predominating in the archaeo-faunal assemblage, and lesser occurrences of steppe rhinoceros (Stephanorhinus hemitoechus) and horse (Equus ferus). The presence of mixed feeders such as aurochs (Bos primigenius) and red deer (Cervus elaphus) suggests the existence of limited arboreal patches in the surrounding landscape, although faunal representation may be biased by anthropic hunting selection. The pollen data further corroborate these findings, showing low taxonomic diversity dominated by herbaceous species, along with Pinus and Juniperus, consistent with a steppe landscape interspersed with cold-adapted, sparse woodlands (Baquedano et al. 2023). Taken together, these environmental records indicate that Level 3 formed in a periglacial, high-altitude mountainous environment under extreme climatic conditions.

The factors influencing the fracturing and detachment of bedrock fragments are multifaceted and depend on several conditions. The size and shape of debris are determined by

the thickness of the bedrock, the presence of partitions along bedding planes, the density of fracturing, and the configuration of dissolution surfaces before rock breakage (White and White 2003). These factors must be considered alongside the sedimentological environment in which karstification occurred in the Calvero de la Higuera, with the three distinct carbonate units showing varying resistance to weathering and endokarst processes (Pérez-González et al. 2010; Martín-Perea et al. 2025). At Des-Cubierta Cave, one of the scarps formed by the hardest upper unit is located across the Monumental sector, where the gallery reaches its maximum width (4.30 m). Here, the initial stages of cone formation are evidenced by a high-intensity area of boulders without archaeological material. This process evolved through successive rockfalls, with the addition of clastic materials. The wider space in this section, along with the contact between the middle and upper units (Fig. 7j), suggests that this area was the primary source of rockfall, likely initiating roof collapse and the creation of a skylight, which further destabilized the area and promoted continued rockfalls.

The emergence of this vertical opening to the exterior due to a collapse doline significantly altered the gallery conditions. This skylight allowed light to penetrate the cave, improved ventilation, and created a potential access point for Neanderthals. Notably, no anthropogenic materials were found in the gallery prior to the formation of the cone's initial phase. However, immediately above this first phase, archaeological materials characteristic of Level 3 begin to appear, including the Mousterian lithic industry and the large ungulate skulls, which persist throughout the level. These findings suggest that the skylight could have facilitated Neanderthal access to the gallery, serving as a vertical entry point into the cave. The possibility of the skylight functioning as an anthropic access, however, is not unequivocally established. The height of the gallery ceilings during this period remains uncertain due to the erosive processes that dismantled the roof. Currently, the maximum preserved wall height in the Monumental sector is approximately two meters above the base of Level 3 and 1.5 m above the apex of the cone's initial phase, a distance that would have further decreased as the deposit accumulated. Alternative access points may also have existed, such as a lateral entrance proposed by Martín-Perea et al. (2025) in the narrowing zone of the Monumental sector. This hypothesis is supported by the presence of the lower marl unit in this area, which is particularly susceptible to erosion, and has been eroded since MIS7. However, apart from these inferences, no definitive evidence of a lateral entrance has been identified. Furthermore, spatial analyses have not revealed discernible patterns or voids of archaeological material in this southern area, such as those observed in the central zone, to support this hypothesis.
Grade 1

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

Crania 1

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

Crania 10

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

Crania 28

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

Crania 2

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

Crania 11

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

Crania 16

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

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

Crania 3

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

Crania 12

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

Crania 18

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

Crania 29

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

Crania 4

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

Crania 13

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

Crania 19

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

Crania 31

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

Crania 5

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

Crania 9

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

Crania 15

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

Crania 26

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

Crania 34

Grade 2

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

Crania 8

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

Crania 20

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

Crania 25

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

Crania 14

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

Crania 23

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

Crania 30

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

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

Crania 17

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

Crania 24

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

Crania 32

Grade 3

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

Crania 6

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

Crania 27

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

Crania 7

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

Crania 33

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

Crania 21

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

Crania 35

Grade 4

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

Crania 22

<figure><img src="image_49.png" /></figure>
**Fig. 10** Photographs of the primary fragments of the crania from Level 3, categorized by grades of completeness, illustrating preservation variability across specimens

Following the initial stages of the cone-shaped deposit, the rockfall process continued over time, progressively accumulating debris materials. While morphologies resulting from collapse and rockfall are typically regarded as poorly sorted (Ford and Williams 1991), no clear classification of clastic materials was identified during fieldwork. However, it was the application of geostatistics to the geological materials that revealed differential spatial distributions according to size, suggesting some degree of sorting. This size sorting was apparent in the density maps, highlighting distinct high-intensity areas for medium and small boulders (Figs. 6 and 7). Additionally, density lines along the x- and z-axes showed that peaks of higher density for small boulders correspond with slight decreases in the density of medium-sized ones (Supplementary Information 4). These patterns provide compelling evidence for sorting and the identification of the conical structure, with medium-sized boulders forming the primary sedimentary feature, while smaller boulders likely occupied the interstitial voids between the larger fragments. This distribution could reflect a static sieving effect, where finer elements percolated into available spaces rather than an active sorting mechanism, highlighting the interplay between gravitational deposition and post-depositional rearrangement within the deposit (Bacchi et al. 2014).

## Temporal implications of the conic sedimentary structure

The cone-shaped deposit evidenced by the intensity analysis of medium boulders revealed zones of lower intensity that can be interpreted as hiatuses in sedimentation, reinforcing the interpretation of a multiphsodic accumulation process (Martín-Perea et al. 2025), and pointing to a long-duration formation history clearly ruling out a single, short-term depositional event.

Dating efforts for Level 3 of Des-Cubierta Cave have been challenging, with unsuccessful radiocarbon attempts on bone remains (Baquedano et al. 2023). However, U/Th dating of speleothem S2 provides a maximum date of 135 ka, while a charcoal fragment from Level 2 gives a minimum date of around the middle of MIS3, with a calibrated age range of 43,402 to 52,964 cal BP (84.6%) and extending to 53,154 cal BP (10.9%) (Martín-Perea et al. 2025). Establishing precise timing beyond this range remains difficult. Sasowsky (2007) notes that collapse processes in caves occur on a geological scale, but measurements of rockfall rates inside karst systems remain scarce. Distinctions between short-term frost processes and larger annual cycles

have been made (Matsuoka 2008), but factors like freeze-thaw and thermal stress, relevant in external areas, appear less significant in the more stable, benign conditions inside the cave. Therefore, data from external observations (Matsuoka 2008, 2019; Bertran et al. 2015) may not be directly applicable to the specific karstic context of this study, where geomorphic processes operated more slowly, at least during phases of restricted skylight opening and relative internal stability. Examples of rockfall dating include the Greek site of Klithi, where radiocarbon dates suggest sediment accumulation over 3000 years (16,500–13,500 B.P.) (Bailey and Woodward 1997), and Shield Trap Cave, where a steady rain of pebbles formed a cone over 5,500 years, with radiocarbon dates ranging from 7,200 to 1,800 B.P. (Oliver 1989). While comparisons with Des-Cubierta Cave are cautious, these examples underscore the long-term, recurrent nature of collapse events responsible for forming such rockfall deposits.

Prior to this study, evidence for recurrent Neanderthal activity in the gallery was inferred from the vertical persistence of similar archaeological materials throughout the ~2 m thickness of Level 3—namely ungulate crania, thermoaltered remains, and Mousterian lithic tools associated with their processing (Baquedano et al. 2023). The present geostatistical analysis represents a substantial step forward, allowing the inclusion of the distinctive Level 3 assemblage within a robust formation model that integrates geological and archaeological data. This model reveals a composite process involving successive stages of clast accumulation interspersed with repeated episodes of Neanderthal use of the cave. Viewed through this integrated framework, the repeated rockfall events shaping the cone, together with the inferred sedimentation hiatuses, reinforce the interpretation of a periodic use of the gallery by Neanderthals, emphasizing the sustained and long-term character of their activities within the cave. This recurrent engagement with the confined space suggests that the introduction of crania formed part of a repeated, culturally motivated behaviour—a transmitted practice extending over an undetermined but prolonged period—reflecting the deliberate and recurrent selection of the widest part of the narrow gallery for a specific, non-subsistence purpose, the meaning of which remains uncertain.

## Implications for crania preservation

The recognition of sedimentary structures, such as the cone-shaped deposit at Des-Cubierta Cave, is critical for the dissection of palimpsests. As Bertran et al. (2019) note, identifying sedimentary structures during fieldwork is essential for taphonomic studies, but latent structures can also be revealed post-excavation using advanced analytical
<figure><img src="image_50.png" /><figcaption>**Fig. 11** Distribution of crania by grades of completeness in plan (a) and section view (b), overlaid on the medium boulder KDE map. Crania avoid the apex of the cone-shaped sedimentary structure, with better-preserved specimens (Grades 3 and 4) clustering in the periph-</figcaption></figure>

methods like geostatistics. At Des-Cubierta Cave, the morphology of the cone challenges conventional expectations of archaeological layer organization, where materials are often assumed to distribute along sub-horizontal or stratified surfaces. Instead, depositional dynamics suggest that archaeological materials may have accumulated in a pattern conforming to the cone's geometry, following its sloping and concave morphology, a factor that should be considered in future attempts to dissect the palimpsest of Unit 3. Future research on Level 3 will integrate this perspective offering a more nuanced understanding of Neanderthal behavior within this dynamic depositional context.

The growth and morphology of the cone were strongly influenced by the gallery's shape and geomorphological features. Geostatistical analyses of medium-sized boulders in the horizontal plane identified the cone's apex in the Monumental area, near the east wall, with growth extending westward, creating a cone-shaped pattern. However, the morphology observed is slightly modified by the subsequent formation of Unit H, which affected the eastern section of the sedimentary cone, adjacent to the east wall. The paleosurface beneath the cone sloped slightly southward, favoring cone growth in that direction. However, the gallery's narrowing immediately south of the Monumental sector, from 4 m to 2.7 m wide, acted as a physical barrier, limiting further expansion. This southern zone was also influenced by the S1 speleothem,

ery of the cone within the gallery's central and widest area, while more fragmented specimens occur in the narrowing zone. For details of the numbered crania see Supplementary Information 5

formed over the entire study area during MIS 7 and subsequently fragmented after gravitational collapse caused by erosion of underlying materials (Level 7), likely during a significant interglacial-glacial transition (Baquedano et al. 2023; Martín-Perea et al. 2025). While the S1 speleothem line is visible throughout the Level 3 area, its collapsed fragments were distributed differently: deeper within sandy-clay deposits in the Monumental area (Levels 5 and 6), and near their original position but verticalized and imbricated in the southern zone (Supplementary Information 1). In this southern area, some in situ remnants of S1 and its fragments created a physical barrier for low- and mid-elevation materials, resulting in a cuvette-like morphology that enabled the accumulation of clastic and mixed archaeological materials. However, no clear spatial patterns were identified in either clasts or artifacts. Above the cone shape, a high-intensity line sloping southward was detected, though its interpretation remains challenging due to significant erosional processes during the Upper Pleistocene. These processes reshaped the karstic system, dismantling the hill and the roofs of Des-Cubierta Cave. The rectilinear, southward-sloping pattern observed in the geostatistical analysis may represent a truncated morphology, visible only in part due to this extensive erosion. Alternatively, it may reflect the final collapse of the remaining roof, marking the last stage of structural instability in the cave's formation history in this section of the gallery.
The formation of the cone significantly impacted on the archaeological remains in the gallery, causing fracturing, limited transport, and some redistribution. Earlier data ruled out exogenous entry of artifacts, with taphonomic analysis revealing low degrees of rounding (14.9%) and polishing (17.7%) on bone surfaces (Baquedano et al. 2023). Similarly, lithic artifacts lacked surface alterations indicative of transport, as noted in the lithic analysis conducted by co-author B.M. Bone refits in Level 3 offered a direct means to evaluate displacement, showing a predominance of short transport distances, both in plan and section. These results primarily reflect the behavior of medium and large fragments, as the low proportion of small fragment refits introduces potential bias. Smaller fragments could have been subject to gravitational movement, sieve effects, or slightly longer transport distances.

The highest refit densities were concentrated in the narrowing zone of the gallery and at higher elevations within Level 3. These patterns suggest that the structural characteristics of the cave and the cone's growth strongly influenced material distribution. The narrow gallery, measuring between 2 and 4.5 m wide and no more than 10 m long, exhibited a paleosurface with a slight southern slope that evolved into a concave morphology over time. The apex of the cone, located near the eastern wall, grew centrifugally but was constrained northward by the slope and southward by the narrowing of the gallery walls. The in situ remains of the S1 speleothem at the southern limit (Supplementary Information 1) acted as an additional natural barrier, restricting material movement and accumulation.

While rockfall dominated as the primary sedimentary process, other formation mechanisms contributed to deviations in refit patterns. Erosive processes during the Upper Pleistocene dismantled karstic roofs and affected the upper levels of Level 3, while Holocene edaphization further disrupted materials at higher elevations. These diagenetic processes likely account for the longer transport distances observed in some instances. Refitting directionality showed a weak preference for southeast movement, consistent with the morphology of the gallery and the cone's growth, but the overall prevalence of short distances limits the interpretive value of these trends.

The spatial pattern of the crania, categorized by their degrees of completeness, highlights the impact of sedimentary and diagenetic processes on their preservation. Introduced into the gallery with anthropically modified morphology, the skulls underwent fragmentation primarily due to successive rockfalls. Field observations already indicated differences in preservation: crania in the Monumental sector were generally more complete, retaining recognizable morphology, while those in the southern zone were more fragmented, complicating their identification.

The morphology of the sedimentary cone appears to have played a hierarchical role in organizing the gallery space, influencing the distribution and preservation of the crania. Notably, no crania—regardless of their preservation state—were found at the apex of the cone. Instead, crania were distributed around the cone’s periphery, with better-preserved examples located to the north and west. These areas, distant from the cone’s primary growth area and the narrowing of the gallery, seem to have experienced fewer destructive forces. Conversely, the southern zone, characterized by the cone’s southward growth and constrained by the narrowing gallery walls and the fragmented S1 speleothem, showed poorer preservation levels. This region received a higher influx of clasts, further amplified by the slope and the presence of a southern step formed by S1’s collapse.

Vertical positioning within Level 3 also influenced cranial preservation. Initially, the roof protected the deposit from external agents, while boulders—although fragmenting materials upon impact—later, as the roof gradually lost its integrity, acted as a stabilizing matrix. Over time, the complete dismantling of the roof and subsequent erosive and edaphic processes exposed the uppermost layers to increased weathering and deterioration, contributing to longer displacement distances and reduced cranial integrity in the highest parts of the deposit.

The restoration process revealed additional insights into preservation patterns. Cranial fragments from the southern zone were smaller and frequently fire-affected (Villaescusa et al. 2022), complicating reconstruction and often preventing the restoration of original morphologies. In contrast, crania from the Monumental sector retained a high degree of completeness, often preserving the neurocranium, splanchnocranium, and cranial appendages (where applicable). These observations confirm the initial field assessment and emphasize the role of gallery morphology and cone dynamics in shaping the distribution and preservation of the remains.

Overall, the sedimentary cone functioned as both a physical and spatial organizing element within the gallery. Its growth influenced material deposition and preservation, while the loss of the roof introduced a temporal dimension to preservation patterns, exposing materials in higher levels to greater post-depositional alteration. This interaction of sedimentary and erosive processes underscores the complexity of Level 3’s formation and its implications for archaeological material integrity.

## Behavioural implications

As emphasized at the beginning of this paper, the importance of understanding the formative histories of archaeological deposits—particularly in contexts where complex
post-depositional dynamics may obscure behavioural patterns—is essential. In Level 3 of Des-Cubierta Cave, this approach is especially crucial, given the exceptional assemblage of at least 35 individuals of large ungulates with defensive cranial appendages showing a standardized anatomical pattern and clear traces of anthropogenic modification (cut and percussion marks) on several well-preserved specimens.

Across Europe, Middle Palaeolithic sites display a remarkable behavioural diversity, from specialized hunting and small-game exploitation to pigment use, bird-claw collection, the construction of possible ritual structures such as Bruniquel, and funerary practices (Romagnoli et al. 2022). Within this broad spectrum, however, the deliberate treatment or accumulation of animal crania remains rare and unevenly documented, often appearing in unusual archaeological contexts such as stone alignments, confined niches, or associations with human remains. At the Grotte du Prince (Italy), de Villeneuve (1906) described in the Foyer B level an accumulation that included several ibex skulls, four *massacres* (frontal portions of bovine crania with attached horns), and exceptionally preserved *Cervus elaphus* antlers, as well as a juvenile elephant skull and a rhinoceros mandible intentionally placed on a natural ledge. At Lazaret (France), de Lumley et al. (2018) reported two caprine crania with attached horns incorporated into a block alignment that delimited a small lateral recess of the cave. In Italy, *Bison priscus* and cervid crania were recorded in Grotta Grande (layers 6 and 8) in association with stone alignments and stalagmitic surfaces (Ronchitelli et al. 2011), and at Molare Rock Shelter, where a *Bison priscus* skull was found in a niche beneath the rock wall (layer 43) and seven additional cranial portions with horns and antlers were discovered near a semicircular stone alignment (layer 56) (Peretto et al. 2004; Spagnolo et al. 2020). In a different but related sphere, cranial remains of ungulates and carnivores have also been associated with Neanderthal burials at Amud (Israel) (Rak et al. 1994), Teshik-Tash (Uzbekistan) (Movius 1953a, b), and Regourdou (France) (Bonifay and Vandermeersch 1962), not without controversy (Gargett 1989; Pelletier et al. 2017).

Within this scattered and heterogeneous body of evidence, the case of Des-Cubierta Cave stands out for the size and recurrence of the behaviour, as well as for the absence of evidence of Neanderthal dwelling in the gallery. The number of crania recovered—at least 35—and their repeated deposition within the same confined gallery over an extended period underscore the sustained and reiterated nature of this practice. The results of the present geostatistical analysis demonstrate that these depositions were not the product of a single event but of successive episodes embedded within a long-term process of gallery use, offering unique evidence

of a structured, transmitted, and culturally motivated practice centred on horned herbivore crania.

## Conclusions

The spatial and geostatistical analysis of geological and archaeological materials, combined with the study of the bone refits and crania preservation, has provided critical insights into the natural and anthropogenic processes that shaped Level 3 of Des-Cubierta Cave. Specifically, this study has revealed:

1. Geological and archaeological materials exhibit distinct spatial patterns, confirming separate formation dynamics for each material type.

2. The absence of archaeological materials within clast-rich areas indicates an initial phase of debris cone development from rockfall processes. This phase was followed by Neanderthal activity, as evidenced by the overlying archaeological materials.

3. The sedimentary cone-shaped morphology, determined through size-based geostatistical analysis of clasts, influenced the spatial organization of the gallery and the paleosurface available during Neanderthal use.

4. The identification of sedimentary hiatuses within the cone structure reveals distinct phases of growth, contributing to a refined understanding of the geological and anthropogenic temporality of Level 3.

5. Bone refits, especially medium and large fragments, show limited movement influenced by the gallery's morphology and sedimentary processes, including rockfalls, erosion, and edaphic alterations.

6. Crania completeness varies with proximity to the sedimentary cone and altitudinal position. Crania near the apical area on the horizontal plane exhibit better integrity, while those in southern zones and higher elevations are more fragmented due to gravitational, erosive and edaphic processes.

7. Two distinct zones of material preservation were identified: the central zone, characterized by the widest part of the gallery and the sedimentary cone, where cranial remains are better preserved and bone movement is minimal; and the southern area, defined by the gallery narrowing and imbricated speleothem fragments, creating an irregular surface with poorer preservation of cranial remains.

This study demonstrates the utility of geostatistical methods in complementing traditional archaeological approaches to investigate site formation processes. The sedimentary dynamics identified in Des-Cubierta Cave are key to
understanding the distribution and preservation of archaeo-
logical materials. These findings will be pivotal in future
studies, particularly those focused on disentangling the
temporality of the archaeological record and dissecting
the palimpsest, processes inherently shaped by the deposi-
tional events and sedimentary structures documented in this
research.

The integration of geospatial with geoarchaeological pers-
spectives enables a more comprehensive interpretation of
spatial and temporal patterns in karst environments, provid-
ing a comparative basis for evaluating sites with analogous
sedimentary histories.

Beyond its methodological and geoarchaeological con-
tributions, this study provides a framework for interpreting
the recurrent use of Des-Cubierta cave by Neanderthals.
The integration of geological, spatial, and taphonomic data
demonstrates that the accumulation of large herbivore cra-
nia was not a single depositional event, but rather the result
of repeated episodes embedded within a long-term process
of gallery use. This sustained and reiterated behaviour high-
lights the structured and transmitted nature of this practice,
adding a significant piece to the broader discussion on the
complexity and symbolic potential of Neanderthal cultural
expressions.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s12520-025-02382-5.

**Acknowledgements** The authors are grateful to the Pinilla del Valle team, including all students and volunteers who participated in the fieldwork. Particular thanks are due to the restoration team, under the direction of María Cruz Ortega, and comprising Aroa Serrano, Laura Gómez Tejerina, Bárbara Martín, Sofía de León, Andrea Díaz-Cortés, Lucía Hernández-Vivanco, Elena Cubedo Izquierdo, Laura Gómez Morgado, Alejandro Velázquez, José Ángel Correa, Javier Casado, as well as all other restorers who contributed during the excavation campaigns. We thank Javier Trueba (Madrid Scientific Films), Mario Torquemada (MAPCM), and Alfonso Dávila for their photographic work, and the latter also for providing the site's photogrammetric data. We are also grateful to José Antonio Vallejo for his constant support during the excavations. L.V. and S.G-S. were beneficiaries of predoc-toral contracts under the programme H2019/HUM-5840, provided by the Dirección General de Investigación e Innovación Tecnológica de la Comunidad de Madrid (co-funded by the European Social Fund). L.V. also received support from the 2023 Research Grants of the Fundación Banco Sabadell. M.G.-P. and L.V. are supported by a research contract within the PHS-2024/PH-HUM-469 programme of the Comunidad de Madrid. LC-S is funded by the European Union's Horizon Europe pro-gramme under the Marie Skłodowska-Curie Postdoctoral Fellowship grant agreement No. 101152206. A.M is funded by a post-doctoral research grant by the Fyssen Foundation. This article is a contribution to the VALLENEANDERTAL-CM project (PHS-2024/PH-HUM-469), which is part of the R&D projects call among research groups focused on Human and Social Processes, funded by the Community of Madrid. The Community of Madrid funds the excavations at the Pinilla del Valle sites through its Archaeological and Paleontological Museum.

The Grupo Mahou-San Miguel and Canal Isabel II are thanked for their sponsorship, as are the Dirección General de Juventud and the staff of the Los BatanesYouth Hostel for making it possible for the team to stay at their facilities during the excavation campaigns. We also thank the Sierra de Guadarrama National Park and the Pinilla del Valle City Council for their collaboration, along with the many people who, over the years, have been part of the excavation team. We thank the Dirección General de Patrimonio Cultural for authorizing the permission to carry out the excavations.

**Authors' contributions** L.V.: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review and editing. E.B.: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing – review & editing. D.M-P.: Conceptualization, Data curation, Investigation, Resources, Visualization (Fig.1 y 2), Writing – review and editing. B.M.: Data curation, Investigation, Resources (lithic materials data), Writing – review and editing. M.A.G-P.: Conceptualization, Data curation, Investigation, Resources (cranial refitting and crania completeness data), Visualization (Fig.10), Writing – review and editing. L.C-S.: Conceptualization, Formal analysis, Methodology, Supervision, Writing – review and editing. A.I.O.: Data curation, Investigation, Resources, Writing – review and editing. R.H.: Conceptualization, Data curation, Investigation, Resources, Writing – review and editing. C.L.: Data curation, Funding acquisition, Investigation, Resources, Writing – review & editing. M.C.O.: Data curation, Investigation, Resources, Writing – review and editing. S.G-S.: Conceptualization, Data curation, Investigation, Resources, Visualization (Fig.1, Fig.11), Writing – review and editing. A.M.: Data curation, Investigation, Writing – review and editing. N.G.: Data curation, Investigation, Writing – review and editing. D.J.A-L.: Data curation, Investigation, Writing – review and editing. R.G-G. Data curation, Investigation, Writing – review and editing. L.R.: Data curation, Investigation, Writing – review and editing. A.P-G.: Investigation, Project administration, Writing – review and editing. J.L.A.: Investigation, Project administration, Writing – review and editing. This research was carried out within the framework of a research project co-directed by E.B., A.P-G., and J.L.A. All authors approved the final version of the manuscript and agree to be held accountable for the content therein.

**Data availability** Data, code, and visual outputs are openly available in our research compendium (https://osf.io/frzts/ DOI https://doi.org/10.17605/OSF.IO/FRZTS), following the reproducible research framework outlined by Marwick (2017).

Declarations

Competing interests The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
## References

Abrunhosa A (2020) Proveniência das matérias-primas e a sua relação com a tipologia: estudo do conjunto lítico Moustierense dos sítios arqueológicos do Calvero de la Higuera (Madrid, Espanha). Dis- sertación, Universidade do Algarve. http://hdl.handle.net/10400. 1/16724

Abrunhosa A, Pereira T, Márquez B et al (2019) Understanding nean- derthal technological adaptation at Navalmaíllo rock shelter (Spain) by measuring lithic Raw materials performance variabil- ity. Archaeol Anthropol Sci 11:5949–5962. https://doi.org/10.100 7/s12520-019-00826-3

Agostinelli C, Lund U (2023) R package 'circular': Circular Statistics (version 0.5-1) https://doi.org/10.32614/CRAN/package.circular

Alférez F, Molero G, Maldonado E et al (1982) Descubrimiento Del primer Yacimiento Cuaternario (Riss-Würm) de vertebrados Con Restos Humanos En La provincia de Madrid (Pinilla Del Valle). Coloquios De Paleontología 37:15–32

Alperson-Afil N (2017) Spatial analysis of fire: archaeological approach to recognizing early fire. Curr Anthropol 58:S258– S266. https://doi.org/10.1086/692721

Álvarez-Fernández A, Márquez B, García-González R, et al (2025) Experimental use-wear on non-flint tools: Unraveling butcher- ing processes at the Cueva Des-Cubierta site (Pinilla del Valle, Spain). J Archaeol Sci Rep 62:105056. https://doi.org/10.1016/j .jasrep.2025.105056

Álvarez-Lao D, Arsuaga JL, Baquedano E, Pérez-González A (2013) Last interglacial (MIS 5) ungulate assemblage from the central Iberian peninsula: the camino cave (Pinilla Del Valle, Madrid, Spain). Palaeogeogr Palaeoclimatol Palaeoecol 374:327–337. h tps://doi.org/10.1016/j.palaeo.2013.01.025

Arilla M, Rosell J, Blasco R (2020) A neo-taphonomic approach to human campsites modified by carnivores. Sci Rep 10:1–15. https://doi.org/10.1038/s41598-020-63431-8

Arriaza MC, Huguet R, Laplana C et al (2017) Lagomorph predation represented in a middle palaeolithic level of the Navalmaíllo rock shelter site (Pinilla Del Valle, Spain), as inferred via a new use of classical taphonomic criteria. Quatern Int 436:294–306. https://d oi.org/10.1016/j.quaint.2015.03.040

Arsuaga JL, Baquedano E, Pérez-González A et al (2010) El Yacimiento arqueopaleontológico Del Pleistoceno superior de La Cueva Del camino En El Calvero de La Higuera (Pinilla Del Valle, Madrid). Actas de La 1ª Reunión de científicos sobre cubiles de Hiena (y Otros grandes carnívoros) En Los Yacimientos arqueológicos de La Península Ibérica. Alcalá de Henares, pp 421–442

Arsuaga JL, Baquedano E, Pérez-González A et al (2012) Understand- ing the ancient habitats of the last-interglacial (late MIS 5) nean- derthals of central iberia: paleoenvironmental and taphonomic evidence from the Cueva Del camino (Spain) site. Quatern Int 275:55–75. https://doi.org/10.1016/j.quaint.2012.04.019

Arteaga-Brieba A (2024) Among stones and bison. The lithic assem- blage of Gran Dolina TD10.2 (Atapuerca). Technological and spatial implications of a specialised Middle Pleistocene kill- butchering site. Dissertation, Universitat Rovira i Virgil. https://h dl.handle.net/20.500.11797/TDX4377

Arteaga-Brieba A, Courtenay LA, Cobo-Sánchez L et al (2023) An archaeostratigraphic consideration of the Gran Dolina TD10.2 cultural sequence from a quantitative approach. Q Sci Rev 309:108033. https://doi.org/10.1016/j.quascirev.2023.108033

Bacchi V, Recking A, Eckert N et al (2014) The effects of kinetic sort- ing on sediment mobility on steep slopes. Earth Surf Proc Land 39:1075–1086. https://doi.org/10.1002/esp.3564

Baddeley A, Rubak E, Turner R (2016) Spatial point Patterns. Meth- odology and applications with R. CRS Press Taylor & Francis Group, New York

Baddeley A, Turner R (2005) Statpat: An R package for analyzing spatial point pattens. J Stat Softw 12(6):1–42 https://doi.org/10. 18637/jss.v012.i06

Bailey G, Woodward J (1997) The klithi deposits: Sedimentology, stratigraphy and chronology. In: Bailey GN (ed) Klithi: palaeo- lithic settlement and quaternary landscapes in Northwest Greece. McDonald Institute, Cambridge, pp 61–94

Baquedano E, Arsuaga JL, Pérez-González A et al (2016a) The Des- Cubierta Cave (Pinilla del Valle, Comunidad de Madrid, Spain): a Neanderthal site with a likely funerary/ritualistic connection. In: European Society for the study of Human Evolution, 6th Annual Meeting. Madrid (Spain), p 41 https://www.eshe.eu/static/eshe/p eshe/PESHE_2016.pdf

Baquedano E, Arsuaga JL, Pérez-González A et al (2023) A symbolic neanderthal accumulation of large herbivore Crania. Nat Hum Behav 7:342–352. https://doi.org/10.1038/s41562-022-01503-7

Baquedano E, Laplana C, Arsuaga JL et al (2016b) Selection of cave shelter by neanderthals (Homo neanderthalensis) and spotted hyaenas (Crocuta crocuta) at the Calvero de la Higuera (Pinilla del Valle, Madrid, Spain). Arpi Arqueologia y Prehistoria del Interior peninsular 04 Extra:5–19 https://cir.cenieh.es/handle/ 20.500.12136/1301?mode=full

Baquedano E, Márquez B, Pérez-González A et al (2012) Neander- tales en el valle del Lozoya: los yacimientos paleolíticos del Cal-vero de la Higuera (Pinilla del Valle, Madrid). Mainake XXXIII: 83–100 https://hdl.handle.net/20.500.14352/44213

Bargalló A, Gabucio MJ, Rivals F (2016) Puzzling out a palimpsest: testing an interdisciplinary study in level O of abric Romaní. Qua- tern Int 417:51–65. https://doi.org/10.1016/j.quaint.2015.09.066

Bel MÁ (2022) Lithic taphonomy and the use of the Spatial infor- mation: assessing palaeolithic sequence of Cova de les cendres (Teulada-Moraira, Spain). Archaeol Anthropol Sci 14:173. https:// doi.org/10.1007/s12520-022-01642-y

Benito-Calvo A, de la Torre I (2011) Analysis of orientation patterns in Olduvai bed I assemblages using GIS techniques: implications for site formation processes. J Hum Evol 61:50–60. https://doi.or g/10.1016/j.jhevol.2011.02.011

Bertran P, Beauval C, Boulogne S et al (2015) Experimental archaeol- ogy in a mid-latitude periglacial context: insight into site forma- tion and taphonomic processes. J Archaeol Sci 57:283–301. https ://doi.org/10.1016/j.jas.2015.02.039

Bertran P, Hétu U B, Texier J-P, Van Steijn H (1997) Fabric character- istics of subaerial slope deposits. Sedimentology 44:1–16. https:// doi.org/10.1111/j.1365-3091.1997.tb00421.x

Bertran P, Klaric L, Lenoble A et al (2010) The impact of perigl- acial processes on palaeolithic sites: the case of sorted patterned grounds. Quatern Int 214:17–29. https://doi.org/10.1016/j.quain t.2009.10.021

Bertran P, Lenoble A, Todisco D et al (2012) Particle size distribu- tion of lithic assemblages and taphonomy of palaeolithic sites. J Archaeol Sci 39:3148–3166. https://doi.org/10.1016/j.jas.2012 .04.055

Bertran P, Todisco D, Bordes J et al (2019) Perturbation assessment in archaeological sites as part of the taphonomic study: a review of methods used to document the impact of natural processes on site formation and archaeological interpretations. Paléo 30–1:52–75. https://doi.org/10.4000/paleo.4378

Bevan A (2012) Spatial methods for analysing large-scale artefact inventories. Antiquity 86:492–506. https://doi.org/10.1017/S000 3598X0006289X

Bevan A (2020) Spatial Point Patterns and Processes. In: Mark Gill- ings, Hacigüzeller P, Lock G (eds) Archaeological spatial analy- sis. a methodological guide. Routledge, Oxon and New York, pp 60–76

Binford LR, Mills MGL, Stone NM (1988) Hyena scavenging behavior and its implications for the interpretation of faunal
assemblages from FLK 22 (the Zinj floor) at Olduvai gorge. J
Anthropol Archaeol 7:99-135. https://doi.org/10.1016/0278-416
5(88)90011-6

Blain H-A, Laplana C, Sánchez-Bandera C et al (2022) A warm and humid paleoecological context for the neanderthal mountain settlement at the Navalmaillo rockshelter (Iberian central System, Madrid). Q Sci Rev 293:107727. https://doi.org/10.1016/j.quascirev.2022.107727

Blain H-A, Laplana C, Sánchez-Bandera C et al (2025) New paleo-
ecological insights for the late pleistocene neanderthal mountain
occurrence at Buena Pinta cave (Iberian central System, Pinilla-
del-Valle, Madrid, Spain). Q Sci Rev 362:109355. https://doi.org
/10.1016/j.quascirev.2025.109355

Blott SJ, Pye K (2012) Particle size scales and classification of sedi-
ment types based on particle size distributions: review and rec-
ommended procedures. Sedimentology 59(7):2071–2096. https://
doi.org/10.1111/j.1365-3091.2012.01335.x

Bonifay E, Vandermeersch B (1962) Dépôts rituels d'ossements d'ours
Dans Le Gisement moustérien du regourdou (Montignac, Dordo-
gne). C R Hebd Acad Sci 255:1635-1636

Canals A, Vallverdú J, Carbonell E (2003) New archaeo-stratigraphic data for the TD6 level in relation to homo antecessor (Lower Pleistocene) at the site of Atapuerca, north-central Spain. Geochrono-chaeology 18:481–504. https://doi.org/10.1002/gea.10071

Cobo-Sánchez L (2020) Taphonomic and spatial study of the archaeo-
logical site DS from Bed I in Olduvai Gorge (Tanzania). Disserta-
tion, Universidad Complutense de Madrid. https://hdl.handle.net
/20.500.14352/11463

Cobo-Sánchez L, Aramendi J, Domínguez-Rodrigo M (2014) Orientation patterns of wildebeest bones on the lake Masek floodplain (Serengeti, Tanzania) and their relevance to interpret anisotropy in the Olduvai lacustrine floodplain. Quaternary International 322-323:277-284. https://doi.org/10.1016/j.quaint.2013.07.130

Cziesla E (1990) On Refitting of Stone Artifacts. In: Cziesla E, Eickhoff S, Arts N, Winter D (eds) The big puzzle: international symposium on refitting stone artefacts, monrepos, 1987. Holos, Bonn, pp 9–44

de la Peña P, Thomas M, Molefyane TR (2022) Particle size distri-
bution: an experimental study using Southern African reduction
methods and Raw materials. PLoS ONE 17:e0278867. https://doi
.org/10.1371/journal.pone.0278867

de la Torre I, Benito-Calvo A (2013) Application of GIS methods to
retrieve orientation patterns from imagery; a case study from beds
I and II, Olduvai gorge (Tanzania). J Archaeol Sci 40:2446–2457.
https://doi.org/10.1016/j.jas.2013.01.004

de la Torre I, Benito-Calvo A, Proffitt T (2018) The impact of hydraulic
processes in Olduvai beds I and II, Tanzania, through a particle
dimension analysis of stone tool assemblages. Gearchaeology
33:218–236. https://doi.org/10.1002/gea.21629

de Lumley H, Barrier P, Cauche D et al (2018) Comportements de
subsistance et mode de Vie des homo erectus. In: De Lumley
M-A (ed) Les restes Humains fossiles de La Grotte du Lazaret.
Nice, Alpes-Maritimes. France. CNRS Editions, Paris, Paris, pp
155-215

Deschamps M, Zilhão J (2018) Assessing site formation and assem-
blage integrity through stone tool refitting at Gruta Da Oliveira
(Almonda karst system, Torres Novas, Portugal): A middle paleo-
lithic case study. PLoS ONE 13:e0192423. https://doi.org/10.137
1/journal.pone.0192423

de Villeneuve L (1906) Les Grottes de Grimaldi (Baoussé-Roussé). Imprimerie de Monaco, Monaco

Dibble HL, Aldeias V, Goldberg P et al (2015) A critical look at evi-
dence from La Chapelle-aux-Saints supporting an intentional
neandertal burial. J Archaeol Sci 53:649–657. https://doi.org/10
.1016/j.jas.2014.04.019

Diez-Martín F, Cobo-Sánchez L, Baddeley A et al (2021) Tracing the
Spatial imprint of Oldowan technological behaviors: A view from
DS (Bed I, Olduvai Gorge, Tanzania). PLoS ONE 16:e0254603.
https://doi.org/10.1371/journal.pone.0254603

Domínguez-Rodrigo M, Barba Egido R, Egeland CP (2007) Decon-
structing olduvai: a taphonomic study of the bed i sites. vertebrate
paleobiology and paleoanthropology series. Springer Nether-
lands, Dordrecht https://doi.org/10.1007/978-1-4020-6152-3

Domínguez-Rodrigo M, Bunn HT, Pickering TR et al (2012) Autoch-
thony and orientation patterns in Olduvai bed I: a re-examination
of the status of post-depositional biasing of archaeological assem-
blages from FLK North (FLKN). J Archaeol Sci 39:2116–2127. h
ttps://doi.org/10.1016/j.jas.2012.02.027

Domínguez-Rodrigo M, Cobo-Sánchez L (2017) The Spatial pattern-
ing of the social organization of modern foraging homo sapiens:
A methodological approach for Understanding social organiza-
tion in prehistoric foragers. Palaeogeogr Palaeoclimatol Palaeo-
ecol 488:113–125. https://doi.org/10.1016/j.palaeo.2017.06.008

Domínguez-Rodrigo M, Uribelarrea D, Santonja M et al (2014) Autochthonous anisotropy of archaeological materials by the action of water: experimental and archaeological reassessment of the orientation patterns at the Olduvai sites. J Archaeol Sci 41:44–68. https://doi.org/10.1016/j.jas.2013.07.025

Falcucci A, Giusti D, Zangrossi F et al (2024) Refitting the context: A reconsideration of cultural change among early homo sapiens at fumane cave through blade break Connections, Spatial Taphonomy, and lithic technology. J Paleolithic Archaeol 8:2. https://doi.org/10.1007/s41982-024-00203-0

Fisher JW (1995) Bone surface modifications in zooarchaeology. J
Archaeol Method Theory 2:7-68. https://doi.org/10.1007/BF02
228434

Ford D, Williams P (1991) Karst geomorphology and hydrology. Geographical J 157:87. https://doi.org/10.2307/635167

Fraile-Márquez C, Diez-Martín F, Duque-Martinez J et al (2022) Fac-
ing the palimpsest conundrum: an archaeo-stratigraphic approach
to the intra-site analysis of SHK extension (Bed II, Olduvai
Gorge, Tanzania). Archaeol Anthropol Sci 14:230. https://doi.or
g/10.1007/s12520-022-01691-3

Galindo-Pellicena MA, Arsuaga JL, Laplana C et al (2019) Distinguishing between Bos and Bison petrous bones. A case study: bovines from the Des-Cubierta cave (Pinilla Del Valle, Madrid). Span J Palaeont 34(2):257–268. https://doi.org/10.7203/sjp.34.2.16115

García-Moreno A, Hutson JM, Villaluenga A et al (2023) Connect-
ing bones at Schöningen 13II-4 spear horizon: an analysis of site
formation and human activity through faunal refitting. Archaeol
Anthropol Sci 15:178. https://doi.org/10.1007/s12520-023-0187
6-4

García-Moreno A, Smith GM, Kindler L et al (2016) Evaluating the incidence of hydrological processes during site formation through orientation analysis. A case study of the middle palaeolithic Lake-land site of Neumark-Nord 2 (Germany). J Archaeol Science: Rep 6:82–93. https://doi.org/10.1016/j.jasrep.2016.01.023

Gargett RH (1989) Grave shortcomings: the evidence for neandertal
burial. Curr Anthropol 30:157-190

Giusti D, Arzarello M (2016) The need for a taphonomic perspective
in Spatial analysis: formation processes at the early pleistocene
site of Pirro Nord (P13), Apricena, Italy. J Archaeol Science: Rep
8:235–249. https://doi.org/10.1016/j.jasrep.2016.06.014

Goder-Goldberger M, Gilead I, Paixão E et al (2025) New insights on an old excavation: Re-visiting the late middle palaeolithic site of far’ah II, North-western Negev, Israel. J Paleolithic Archaeol 8:13. https://doi.org/10.1007/s41982-025-00212-7

Goldberg P, Macphail RI (2005) Practical and theoretical geoarchaeol-
ogy. Wiley, Oxford
Hazelton ML, Davies TM (2009) Inference based on kernel estimates of the relative risk function in geographical epidemiology. *Biom* J 51:98–109. https://doi.org/10.1002/bimj.200810495

Hofman JL (1986) Vertical movement of artifacts in alluvial and stratified deposits. Curr Anthropol 27:163–171. https://doi.org/10.1086/203414

Isaac GL (1983) Bones in contention: competing explanations for the juxtaposition of early pleistocene artifacts and faunal remains. In: Clutton-Brock J, Grigson C (eds) Animals and archaeology: 1. Hunters and their prey. B.A.R., Oxford

Laplana C, Sevilla P, Arsuaga JL et al (2015) How Far into Europe did Pikas (Lagomorpha: Ochotonidae) go during the pleistocene? New evidence from central Iberia. PLoS ONE 10:e0140513. http://doi.org/10.1371/journal.pone.0140513

Laplana C, Sevilla P, Blain H-A et al (2016) Cold-climate rodent indicators for the late pleistocene of central Iberia: new data from the Buena Pinta cave (Pinilla Del Valle, Madrid Region, Spain). Comptes Rendus - Palevol 15:696–706. https://doi.org/10.1016/j.crpv.2015.05.010

Lee Lyman R (1985) Bone frequencies: differential transport, in situ destruction, and the MGUI. J Archaeol Sci 12:221–236. https://doi.org/10.1016/0305-4403(85)90022-6

Lenoble A, Bertran P (2004) Fabric of palaeolithic levels: methods and implications for site formation processes. J Archaeol Sci 31:457–469. https://doi.org/10.1016/j.jas.2003.09.013

Lloveras L, Moreno-García M, Nadal J, Zilhão J (2011) Who brought in the rabbits? Taphonomical analysis of mousterian and solutean leporid accumulations from Gruta do Caldeirão (Tomar, Portugal). J Archaeol Sci 38:2434–2449. https://doi.org/10.1016/j.jas.2011.05.012

Lund U, Agostinelli C (2025) CircStats: Circular Statistics, from Topics in Circular Statistics (2001) (version 0.2-7) https://doi.org/10.32614/CRAN/package.CircStats

Luzón C, Yravedra J, Courtenay LA et al (2021) Taphonomic and Spatial analyses from the early pleistocene site of Venta micena 4 (Orce, Guadix-Baza Basin, Southern Spain). Sci Rep 11:1–18. https://doi.org/10.1038/s41598-021-93261-1

Martín-Perea DM (2021) Tafonomía de los yacimientos paleontológicos de Batallones-3 y Batallones-10 (Vallesiense Superior, Madrid, España). Dissertation, Universidad Complutense de Madrid. https://hdl.handle.net/20.500.14352/3641

Martín-Perea DM, Baquedano E, Arsuaga JL et al (2025) Stratigraphy and geomorphology of Des-Cubierta cave (Pinilla Del Valle, Madrid, Spain): geological insights into a neanderthal symbolic accumulation of large Crania. J Quat Sci 1–19. https://doi.org/10.1002/jqs.3722

Martín-Perea DM, Courtenay LA, Domingo MS, Morales J (2020) Application of artificially intelligent systems for the identification of discrete fossiliferous levels. PeerJ 2020:1–25. https://doi.org/10.7717/peerj.8767

Marwick B (2017) Computational reproducibility in archaeological research: basic principles and a case study of their implementation. J Archaeol Method Theory 24:424–450. https://doi.org/10.1007/s10816-015-9272-9

Matsuoka N (2008) Frost weathering and Rockwall erosion in the southeastern Swiss alps: Long-term (1994–2006) observations. Geomorphology 99:353–368. https://doi.org/10.1016/j.geomorph.2007.11.013

Matsuoka N (2019) A multi-method monitoring of timing, magnitude and origin of rockfall activity in the Japanese alps. Geomorphology 336:65–76. https://doi.org/10.1016/j geomorph.2019.03.023

McPherron SJP (2005) Artifact orientations and site formation processes from total station proveniences. J Archaeol Sci 32:1003–1014. https://doi.org/10.1016/j.jas.2005.01.015

McPherron SP (2018) Additional statistical and graphical methods for analyzing site formation processes using artifact orientations.

PLoS ONE 13(1):e0190195. https://doi.org/https://doi.org/ https://doi.org/10.1371/journal.pone.0190195

Mendez-Quintas E, Panera J, Altamura F et al (2019) Gombore II (Melka Kunture, Ethiopia): A new approach to formation processes and Spatial patterns of an early pleistocene Acheulean site. J Archaeol Sci 108:104975. https://doi.org/10.1016/j.jas.2019.104975

Merino-Pelaz A, Cobo-Sánchez L, Organista E et al (2024) Unraveling the Spatial imprint of hominin and carnivore accumulations in early pleistocene African sites. Archaeol Anthropol Sci 16:128. https://doi.org/10.1007/s12520-024-02020-6

Mielgo C, Huguet R, Laplana C et al (2024) Intra-site Spatial approaches based on taphonomic analyses to characterize assemblage formation at pleistocene sites: a case study from Buena Pinta cave (Pinilla Del Valle, Madrid, Spain). Archaeol Anthropol Sci 16:1–28. https://doi.org/10.1007/s12520-023-01913-2

Méndez-Quintas E, Santonja M, Pérez-González A et al (2022) Exploring the formation processes on open-air palaeolithic sites: A late middle pleistocene Acheulean assemblage at Arbo site (Miño river basin, Spain). J Archaeol Science: Rep 43:103453. https://doi.org/10.1016/j.jasrep.2022.103453

Moclán A, Cobo-Sánchez L, Domínguez-Rodrigo M et al (2023a) Spatial analysis of an early middle palaeolithic kill/butchering site: the case of the Cuesta de La Bajada (Teruel, Spain). Archaeol Anthropol Sci 15:91. https://doi.org/10.1007/s12520-023-01792-7

Moclán A, Domínguez-Rodrigo M, Huguet R et al (2024) Deep learning identification of anthropogenic modifications on a carnivore remain suggests use of hyena pelts by neanderthals in the Navalmaíllo rock shelter (Pinilla Del Valle, Spain). Q Sci Rev 329:108560. https://doi.org/10.1016/j.quascirev.2024.108560

Moclán A, Huguet R, Márquez B et al (2020) Identifying the bonebreaker at the Navalmaíllo rock shelter (Pinilla Del Valle, Madrid) using machine learning algorithms. Archaeol Anthropol Sci 12:46. https://doi.org/10.1007/s12520-020-01017-1

Moclán A, Huguet R, Márquez B et al (2021) A neanderthal hunting camp in the central system of the Iberian peninsula: A Zooarchaeological and taphonomic analysis of the Navalmaíllo rock shelter (Pinilla Del Valle, Spain). Q Sci Rev 269:107142. https://doi.org/10.1016/j.quascirev.2021.107142

Moclán A, Huguet R, Márquez B et al (2023b) Identifying activity areas in a neanderthal hunting camp (the Navalmaíllo rock shelter, Spain) via Spatial analysis. Archaeol Anthropol Sci 15:44. https://doi.org/10.1007/s12520-023-01746-z

Morin E, Tsanova T, Sirakov N et al (2005) Bone refits in stratified deposits: testing the chronological grain at Saint-Césaire. J Archaeol Sci 32:1083–1098. https://doi.org/10.1016/j.jas.2005.02.009

Morrow TM (1996) Lithic refitting and archaeological site formation processes. In: Odell GH (ed) Stone Tools. Theoretical insights into human prehistory. Springer Science+Business Media, New York, pp 345–373

Movius HL (1953a) The mousterian cave of Teshik-Tash, southeastern Uzbekistan, central Asia. Bull Am School Prehistoric Res 17:11–71

Movius HL (1953b) Palaeolithic and Mesolithic Sites in Soviet Central Asia. Proceedings of the American Philosophical Society 97:383–421

Márquez B, Mosquera M, Pérez-González A et al (2013) Evidence of a Neanderthal-Made Quartz-Based technology at Navalmaíllo rockshelter (Pinilla Del Valle, Madrid Region, Spain). J Anthropol Res 69:373–395. https://doi.org/10.3998/jar.0521004.0069.306

Oliver JS (1989) Analogues and site context: bone damages from shield trap cave (24CB91), carbon County, Montana, USA. In: Bonnichsen R, Sorg MH (eds) Bone modification. Center for the Study of the First Amerioane, Maine, pp 73–98
Panera J, Rubio-jara S, Domínguez-Rodrigo M et al (2019) Assessing functionality during the early Acheulean in level TKSF at Thiongo Korongo site (Olduvai Gorge, Tanzania). Quatern Int 526:77-98. https://doi.org/10.1016/j.quaint.2019.09.013

Pelletier M, Royer A, Holliday TW et al (2017) Rabbits in the grave! Consequences of bioturbation on the neandertal burial at regour-dou (Montignac-sur-Vézère, Dordogne). J Hum Evol 110:1-17. https://doi.org/10.1016/j.jhevol.2017.04.001

Peretto C, Biagi P, Boschian G et al (2004) Living-floors and structures from the lower paleolithic to the bronze age in Italy. Coll Antropol 28:63-88

Peters C, van Kolfschoten T (2020) The site formation history of Schöningen 13II-4 (Germany): testing different models of site formation by means of Spatial analysis, Spatial statistics and orientation analysis. J Archaeol Sci 114:105067. https://doi.org/10.1016/j.jas.2019.105067

Petraglia MD, Akoshima K, Straus LG (1994) Interpreting the formation at the Abri dufaure: an upper paleolithic site in Southwestern France. J Anthropol Archaeol 13:139-151. https://doi.org/10.1009/6/jaar.1994.1009

Petraglia MD, Potts R (1994) Water flow and the formation of early pleistocene artifact sites in Olduvai Gorge, Tanzania. J Anthropol Archaeol 13:228-254. https://doi.org/10.1006/jaar.1994.1014

Pollarolo L, Wilkins J, Kuman K, Galletti L (2010) Site formation at Kudu koppie: A late earlier and middle stone age site in Northern Limpopo Province, South Africa. Quatern Int 216:151-161. https://doi.org/10.1016/j.quaint.2009.11.019

Pérez-González A, Baquedano E, Arsuaga JL et al (2010) Aproximación geomorfológica a los yacimientos del Pleistoceno Superior del Calvero de la Higuera en el Valle Alto del lozoya (Sistema Central español, Madrid). In: Actas de la 1ª Reunión de científicos sobre cubiles de hiena (y otros grandes carnívoros) en los yacimientos arqueológicos de la Península Ibérica. pp 403-420

Rabuñal Gayo JR (2021) Análisis de la distribución espacial y estudio morfo-técnico de la industria lítica de los yacimientos mesolíticos de El Arenal de la Virgen y Casa Corona (Villena, Alicante). Dissertation, Universidad de Alicante http://hdl.handle.net/10045/142527

Rabuñal JR, Gómez-Puche M, Polo-Díaz A, Fernández-López de Pablo J (2023) Unraveling early holocene occupation patterns at El arenal de La Virgen (Alicante, Spain) open-air site: an integrated palimpsest analysis. Archaeol Anthropol Sci 15:108. https://doi.org/10.1007/s12520-023-01805-5

Rak Y, Kimbel WH, Hovers E (1994) A neandertal infant from Amud Cave, Israel. J Hum Evol 26:313-324. https://doi.org/10.1006/jhrev.1994.1019

Rapp GJ, Hill CL (2006) Geoarchaeology: the Earth-Science approach to archaeological interpretation, 2nd edn. Yale University Press. New Haven and London

R Core Team (2024) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Ripley BD (1979) Tests of “Randomness” for Spatial Point Patterns. J R Stat Soc Se B: Stat Methodol 41:368-374. https://doi.org/10.111/j.2517-6161.1979.tb01091.x

Romagnoli F, Rivals F, Stefano B (2022) Updating Neanderthals. Understanding behavioural complexity in the late middle Palaeolithic, 1st edn. Academic, London

Ronchitelli A, Boscato P, Surdi G et al (2011) The Grotta Grande of Scario (Salerno, Italy): archaeology and environment during the last interglacial (MIS 5) of the mediterranean region. Quatern Int 231:95-109. https://doi.org/10.1016/j.quaint.2010.07.006

RStudio C (2024) R Studio: integrated development for R. RStudio, Boston. https://posit.co/download/rstudio/

Saladié P, Rodríguez-Hidalgo A, Domínguez-Rodrigo M et al (2021) Dragged, lagged, or undisturbed: reassessing the autochthony of

the hominin-bearing assemblages at Gran Dolina (Atapuerca, Spain). Archaeol Anthropol Sci 13:65. https://doi.org/10.1007/s12520-021-01303-6

Sasowsky ID (2007) Clastic sediments in Caves – Imperfect recorders of processes in karst. Acta Carsologica 36:143-149. https://doi.org/10.3986/ac.v36i1.216

Sañudo P, Blasco R, Fernández Peris J (2016) Site formation dynamics and human occupations at bolomar cave (Valencia, Spain): an archaeostratigraphic analysis of levels I to XII (100 e 200 ka). Quatern Int 417:94-104. https://doi.org/10.1016/j.quaint.2015.09.044

Schoville BJ (2014) Testing a taphonomic predictive model of edge damage formation with middle stone age points from pinnacle point cave 13B and die Kelders cave 1, South Africa. J Archaeol Sci 48:84-95. https://doi.org/10.1016/j.jas.2013.10.002

Smith BA, Davies TM, Higham CFW (2015) Spatial and social variables in the Bronze Age Phase 4 cemetery of Ban Non Wat, Northeast Thailand. J Archaeol Sci Rep 4:362-370. https://doi.org/10.1016/j.jasrep.2015.10.003

Sánchez-Romero L, Benito-Calvo A, Marín-Arroyo AB et al (2020) New insights for Understanding Spatial patterning and formation processes of the neanderthal occupation in the Amalda I cave (Gipuzkoa, Spain). Sci Rep 10:8733. https://doi.org/10.1038/s41982-020-65364-8

Sossa-Ríos S, Mayor A, Hernández CM et al (2022) Multidisciplinary evidence of an isolated neanderthal occupation in abric Del pastor (Alcoi, Iberian Peninsula). Sci Rep 12:15883. https://doi.org/10.1038/s41598-022-20200-z

Sossa S, Alejandro R, Laura M et al (2025) Knapping ... sleeping and consuming ? Spatial variability in the High – Resolution neanderthal context of abric Del pastor (Alcoi, Eastern Iberia). J Palaeolithic Archaeol 8:23. https://doi.org/10.1007/s41982-025-0022-2-5

Spagnolo V, Aureli D, Martini I et al (2020) Short and close in time: overlapped occupation from the layer 56 of the molare rock shelter (Southern Italy). Archaeol Anthropol Sci 12:92. https://doi.org/10.1007/s12520-020-01037-x

Spagnolo V, Crezzini J, Falguères C et al (2024) Grotta Grande (southern Italy). Disentangling the neandertal and carnivore interaction in a short-term palimpsest at the last glacial onset (116-109 ka). Q Sci Rev 331:108628. https://doi.org/10.1016/j.quascirev.2024.108628

Spagnolo V, Marciani G, Aureli D et al (2016) Between hearths and volcanic ash: the SU 13 palimpsest of the Oscurusciuto rock shelter (Ginosa - Southern Italy): analytical and interpretative questions. Quatern Int 417:105-121. https://doi.org/10.1016/j.quaint.2015.11.046

Spagnolo V, Marciani G, Aureli D et al (2019) Neanderthal activity and resting areas from stratigraphic unit 13 at the middle Palaeolithic site of Oscurusciuto (Ginosa - Taranto, Southern Italy). Q Sci Rev 217:169-193. https://doi.org/10.1016/j.quascirev.2018.06.024

Stein JK (2001) A review of site formation processes and their relevance to geoarchaeology. Earth sciences and archaeology. Springer US, Boston, MA, pp 37-51

Todd LC, Standford DJ (1992) Application of conjoined bone data to site structural studies. In: Hofman JL, Enloe JG (eds) Piecing Together the Past: Applications of Refitting Studies to Archaeology. BAR International Series 578, Oxford, pp 21-35

Vaquero M, Fernández-Laso MC, Chacón MG et al (2017) Moving things: comparing lithic and bone refits from a middle paleolithic site. J Anthropol Archaeol 48:262-280. https://doi.org/10.1016/j.jaa.2017.09.001

Villaescusa L, Cobo-Sánchez L, Márquez B et al (2022) New contributions from the spatial statistical analysis of archaeological remains to the formation history of Level 3 at the Neanderthal
site Cueva Des-Cubierta Cave in Pinilla del Valle, Spain. Paleo
Anthropology, (2). p 182 https://doi.org/10.48738/2022.iss2.809

Villa P (1982) Conjoinable pieces and site formation processes. Am Antiq 47:276-290. https://doi.org/10.2307/279901

Villa P (2004) Taphonomy and stratigraphy in European prehistory taphonomy and stratigraphy in European prehistory. Before Farming 1:1–20. https://doi.org/10.3828/bfarm.2004.1.1

White WB (2007) Cave sediments and paleoclimate. J Cave Karst Stud 69:76-93

White WB, White EL (2003) Gypsum wedging and cavern break-
down: studies in the mammoth cave system, Kentucky. J Cave
Karst Stud 65:43–52

Wickham H (2016) ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org

Woodward JC, Goldberg P (2001) The sedimentary records in mediterranean rockshelters and caves: archives of environmental change. Geoarchaeology - Int J 16:327-354. https://doi.org/10.1002/geo.1007

Zilhão J, Banks WE, D'Errico F, Gioia P (2015) Analysis of site forma-
tion and assemblage integrity does not support attribution of the
Uluzzian to modern humans at Grotta Del Cavallo. PLoS ONE
10:e0131181. https://doi.org/10.1371/journal.pone.0131181

Zilio L, Hammond H, Karampaglidis T et al (2021) Examining neander-
thal and carnivore occupations of Teixoneres cave (Moià, Barcelona,
Spain) using archaeostratigraphic and intra-site Spatial analysis. Sci
Rep 11:4339. https://doi.org/10.1038/s41598-021-83741-9

**Publisher's note** Springer Nature remains neutral with regard to juris-
dictional claims in published maps and institutional affiliations.

Authors and Affiliations

Lucía Villaescusa¹,² · Enrique Baquedano²,³ · David M. Martín-Perea⁴ · Belén Márquez² ·

M. Ángeles Galindo-Pellicena¹,² · Lucía Cobo-Sánchez⁵ · Ana Isabel Ortega⁶ · Rosa Huguet⁷,⁸,⁹ · César Laplana¹,² ·

M. Cruz Ortega¹⁰ · Sandra Gómez-Soler²,¹¹,¹² · Abel Moclán³,¹³ · Nuria García¹⁰,¹⁴ · Diego J. Álvarez-Lao¹⁵ ·

Rebeca García-González¹⁶ · Laura Rodríguez¹⁶,¹⁷ · Alfredo Pérez-González³ · Juan Luis Arsuaga¹⁰,¹⁴

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

Lucía Villaescusa
lucia.villaescusa.fernandez@gmail.com

1 Arqueo-Paleo research group, Departamento Geología, Geografía y Medio Ambiente, Universidad de Alcalá, Alcalá de Henares, Spain

2 MARPA - Museo Arqueológico y Paleontológico de la Comunidad de Madrid, Alcalá de Henares, España

³ Institute of Evolution in Africa – IDEA, Madrid, Spain

⁴ Institute of Vertebrate Paleontology and Paleoanthropology, Chinese Academy of Sciences, Beijing, China

5 ICArEHB – The interdisciplinary Center for Archaeology and Evolution of Human Behaviour, Universidade do Algarve, Campus de Gambelas, Faro 8005-139, Portugal

6 Real Academia Burgense de Historia y Bellas Artes, Institución Fernán González, Burgos, España

⁷ IPHES-CERCA, Institut Català de Paleoecologia Humana i Evolució Social, Tarragona, España

8 Departament d'Història i Història de l'Art, University Rovira i Virgili, Tarragona, España

9 Unit associated to CSIC, Departamento de Paleobiología, Museo Nacional de Ciencias Naturales, Madrid, Spain

¹⁰ UCM-ISCIIII Research Center of Human Evolution and Behavior, Madrid, Spain

11 Archaeological Micromorphology and Biomarkers Laboratory (AMBi Lab), Instituto Universitario de Biologíaica "Antonio González", La Laguna University, Tenerife, Spain

12 Department of History and Philosophy, University of Alcalá de Henares, Madrid, Spain

13 PALEVOPRIM Lab, UMR 7262, CNRS and Université de Poitiers, Potiers, France

14 Department of Geodynamics, Stratigraphy and Palaeontology, Faculty of Geology, Complutense University of Madrid, Madrid, Spain

15 Departamento de Geología, Universidad de Oviedo, Oviedo, España

16 Laboratorio de Evolución Humana, Facultad de Humanidades y Comunicación, Universidad de Burgos, Burgos, España

17 Área de Antropología Física, Facultad de Cc. Biológicas y Ambientales, Universidad de León, León, España