Introduction

The problem of dissecting palimpsests in archaeological sites (Bailey 2007; Henry 2012) and the consequent objective of identifying indicators of diachrony in the archaeological record have become central issues in current research.

In the framework of palimpsest dissection, the stratigraphic association of faunal assemblages and lithic industry allows us to identify patterns in their composition. However, the field layers in Pleistocene palimpsests deposits may appear coherent during excavation but can often belie the underlying complexity of the archaeological record.

Processes of cave and rock shelter formation, such as erosion or sedimentation rates, can make it difficult to identify and define the periods of occupation that occurred in the same geological space and time or stratigraphic unit. But also, there are excavation-related factors like the “analytical lumping” (Discamps et al. 2019), which refers to the inadequate separation of distinct assemblages during or after excavation and analysis.

Discamps et al. (2023) proposed the role of Post-Excavation Stratigraphies (PES) as a promising approach to critically review and redefine field layers, incorporating spatial data, and considering the dynamic nature of archaeological assemblages.

In this sense, it is important to bear in mind that taphonomic processes, such as mass movements or water flows, can cause the displacement of materials from their original position of abandonment and contribute to the generation of mixed associations (Voorhies 1969; Enloe 2006; Goldberg and MacPhail 2006; Domínguez-Rodrigo et al. 2018). In this way, elements from different biological entities or from different time periods will be found together, which may lead to misleading associations and interpretations of archaeological assemblages.

Thus, it is important to identify the temporal dimension of past human behaviour, together with the processes of formation and transformation of archaeo-sedimentary deposits, in order not to condition the analysis and interpretation of the archaeological records and their relationship with the sediments that contain them (Leroi-Gourhan and Brézillon 1972; Farrand 1975; Binford 1981b; Ferring 1986; Schifer 1987; Shott 1998; Straus et al. 2001; Holdaway and Wandsnider 2008; Finlayson et al. 2008; Lucas 2012). In order to solve these problems, the theoretical bases of archaeo-stratigraphy were established in the 1980s.

Archaeo-stratigraphy attempts to isolate different events or sublevels related to human occupation within a single geological deposit (Canals 1993, 1996; Canals et al. 2003). The main objective is to identify and order the events of accumulation and abandonment of materials involved in the formation of a particular site (Goldberg and MacPhail 2006; Vaquero et al. 2012; Machado et al. 2013; Mallol et al. 2013).

Several studies have suggested the diachronic relationship within lithostratigraphic units and the need to apply high-resolution methodology for the analysis of these units (Carr 1987; Audouze and Enloe 1997; Shott 2008; Vaquero 2008, 2013; Henry 2012; Vaquero et al. 2012, 2015; Sañudo et al. 2012; Rosell et al. 2012a; Machado et al. 2015; Mallol y Mentzer 2015; Modolo and Rosell 2016; Romagnoli et al. 2018; Martín-Perea et al. 2020; Arteaga-Brieba et al. 2023).

In recent years, different approaches have been used to improve the temporal resolution of archaeological analyses. One of the most common is the use of archaeo-stratigraphy in combination with the refitting, spatial and trace study of the lithic industry (e.g. Canals et al. 2003; Bourguignon et al. 2008; Vaquero 2008; Vaquero et al. 2012, 2017; Machado et al. 2013, 2019; Sumner and Kuman 2014; Romagnoli and Vaquero 2016; Sañudo et al. 2016; Nerudová and Neruda 2017; López-Ortega et al. 2017, 2019; Takakura 2018; Mosquera et al. 2018; Moncel et al. 2021; Hammond et al. 2022; Arteaga-Brieba et al. 2023; Moclán et al. 2023a, 2023b).

Archaeo-stratigraphy also can be combined with zooarchaeological, taphonomic, technological, raw material, and dental wear data (Bleed 2002; Chacón et al. 2015; Pérez-Diaz and López-Sáez 2021; Malinski-Buller et al. 2021; Pieruccini et al. 2022). The integration of faunal remains as part of archaeo-stratigraphic studies has demonstrated the potential of using taxonomic and taphonomic analysis of faunal remains to narrow down the time scale of occupation levels at archaeological sites (Rivals et al. 2009; Schoville and Otárola 2014; Pérez et al. 2015; Rodríguez et al. 2016; Gabucio et al. 2018; Saladié et al. 2018). Some works also analyse bone refitting with spatial and archaeo-stratigraphic perspectives (e.g. Rosell et al. 2012a, b; Chacón et al. 2015; Pérez et al. 2015, 2017; Gabucio et al. 2018; Machado and Pérez 2016; Fernández-Laso et al. 2020).

Site formation and taphonomic processes affecting assemblages can be significant factors in distinguishing between biological and geological time scales. (i.e. Goldberg and Berna 2010; Vallverdú and Courty 2012; Henry 2012; Mallol et al. 2012; Vallverdú 2013a, b; Tsatskin and Zaidner 2014; Polo et al. 2016; Spagnolo et al. 2016). Therefore, we highlight the importance of archaeo-stratigraphic analyses based on the taphonomic analysis of skeletal remains to analyse geologically homogeneous sedimentary units, which includes the definition of analytical units that allow us to explore the formation of accumulations in as much detail as possible.

The Upper Pleistocene levels of Buena Pinta Cave have been interpreted as a hyena den (Huguet et al. 2010; Baquedano et al. 2012, 2016). However, the presence of lithic elements in the western part of the site (level 23) together with bone remains with apparent anthropogenic fracture (Baquedano et al. 2012) suggests a carnivore accumulation with the presence of sporadic human activity (Baquedano et al. 2012, 2016).

Until now, level 23 of Cueva de la Buena Pinta has been defined during fieldwork as a homogeneous package, and the archaeological remains have been subsequently assigned to this field layer. However, the heterogeneity of taphonomic alterations observed during the study of the faunal materials and the presence of pebbles and lithic industry suggest the palimpsestic nature of level 23, including geological and/or biological disturbances.

The main objective of this work is to evaluate the integrity of level 23. Thus, following the PES idea of Discamps et al. (2023), we re-evaluate the western part of the site of Buena Pinta Cave (Pinilla del Valle, Madrid) using a taphonomic perspective to carry out intra-site spatial analysis of the bone remains, coprolites, and lithic remains. The integration of quantitative and qualitative data allows us to demonstrate in a practical way the need to reassess the field layers and interpret spatial patterns in the distribution of the archaeological remains thanks to taphonomy.

Buena Pinta Cave

Buena Pinta Cave was discovered in 2003 and has been excavated continuously until 2022, except in 2020 and 2021 due to the COVID-19 pandemic. It is part of the Calvero de la Higuera archaeo-palaeontological complex (Fig. 1A) near Pinilla del Valle (altitude 1114 m a.s.l.) in the upper reaches of the Lozoya river valley in the Sierra de Guadarrama (in the Sistema Central) about 55 km north of Madrid (Pérez-González et al. 2010; Huguet et al. 2010; Arsuaga et al. 2010, 2012; Baquedano et al. 2012, 2016, 2021, 2023; Márquez et al. 2013, 2016, 2017; Abrunhosa et al. 2020; Laplana et al. 2015a, b; Arriaza et al. 2017; Moclán et al. 2018, 2020, 2021, 2023b). The caves and rock shelters in this complex (Fig. 1B) are associated with cavities that developed in Upper Cretaceous dolomitic rocks, subsequently occupied by Neanderthals during isotopic stages 5, 4, and 3 (Pérez-González et al. 2010).

Fig. 1
figure 1

Location of Calvero de la Higuera. A Map of the Iberian Peninsula and topographic map (IGN) showing the location of the sites. B Orthophoto of the Neanderthal Valley sites (1, Camino Cave; 2, Navalmaíllo Rockshelter; 3, Buena Pinta Cave; 4, Des-Cubierta Cave). C Aerial view of the sites to the east ©Javier Trueba (Madrid Scientific Films). D Stratigraphic column of the western part of Buena Pinta Cave

Full size image

Buena Pinta Cave (Fig. 1C) is a cavity of phreatic origin (formed below the water table) which is divided into two main units (Baquedano et al. 2021). A Holocene deposit of approximately 1.80 m in the upper part of the sequence blocks the entrance to the gallery where the remains of a Bronze Age burial were recovered. An Upper Pleistocene deposit in the lower part of the sequence is composed of several levels. In this lower unit, comprised of levels 2 to 5, level 3 (Fig. 2) was dated by TL at 63.4 ± 5.5 ka BP (Pérez-González et al. 2010). Based on the biochronological interpretation of the microfaunal assemblage from levels 2 to 5, in conjunction with this data, it can be inferred that the sequence corresponds to the middle of the Late Pleistocene, specifically ranging from MIS4 to early MIS3 (Laplana et al. 2016). Although the stratigraphy of the entire cave has not yet been well defined (in progress), a great complexity of facies has been revealed.

Fig. 2
figure 2

Top: plan of Buena Pinta Cave with the faunal materials excavated between 2003 and 2019. Bottom: XZ section of the faunal remains. In blue the faunal remains from level 23. In green, those from level 3. In black, the fauna of the remaining levels

Full size image

A partial zooarchaeological and taphonomical study of the fossil assemblage from level 3 of the cave (Fig. 2), together with the presence of numerous coprolites, has identified this level as a hyena den (Huguet et al. 2010; Baquedano et al. 2012, 2016). Some evidence of the presence of hominins has also been documented in level 3, such as lithic pieces and two molars attributed to Homo neanderthalensis (Baquedano et al. 2010; Huguet et al. 2010).

The lithic remains are mainly concentrated in the western part of the site (level 23). Most of the pieces are made of quartz, and unretouched flakes and fragments are abundant (Baquedano et al. 2012). Among the retouched pieces, denticulates dominate, and the cores are knapped unifacially, bifacially, or trifacially, usually in a unipolar longitudinal manner.

Geologically (Fig. 1D), the archaeological assemblage of level 23 is embedded in a 1-m-thick silt deposit. Towards the bottom of this silt deposit, there are abundant very small and small dolostone boulders (following Blott and Pye 2012), which are sub-rounded and sub-angular, together with sub-rounded and rounded very small and small gneiss, granite, quartzite, and slate boulders. The base of this deposit is erosive and diffuse, where silts and clasts are mixed with underlying sediments and small to medium-sized dolostone boulders and speleothem fragments. The underlying unit is a 30-cm-thick reddish-brown clay deposit with scarce very small and degraded dolostone boulders. During fieldwork, archaeological materials recovered from this unit were preliminarily assigned to level 23 given its meagre extent in the excavation.

Level 23 was defined during fieldwork limited to what can be observed in situ. The area of clay deposit was very limited, and it was not possible to know whether it really corresponded to a different level or to an intrusion of clays within level 23. Therefore, tentatively all the materials collected in the field were ascribed as level 23.

Materials and methods

Materials

In this work, a total of 2915 faunal remains assigned in the field to level 23 (Fig. 2) have been examined. All this material was collected during the 2009–2019 excavation campaigns and includes all the mapped remains. All faunal remains greater than 2 cm along their longest axis were recorded, as well as smaller identifiable remains of taphonomic and/or palaeontological interest (e.g. bone flakes, teeth). All lithic remains and dolostone blocks larger than 20 cm were also mapped.

Methods

Zooarchaeology and taphonomy

All skeletal remains have been measured in tenths of millimetres with a digital calliper. They have been identified anatomically and taxonomically using osteological atlases (e.g. Pales and Lambert 1971; Schmid 1972; Hillson 2005) and the comparative anatomy collection of the Institut Catalá de Paleoecologia Humana i Evolució Social (IPHES).

Three death age groups have been considered: immature, adult, and senile, taking into account the degree of epiphyseal fusion in the long bones, cortical development, growth of the dental crown, identification of deciduous or permanent teeth, and the wear of the occlusal surface of the teeth (Klein and Cruz-Uribe 1984; Barone 1986).

The number of remains (NR), the number of identified specimens (NISP), the minimum number of elements (MNE), and the minimum number of individuals (MNI) were calculated according to Lyman (1994a, 2008), Klein and Cruz-Uribe (1984), Brain (1969), Pickering (2002), and Yravedra and Domínguez-Rodrigo (2009).

The analysis of bone surface modifications was carried out with a Euromex StereoBlue SB.1902-P, 6.7x and 45x binocular loupe. A HIROX KH-8700 digital microscope was also used when more detailed analysis was required.

The state of preservation of the bone surfaces has been determined on a scale of six grades: (0) 0%, (1) 1–25%, (2) 25–50%, (3) 50–75%, (4) 75–99%, and (5) 100% cortical preservation (Moclán et al. 2021). Anthropogenic modifications such as cut and percussion marks have been identified (Potts and Shipman 1981; Blumenschine and Selvaggio 1988; Blumenschine 1995; Galán and Domínguez-Rodrigo 2013). Carnivore activity has been identified by the presence of tooth marks such as pits, punctures, scores or furrowing, and by digested remains (Haynes 1980, 1983; Binford 1981a; Selvaggio 1994; Pickering 2001). In the case of carnivores, according to Fernández-Jalvo and Andrews (2016), enzymes and acids generate rounding and polishing in specific areas of the bone as well as sharpening of the edges if the bone is digested. In addition, in the case of polishing by carnivores, when it appears in specific areas of the bone, it is generally associated with other alterations produced by these agents (e.g. gnawing, furrowing).

Weathering was analysed according to Behrensmeyer (1978), which ranges from stage 0 (not weathered, exposed less than 1 year before burial) to stage 5 (extremely weathered, exposed for 15–30 years before burial). Water-related modifications have been identified, such as abrasion and cortical loss (Fernández-Jalvo and Andrews 2016). Water affects all the remain and generates rounding in the whole bone or in the fracture zones, with polishing all over the altered surface. Rounding and polishing have been grouped into three different grades, according to Bouchud (1974), Shipman and Rose (1983), and Cáceres (2002): (1) observable microscopically, (2) observable macroscopically, and (3) significantly affects the whole bone.

Other biostratinomic and fossil-diagenetic alterations have been documented as either present or absent, such as trampling (Behrensmeyer et al. 1986; Olsen and Shipman 1988) or biochemical marks (Morlan 1980; Andrews 1995; Domínguez-Rodrigo and Barba 2006; Backwell et al. 2012). The presence/absence of calcareous and carbonate concretions has also been taken into account. Manganese oxides (Courty et al. 1989; Coard 1999) have been documented in terms of their presence/absence, and their dispersion has been recorded according to Moclán (2023).

Archaeo-stratigraphy

For the archaeo-stratigraphic study, we considered all materials mapped in level 23 of Buena Pinta Cave. Of the 3872 remains mapped, 2915 are bone remains, 665 lithic remains, 79 coprolites, and 213 dolostone blocks.

The zooarchaeological and taphonomic characteristics of the mapped bone remains have been used to carry out the archaeo-stratigraphic analysis. Lithic remains data have been incorporated as they provide basic information on the human presence in the cave and the patterns of anthropic occupation of the cave. For this purpose, the raw material, the origin of the material (i.e. natural or modified), and the spatial coordinates of the remains have been taken into account.

The spatial coordinates were recorded using a motorized Leica total station (model TCRP 1205 R400). Three-dimensional geographic information system (GIS) software models (with ArcGIS® ArcMap 10.5) were used to perform Kernel density analyses and spatial projections. All materials have been recorded as individualized elements in a three-dimensional diagram representing the deposit using Cartesian coordinates.

From the spatial data, different plans and sections have been created to facilitate visualization. The projections presented will focus on the XZ and YZ sections. These maps have been used to determine whether the distribution of archaeological materials was uniform throughout the site or whether they were concentrated in specific sectors according to the taphonomic alterations identified.

Statistical analyses

The materials included in this work have been analysed from a spatial-statistical point of view, for which the statistical software R (R Core Team 2022) has been used, along the lines of what has been proposed in other recent works (e.g. Diez-Martín et al. 2021; Domínguez-Rodrigo et al. 2017, 2018; Domínguez-Rodrigo and Cobo-Sánchez 2017; Luzón et al. 2021; Marín et al. 2019; Saladié et al. 2021).

The R “spatstat” library (Baddeley and Turner 2005) has been used to spatially analyse the fossil record. The Kernel Smoothed Intensity function, which identifies those areas where intensity is highest by means of a heat map, was used. For the calculation of the intensity maps, Likelihood cross-validation methods were used. To identify those areas where the material is statistically significant, the “scanLRTS” command based on “poisson” probability has been used (Baddeley et al. 2015). For this analysis, the alternative hypothesis “greater” has been used. To further refine the analysis, a hot spot analysis has been carried out, taking into account those areas in which the probability of the presence of archaeological material is significant at 99%.

The libraries “ggplot2” (Wickham 2016) and “patchwork” (Lin Pedersen 2020) have been used to obtain graphs showing the vertical distribution of materials.

The “spatstat” library (Baddeley and Turner 2005) has been used to work with a spatial pattern of points (SPP), both single category (i.e. unmarked patterns) and multi-category (e.g. uncut cobbles vs. lithic industry), known as “marked spatial patterns” (Baddeley et al. 2016).

For unmarked SPPs, the homogeneity of the samples was tested with a χ2 test, dividing the analysed surface into 5 × 5, 8 × 8, 10 × 10, and 12 × 12 grids (see Moclán et al. 2023a, b). Given the inhomogeneity of all samples (χ2 p value ≤ 0.05), a Studentised Permutation test using an inhomogeneous K-function and a K-scaled function was used to characterise the inhomogeneity of the SPPs. To carry out this analysis, the spatial window (i.e. excavation area) has been tessellated in a 3 × 4 grid (see Moclán et al. 2023a, b).

The results obtained show that the scaled functions cannot be used to analyse the SPP of level 23’s faunal remains (section XZ, p value = 0.001; section YZ, p value = 0.001). However, the results of the Studentised Permutation test with the inhomogeneous K-function have shown that it is feasible to perform the analysis with the inhomogeneous functions (section XZ, p value = 0.315; section YZ, 0.691).

Two tests were performed in order to identify whether the distribution of the analysed SPPs is random or not. First, Hopkins-Skellam and Clark-Evans tests were carried out, which allow to differentiate between random, cluster, and regular SPPs. However, these tests tend to confuse inhomogeneity with cluster patterns, so their results should be evaluated with caution (Baddeley et al. 2016). Diggle-Cressie-Loosmore-Ford (DCLF) and Maximum Absolute Deviation (MAD) tests were used to identify whether an SPP corresponds to a random model or not. Both tests (MAD and DCLF) have been performed using the inhomogeneous K, L, F, G, and J functions. In addition, to test the null hypothesis of a random model against the alternative hypotheses of cluster or regular patterns, the K, L, F, G, and J functions in their inhomogeneous versions were used, as well as the Pair Correlation Function, given the result of the Studentised Permutation test.

Finally, the “nnclean” function of the “spatstat” library has been used to perform a Nearest-neighbour cleaning test with the faunal SPPs with k values of 10, 20, 50, and 80 to check the type of distribution of the sample in both sections. This test assumes the presence of a cluster pattern, so it should be used in an exploratory way in case the distribution is not of this type with the intention of identifying possible associations of material.

The homogeneity of SPPs with marks have been analysed using a χ2 (Baddeley et al. 2016). Next, the probability with which some marks may appear in a given space compared to others has been analysed using a Spatially-Varying Relative Risk test. To analyse the statistical correlation, K· and L· functions have been used when only two marks are compared, and the Kcross and Lcross functions were used to compare two or more marks simultaneously. These analyses aim to measure the type of spatial correlation between faunal remains, lithic industry, and pebbles. We have used only the inhomogeneous versions of the functions, given the results obtained with the χ2.

Results

Zooarchaeology and taphonomy

A total of 1353 specimens (26.8%) assigned to level 23 have been taxonomically identified and correspond to 25 different species, with deer (%NISP = 16.01) and bovines (%NISP = 13.81) being the most abundant. The taxonomic diversity of this assemblage is high, with 11 species of ungulates (e.g. Stephanorhinus hemitoechus, Equus ferus, Sus scrofa, Capreolus, Rupicapra, proboscidean) and 9 carnivoran species (e.g. Crocuta, Panthera pardus, Canis lupus, Vulpes vulpes). Smaller taxa also appear, such as Oryctolagus cuniculus, Testudines, and Erinaceus sp.

The two proboscidean tusk fragments documented in this fossil assemblage are the only records of this taxon documented in the whole of the Calvero de la Higuera complex so far.

We identified 637 elements that correspond to at least 75 individuals, adults being the most represented (63%). Thirty-six percent of the individuals correspond to immature individuals, with the presence of two perinatal individuals of deer standing out. Finally, 2 senile individuals have been recorded.

In one-third of the sample (n = 966, 33.14%), some alteration related to carnivore activity was observed (Supplementary File 1C and 1D), the most frequent modifications being grooves and pits. At the same time, anthropogenic action in the form of cut and percussion marks was only recorded in 0.9% (n = 26) of the sample.

Large differences are observed in relation to fossil-diagenetic alterations. While 45.93% (n = 1339) of the assemblage has a good preservation of the cortical surface, about 35% (n = 1020) has almost no cortical preserved. Eighty-eight percent (n = 2565) of the assemblage shows no weathering, but some remains are affected in the most advanced degrees (Supplementary File 1F).

A total of 1877 remains (64.39%) show water-related modifications. And 39.37% have suffered cortical loss due to the abrasion and exhibit a bleaching of the cortical (Supplementary File 1E). Then, 34.42% of the sample is rounded, affecting all sizes and all grades (Table 1, Supplementary File 1A and 1B). Polishing, on the other hand, affects 26.21% of the remains, with a wide range of development (Table 1, Supplementary File 1A). Manganese oxide precipitates are observed in 81.05% of the sample (n = 2363) in all weight groups, as well as calcareous concretions (n = 1335, 45.78%). In addition, 13 fragments of the fossiliferous breccia have been documented (Fig. 2), with carbonate concretions.

Table 1 Number of remains according to their degree of rounding and polishing. The % are based on the total number of remains in the sample
Full size table

Spatial statistics

The assemblage is concentrated in 6 m2 and is located in the western part of the site. Although excavation work is still in progress, there are areas with a higher density of material (Fig. 3).

Fig. 3
figure 3

A Projection of the faunal remains (orange) from level 23 of Buena Pinta Cave in the three views. B Section YZ (green). C Section XZ (blue). The dotted line marks the unexcavated areas. D Plan layout in purple (XY)

Full size image

Taking into consideration the inclination with which the sediments have been deposited, we have avoided carrying out spatial analyses in plan, as this would generate an altered image of reality. Therefore, all the analyses of the faunal distribution patterns have been carried out in both XZ and YZ sections.

The intensity maps show several zones of high intensity in the XZ and YZ sections (Fig. 4). These peaks are separated by zones where material density is much lower, suggesting the possible presence of periods of lower material accumulation. This separation is most evident where material accumulation is significant at 99%, with the highest density in the lower part of the sedimentary unit.

Fig. 4
figure 4

Faunal remains from level 23 of Cueva de la Buena Pinta in section XZ (top) and YZ (bottom). From left to right: point pattern scaled to unit square (A, D); Kernel estimate of density, with bandwidth selected by likelihood cross-validation (B, E); p values of the likelihood ratio test statistic where p < 0.01, with logarithmic colour scale (C, F)

Full size image

The results of the analyses carried out on section YZ (Fig. 4 D, E, and F) are more illuminating than those of section XZ, due to the natural westward slope over which the strata have formed. Thus, in Fig. 4F, different accumulations of material can be seen in the lower northern half, whose horizontal continuity is lost towards the south. In addition, two features are worth highlighting: a northward sloping distribution of the remains in the last lower third (Fig. 4E–F) and the contrast between the density of material between the lower right half and the lower left half, where the amount of remains decreases drastically towards the wall of the excavated surface.

After verifying the existence of several areas of high density, the type of distribution of the sample was analysed. In section XZ and YZ, the χ2 test indicates that the level 23 faunal assemblage is inhomogeneous (Supplementary File 2), regardless of whether the analysed surface is divided into 5 × 5, 8 × 8, 10 × 10, or 12 × 12 grids.

The Hopkins-Skellam (A = 0.10193, p value < 2.2e-16) and Clark-Evans (R = 0.92838, p value = 0.002) tests in section XZ and Hopkins-Skellam (A = 0.15439, p value < 2.2e-16) and Clark-Evans (R = 0.94366, p value = 0. 002) in section YZ identify the faunal SPPs as a cluster distribution. However, given the inhomogeneity of the SPPs, we must consider these results as probably erroneous and more in support of the result obtained by the χ2.

The results obtained from the application of the DCLF and MAD tests on the XZ and YZ sections have shown that when calculated using the inhomogeneous K and L functions, the null hypothesis of random pattern can be rejected (p value = 0.025), while the same tests with the G, F, and J functions failed to reject the H0 (p value ≥ 0.05).

Once applied, Fig. 5 shows how the inhomogeneous K and L functions suggest a cluster trend. In both the XZ and YZ sections, the line showing the trend of the SPPs pattern diverges above the inhomogeneous Poisson model and overhangs the confidence interval (envelope) in some sections. In the XZ section, it is more evident, although in the K functions of both sections, a regular distribution is observed up to about 0.15 radius, where the line separates at the top to indicate a clustering trend. The L function shows less clearly the tendency to cluster, especially in the YZ section, as the line is less separated from the envelope and the theoretical Poisson model (Fig. 5). On the other hand, functions F, G, J, and Pair correlation (Fig. 5; Supplementary Files 3–4) do not show conclusive results that contribute to determine whether the distribution of SPPs tends to a positive association between them (cluster) or to segregation.

Fig. 5
figure 5

Pair correlation function and inhomogeneous functions K and L for the spatial point pattern of faunal remains in section XZ (top) and YZ (bottom)

Full size image

As the results seem to indicate a tendency for the distribution to be clustered in some parts of the SPP, a Nearest-neighbour cleaning test has been carried out (Fig. 6). It identifies most of the surface of both sections as clusters, leaving out a fraction of the sample. However, the agglomeration in the lower half of XZ and the presence of certain areas with lower cluster probability in YZ, especially at k values of 10 and 20, are striking (Fig. 6 A and B).

Fig. 6
figure 6

Results obtained by Nearest-neighbour cleaning test on the set of faunal remains (feature) from level 23 of Buena Pinta Cave in section XZ and YZ with k values of 10 (A), 20 (B), 50 (C), and 80 (D)

Full size image

Archaeo-stratigraphy

Based on the intensity analyses and on the statistical functions K and L, we can confirm the presence of faunal remains with a weak positive association in the so-called level 23 of the Cueva de la Buena Pinta together with other faunal remains with a regular distribution (negative association). Given that it is assumed that the nearby specimens are likely to have followed the same principles of formation and alteration (Fernández-López 2000), the spatial analysis of the remains has been carried out based on their taphonomic characteristics in order to know the distribution of the remains and to be able to delimit and characterise these materials.

Although the most abundant alterations observed in the analysed assemblage are manganese oxides (n = 2363, 81.05%) and biochemical dissolutions derived from plant roots (n = 1801, 61.76%), no groupings of bone remains are observed based on their presence/absence. The same is true for trampling or tooth marks (Fig. 7), which are evenly distributed throughout the assemblage. The intensity and probability maps at 99% confidence (Supplementary Files 5 and 7) show several intensity foci, confirming a distribution of remains altered by these taphonomic processes throughout the entire assemblage.

Fig. 7
figure 7

Faunal remains from level 23 of Buena Pinta Cave in section XZ (left) and YZ (right) projected according to the presence (green) and absence (grey) of manganese oxides, root marks, trampling, and tooth marks

Full size image

However, the projection of the fauna based on polishing (Fig. 8) shows different dynamics. Some degree of polishing has been recorded in 26.41% (n = 770) of the remains, with 5.22% of them (grades 2 and 3) mainly concentrated below a depth of 2.8 m.

Fig. 8
figure 8

Faunal remains (black dots) from level 23 of Buena Pinta Cave and their intensity based on the degree of polishing in section XZ (A, B, C, D) and YZ (E, F, G, H)

Full size image

In addition, although there are remains with low degrees of polishing throughout the sequence, another accumulation of polished remains, this time of lesser intensity, can be seen 2.6-m deep.

Based on the degree of polishing suffered by the remains and their location in the sequence (Fig. 8), we propose the existence of two distinct associations: one above − 2.5 m with a predominance of unpolished remains alongside some with grades 1 and 2 (Fig. 8 A, B, E, and F) and another below this level where, despite finding unpolished remains, a greatest number of bones with the highest degrees of alteration are concentrated, especially below – 2.8 m (Fig. 8 C, D, G, and H).

A similar pattern is observed if rounding is taken into account. While undisturbed bones can be seen throughout the entire sequence, rounded bones are especially concentrated in the lower part of the assemblage.

Given the number of remains affected by rounding in the whole assemblage (n = 1352, 46.38% of the sample) and their dispersion, generating a division is more complicated. However, considering the proportion of remains with a high degree of alteration and the areas with a lower density of remains, especially in the YZ section, two sets can be differentiated. One above – 2.5 m where 19.69% (n = 267) of the rounded remains are found (of which 87% are grade 1, Fig. 9 A, B, E, and F) and the other below, where 36.05% (n = 488) of the rounded remains are concentrated, with 15.28% showing grade 2 modifications (Fig. 9 C and G) and 4.50% grade 3 (Fig. 9 D and H).

Fig. 9
figure 9

Faunal remains (black dots) from level 23 of Buena Pinta Cave and their intensity based on the degree of rounding in section XZ (A, B, C, D) and YZ (E, F, G, H)

Full size image

In the area between 50.8 and 51.1 along the Y axis (Fig. 9), a set of highly rounded remains is observed with a diagonal trend, which then stabilises horizontally. These modified remains follow a wedge shape derived from the erosive base of the conglomerate over the underlying unit of reddish-brown clays (Fig. 2).

The remains with calcareous concretions appear throughout the assemblage, affecting 45.78% (n = 619) of the bones. Despite this, 54.7% of the affected remains are concentrated at a depth of 2.8 metres and above (Fig. 10), coinciding with a focus of intensity (Supplementary Files 6 and 8). In the whole assemblage, only 16 bones have carbonate concretions, and they are concentrated in the same sector, at a height of between − 2.3 m and − 2.55 m.

Fig. 10
figure 10

Faunal remains from level 23 of Buena Pinta Cave in section XZ (left) and YZ (right) projected according to the presence or absence of concretions (top) and water-derived modifications (bottom)

Full size image

On the other hand, the bones altered by water follow a similar distribution to the rounded and polished bones, separated into two units with a boundary around − 2.5 m (Fig. 10). Thus, the upper part (between − 1.6 and − 2.5 m) is dominated by unmodified bones, with only 6.74% of the bones altered by water. Finally, there is a lower zone where 20.5% of the sample modified by water is concentrated, especially from a depth of 2.8 m down.

Finally, the weathered bones (n = 318) are distributed throughout the sedimentary unit (Fig. 11), especially those in stage 2. Bones in an early stage of weathering are concentrated in the lower zone, between 2.6 m and 3-m deep. In contrast, the only remains with a longer time of subaerial exposure (stage 4) are located in the upper half of the assemblage, 2.6 m beneath the surface.

Fig. 11
figure 11

Faunal remains from level 23 of Buena Pinta Cave in section XZ (left) and YZ (right) projected according to weathering degree

Full size image

Vertical distribution analysis

If the taphonomic modifications are just a random effect, then we should expect more modified remains in areas with more bones. Based on the density of the whole sample (Fig. 12A), this is the case, for example, of root marks and manganese oxides (Fig. 12B and C). However, another dynamic can be observed in the vertical distribution around − 2.5 m, where the remains unaltered by taphonomic processes, such as the water-related modifications, the formation of calcareous concretions, or rounding and polishing, begin to decrease, and the density of more altered fossils increases towards the bottom of the level.

Fig. 12
figure 12

Vertical distribution of the faunal remains (A) and the density of fossils altered by: root marks (B), manganese oxides (C), water (D), concretion (E), polishing degree (F), and rounding degree (G)

Full size image

The density values of the water-modified and concreted remains are quite similar, with a peak of higher density of affected material at around − 2.8 m (Fig. 12D and E). In the upper half, polished bones are scarce. At around − 2.5 m, the density of more advanced polishing grades (grade 2 and 3) increases significantly, with a peak at around − 2.8 m (Fig. 12F).

For rounded bones, from − 2.5 m downwards, there is again a reversal of values revealing a change in depositional dynamics (Fig. 12G). The remains without evidence of rounding decrease at the same point where fossils with some degree of alteration increase. There is a clear predominance of remains modified by post-depositional processes in the deeper areas of the sequence, drawing a distribution with a maximum peak at around − 2.8 m. This contrasts greatly with the more normally distributed disposition of unmodified or slightly modified bones (Figs. 12F and G).

The vertical distribution of the density of weathered remains follows a completely different pattern to those presented above. Unweathered bones or bones with evidence of early stages of weathering (stages 1, 2, and 3) are more evenly distributed throughout the assemblage (Fig. 13).

Fig. 13
figure 13

Vertical distribution of fossil density according to the degree of weathering

Full size image

Nevertheless, remains with weathering stages 4 and 5 are most abundant between − 2.30 m and − 2.55 m approximately (Fig. 13). This indicates that bones have been exposed a long period of time before burial, which suggests the existence of a sedimentary hiatus.

The vertical analysis of post-depositional processes shows clear distribution patterns and reveals that the processes responsible for these distributions have acted differently throughout the stratigraphic sequence.

The graphs showing the vertical distribution of fossil density suggest a positive spatial association between the remains according to the type and degree of alteration. This concentration is evident in the lower half, from 2.5 m to its maximum peak at around − 2.8 m Z. The most heavily weathered remains are concentrated between 2.30 and 2.55-m deep, which marks a boundary with the upper half that was less altered by fossil-diagenetic processes.

Spatial analysis of coprolites

In the western part of Buena Pinta Cave, 79 coprolites were recovered throughout level 23 (Supplementary File 1G). Section XZ (Fig. 14A) shows how the coprolites are randomly dispersed throughout the level between − 2.6 m and − 3 m. However, around − 2.5 m, a horizontal line of coprolites can be traced. In section YZ (Fig. 14B), a 30-cm-wide accumulation of coprolites is also observed at − 2.5 m.

Fig. 14
figure 14

Spatial projection of mapped coprolites and its intensity analysis in the western part of Buena Pinta Cave in section XZ (A, C) and YZ (B, D)

Full size image

Lithic assemblages and raw material

The spatial analysis of the lithic elements indicates that there is a focus of intensity in the lower part of the assemblage (Fig. 15A, C, and D), between − 2.6 m and − 3 m. Twelve different raw materials have been identified. Most of them come from semi-local source areas, such as gneiss, quartz, flint, sandstone, porphyry, or granite. However, there are also pebbles from distant sources, such as quartzite and slate, originating from outcrops located several kilometres upstream of the Lozoya River.

Fig. 15
figure 15

A Spatial projection of the mapped lithic remains in the western part of Buena Pinta Cave. B Spatial projection of lithic remains eroded or not by water. C, D, E, and F areas where p < 0.01

Full size image

Among the lithic remains, pebbles predominate over lithic tools. On the other hand, 138 remains, both pebbles and tools, show patinas and erosion derived from the action of water (Fig. 15B). However, the highest density areas for remains with and without erosion overlap (Fig. 15E and F), both located in the lower half of the level. This alteration was seen in 23.17% of the lithic tools and in 17.46% of the pebbles.

In the analysis of the type of distribution of the lithic remains as SPPs with various categories or marks, the χ2 test indicates that the assemblage formed by the pebbles and lithic industry is inhomogeneous, regardless of whether the analysed surface is divided into 5 × 5, 8 × 8, 10 × 10, and 12 × 12 grids (Supplementary File 9).

The K· and L· functions (Supplementary File 10) have not yielded any significant results. On the other hand, in section XZ, the separation of the line below the envelope in the Kcross and Lcross functions from a radius of approximately 0.1 m indicates a negative correlation between the lithic remains with and without evidence of anthropic activity (Fig. 16). A similar situation is shown by the results obtained after applying the Kcross and Lcross functions to the sample in YZ section (Fig. 17), suggesting that the marks tend to separate.

Fig. 16
figure 16

Results obtained from the Kcross and Lcross functions on the pebbles and lithic industry from level 23 of Buena Pinta Cave in section XZ

Full size image
Fig. 17
figure 17

Results obtained from the Kcross and Lcross functions on the pebbles and lithic industry from level 23 of Buena Pinta Cave in section YZ

Full size image

The results obtained after the application of relative risk show a similar result in section XZ and YZ. The pebbles are more likely to be found in the lower half of the level (from − 2.6 m down to the bottom of the sequence), while in the upper half, there is a greater mix of lithic industry and pebbles. Thus, lithics are more likely to be found in certain parts of the upper half and pebbles in the rest of the level (Fig. 18).

Fig. 18
figure 18

Spatially varying relative risk results for the pebble and lithic industry samples separately and together with probability at 99% confidence in section XZ and YZ

Full size image

We also have analysed whether or not there is any kind of correlation between the pebbles, lithic industry, and faunal remains. Again, the results of the χ2 test (Supplementary File 11) indicate that sections XZ and YZ are inhomogeneous samples.

The Kcross and Lcross functions have identified some correlation between SPPs at short distances and negative correlation at long distances in both sections (Supplementary Files 12, 13, 14, and 15). However, the negative correlation between fauna and pebbles starts at even shorter distances than with lithics.

Discussion

General spatial distribution

In the western area of Buena Pinta Cave, during excavation work, a 1.2-m-thick sedimentary unit was identified, which was called level 23. However, the taphonomic study of the assemblage has shown that level 23 is not a homogeneous assemblage, but a sedimentary unit formed from different accumulation processes.

Even though the remains have been affected by different alterations, the faunal remains are uniformly distributed throughout the space. On the other hand, the pebbles and lithic industry tend to be concentrated in the lower half, from − 2.6 m down to the base of the unit. The statistical tests show a tendency for the distribution to be clustered in some parts of the SPP, although when the Nearest-neighbour cleaning test is applied, it identifies most of the surface of both sections as clusters.

When SPPs with marks are analysed to characterize the type of correlation between fauna, pebbles, and lithics, they indicate a negative correlation over long distances, especially fauna and pebbles. In other words, the values of one variable tend to increase, while the values of the other decrease, possibly related to different accumulation processes.

Spatial analyses and graphs showing the vertical distribution of density suggest a positive spatial association between the remains according to the type and degree of post-depositional alteration. Thus, remains with rounding and polishing in stage 2 and 3, or modified by water and concretion, increase in the lower half of the assemblage, from − 2.5 m and below, and concentrate around − 2.8 m, reaching their maximum peak (Table 2, Supplementary Files 6 and 8).

Table 2 Number of remains according to their degree of rounding and polishing by dimensions. In brackets, the % based on the total number of rounded or polished remains of each grade
Full size table

The remains that have been less altered by these post-depositional processes are more abundant above − 2.5 m. The presence of faunal remains with carbonate concretions in a very specific area of the level above − 2.5 m has been interpreted as the result of the fracturing of the fossiliferous breccia (Fig. 2) located in the same area where these remains were recovered. Thus, these bones have been reworked from a level with a different origin and included to the rest of the assemblage. Also, the presence of an accumulation of coprolites just above this shift in depositional dynamics (Fig. 19) supports the hypothesis that these remains have a different origin than those accumulated below.

Fig. 19
figure 19

Spatial projection of the materials forming part of each of the proposed levels. Section XZ (top) and section YZ (bottom)

Full size image

It is worth highlighting the concentration of remains with a higher degree of weathering above this change in dynamics, between − 2.30 and − 2.55 m, which is related to a time of lower sedimentation rate. The presence of these intensely weathered remains directly above a horizontal area with the lowest fossil densities suggests that a sedimentary hiatus has occurred. This would mean a paraconformity exists around the − 2.5 m mark, which would divide the level into an upper and a lower archaeo-palaeontological assemblage, with different taphonomic traits.

In terms of lithics, the results show that above − 2.5 m, there is a greater probability of lithic industry than pebbles to be found. Below this level, there is a greater intensity of pebbles and a lower probability of finding lithic industry. In addition, the K and L functions show the tendency to group both types of lithics (natural and knapped). However, to go further, it will be necessary to improve the lithic studies with the same level of stratigraphic resolution as for the faunal remains.

Furthermore, although statistical analyses indicate that fauna and pebbles tend to be separated, there is a spatial coincidence between the highest pebble intensity and faunal remains affected by post-depositional water-related taphonomic processes (i.e. polishing, rounding, concretion), both in the lower half of the sedimentary package irrespective of the section analysed.

The different patterns of bone distribution based on post-depositional alterations suggest the existence of different agents in the formation of the main accumulations. In addition, data reveals short and long-distance displacements caused by natural post-depositional processes, that could have transported, mixed, and/or accumulated the remains (e.g. flowing water).

The results obtained in this work reveal that the so-called level 23 instead corresponds to three distinct assemblages or levels (Fig. 19) formed by different accumulation processes. Thus, from bottom to top of the stratigraphic sequence, the levels proposed in this work are as follows.

Clay deposit (unit 32A): depth between − 3.2 m and − 2.7 m

It is the oldest unit in the complex, formed by cave clays and separated from the overlying level by an erosive surface. The taphonomic study carried out in this work confirms its existence, despite the small area in which it can be observed.

There is a low density of faunal remains and some have been affected by fossil-diagenetic alterations at an early stage of development, mostly related to humid environments. And 42.25% of the remains are well preserved, while 14.29% have not preserved any of the bone surface. Only 10% of the sample shows evidence of weathering (Table 3), all of which are categorized as being in the mildest stages, and 60.50% of the remains have calcareous concretions.

Table 3 Taphonomic modifications based on the % number of remains observed according to depth
Full size table

Other water-related alterations are the loss of cortical and the whitish appearance of the bones, which are observed in 26.89% of the cases. In addition, 46.22% of the total sample shows rounding, and 28.57% of the remains have been polished. In both cases, there are remains affected from the slightest to the highest degree, although the first stages are the most abundant.

No lithic artefacts have been documented in this level.

Lower level (unit 23): depth between 3 m and 2.5 m

This is a crude fining upwards sequence from boulders and gravels to silts, which includes very rounded and polished (stages 2 and 3) bone remains alongside rounded lithics, composed of very heterogeneous raw materials. Additionally, two elephant tusk fragments are found in this level, which are currently the only documented proboscidean remains in the whole of the Calvero de la Higuera complex.

There is a medium preservation of bone cortical surfaces, although about 30% of the remains are poorly preserved. And 9.06% of the remains are weathered (Table 3). Stages 1–2 predominate despite, with almost 2% of cases showing a higher degree of alteration (stages 3, 4, and 5) located especially at the top of the level. Calcareous concretions have been recorded in almost half of the assemblage (58.99%).

Then, 34.51% of the sample shows alterations related to the physical action of water. 47.72% of the remains have evidence of rounding and 34.23% of polishing. In addition, this is the area of the sedimentary package where more than 72% of the rounded bones in stage 2 and 3 are concentrated, as well as more than 76% of the polished bones in these same stages (Table 2), to the point of losing the original morphology of the bone and having its cortical completely polished.

Lithic remains reach maximum intensity peaks in this area, regardless of whether they are natural or knapped and whether or not they have been eroded. Statistically, it is more likely that pebbles, rather than lithic industry, appear in this level. Furthermore, the lack of correlation shown by the K and L functions together with the spatial projections made with GIS indicate natural clasts are concentrated in this level. There are rounded pebbles and very heterogeneous materials, including allochthonous raw materials transported several kilometres downstream. However, future studies will allow us to explain these data in more detail.

Upper level (unit 2/3): depth between 2.5 m and 1.6 m

The taphonomic analysis indicates a change in the depositional environment, separating this assemblage from the previous one with a paraconformity at around − 2.55 m. There is a low percentage of remains affected by post-depositional taphonomic variables related to water.

Polishing is observed in 9.11% of the remains, whereas rounding in 16.03% (Table 3). In the case of these last two alterations, those in grade 1 predominate (Table 2). Calcareous concretions have been recorded in 24.38% of the remains. Remains with carbonate concretions are only present in this assemblage, towards the base of the level. Biochemical marks left by small micro-organisms and acids present in plant roots are very abundant (77.10%), although this is related to the greater proximity of this package to the surface.

At the base of the level (− 2.5 m), there is an accumulation of coprolites, which indicates at least one carnivore occupation. The anthropic impact is particularly evident given the significant presence of lithic industry, which is more likely to appear in this level.

Taphonomy as a spatial tool

Field layers, typically established during fieldwork based on lithological and archaeological criteria, are normally restricted to what is observable within the excavation areas. The issue arises when these layers are often treated as static and rarely re-evaluated post-excavation. However, Discamps et al. (2023), and references therein, clearly demonstrated the interest of critically reviewing field layers due to the possible existence of features of the archaeological assemblages that are impossible to identify without an in-depth analysis of the remains. Here, we have presented a clear case of a unit identified during the excavation of the western part of Buena Pinta Cave that can be used as an example of this topic.

Post-Excavation Stratigraphy (PES) analysis encompasses a broad range of attributes, such as artifact density, changes in proportions of different types of remains, species distributions, bone surface modifications, lithic implement attribution, raw materials, and refits.

The analysis of palimpsests describing spatial structure by looking at variation in artifact density across a level has been widely applied (see Reeves et al. 2019). However, as we demonstrate in spatial statistics section, there is a potential risk of misinterpreting density measurements as direct indicators of behaviour when they might solely reflect taphonomic or geological processes, or even a combination all of them. To avoid this, detail density and statistical analysis based on taphonomic alterations have been carried out.

In Palaeolithic deposits, where faunal remains are among the most abundant elements of the archaeological record, multiple formation processes can modify the “original” spatial organization and representation of remains (Schiffer 1983). Disciplines such as palaeontology or zooarchaeology can identify and characterize the animal species that appear in the sites. However, only taphonomy provides empirical arguments with which to draw inferences about the processes and agents that modified and accumulated the remains at the study site (Andrews 1995: 147).

As Discamps et al. (2023) proposed and following the original definition of the term (Efremov 1940), we manifest the utility of taphonomy not only for identifying the accumulating and modification agents but for a critical evaluation of field layers and assemblage integrity. This work demonstrates the usefulness and feasibility of integrating taphonomic results as part of the archaeo-stratigraphic analysis, especially when the statistical analyses, as in this case, do not show very enlightening results by themselves, due to the irregular morphology of the excavated area and its extension.

While the significance of spatial statistics in archaeological site interpretation has been acknowledged for a considerable time (Hodder and Orton 1976; Whallon 1974), its current application in the analysis of site formation and modification processes can be considered relatively limited today (Domínguez-Rodrigo et al. 2014; Carrer 2015; Domínguez-Rodrigo and Cobo-Sánchez 2017; Cobo-Sánchez 2020; Diez-Martín F et al. 2021; Moclán et al. 2023a, 2023b; Arteaga-Brieba et al. 2023). Nevertheless, we must take into account the possibility that, as we have shown, due to constraints such as the excavation area, the statistical results may become consistent only when unified with those obtained in the spatial analysis of the assemblage based on taphonomic modifications.

On the other hand, as we have seen, investigations assessing the spatial distribution in three dimensions is necessary to identify distinct patterns and create spatially consistent stratigraphic units. Both horizontal and vertical distributions of bone surface modifications, lithic remains, and coprolites are studied in the case of the western part of Buena Pinta Cave, using two-dimensional plots, as well as variations in density maps. But also, vertical spatial data has been used as a key part of the analysis to define the different layers. This type of study based on the third coordinate (“Z” dimension), although uncommon, can be seen in works like McPherron et al. (2005), Anderson and Burke (2008), Giusti and Arzarello (2016), or Reeves et al. (2019).

In summary, this discussion highlights the necessity of re-evaluating field layers. Spatial analysis and methodologies specifically, such as statistics or density analysis, play a pivotal role in understanding these intricate archaeological contexts and offer a pathway to more robust interpretations of human behaviour in the past. It has been demonstrated how the use of taphonomy is crucial for the correct identification of different events in the processes of formation and transformation of archaeological deposits. But also, combining the analysis of lithic remains and coprolites, we demonstrated that offering an alternative perspective of the same sequence, can contribute information that, when amalgamated, enhances our comprehension of the formation processes and stratigraphy.

Conclusions

The revision of the field layers after the excavation process allows us to highly improve the observations made during archaeological field work. In this sense, this work highlights the importance of re-evaluate field layers after excavation using intra-site spatial data and including taphonomic analysis as a tool for a critical evaluation of field layers and assemblage integrity.

Based on taphonomic analyses, with special attention to fossil-diagenetic processes, we have refuted the integrity of the western part of Buena Pinta Cave and consider that the deposit named as level 23 during fieldwork is, in fact, the product of several accumulation processes resulting in three different archaeological levels. Future research will allow us to interpret in a better way the taphonomic importance of both carnivore and anthropic taphonomic agents at Buena Pinta cave in these three new post-excavation identified levels.