Introduction

Different typologies of sites have been identified in both archaeological and ethnographic records. A general trend can be differentiated between residential and non-residential sites. In the former, daily activities are carried out, whereas the latter are characterized by activities related to obtaining resources (Binford 1978; Scott 1980; Shick 1987; Turq 1988, 2000; Féblot-Augustins 1990; Freeman 1991; O’Connell et al. 1991; Kelly 1992; Baena et al. 2008; Voormolen 2008; Yravedra et al. 2012; Brenet 2013; Bárez del Cueto et al. 2016; Rodríguez-Hidalgo et al. 2017). Sites related to the acquisition of animal carcasses (and where the first butchering phases are carried out) are not common in the archaeological record (Domínguez-Rodrigo 2008), especially if we analyse the most ancient periods of the human evolution. As can be seen in the case of Hadza populations, kill sites are extremely variable, depending on the size of carcasses, the number of hunters, and other issues like the number of butchered animals or the carcass transport and/or the acquisition strategies (i.e. hunting, scavenging) (O’Connell et al. 1988a, b, 1992). As an example, it is not the same if the carcass of an animal is completely butchered at the site where the animal died (i.e. kill site), or if it is transported to a near location where it is subsequently butchered (i.e. near-kill location) (O’Connell et al. 1992).

Thus, two different types of kill/butchering sites can be differentiated: those where the animals died and were butchered (these sites are isolated spots in the landscape and they can be considered as real kill sites), and those where the animals were transported to be initially butchered (i.e. near-kill/butchering sites) selected by their spatial/structural characteristics (e.g. shelters, shaded areas). Both types of sites have been proposed in the archaeological record (Domínguez-Rodrigo 2008; Pope et al. 2020); however, the latter are more probable to be preserved due to the possible recurrent use of these locations (O’Connell et al. 1992; Parfitt and Roberts 1999).

As it is often difficult to determine the real typology of a site, the most general term used to describe these sites is kill/butchering sites (or kill-butchery site) (e.g. O’Connell et al. 1992; Miller and Burgett 2000; Domínguez-Rodrigo 2008; Voormolen 2008; Rodríguez-Hidalgo et al. 2017), which encompasses both types of sites. We use this term to avoid overspeculation about the functionality of occupations and to acknowledge the complex nature of animal exploitation in the archaeological record. It is important to note that in both types of sites, butchering activities are carried out, including the initial consumption of some animal resources (O’Connell et al. 1992; Voormolen 2008; Domínguez-Rodrigo et al. 2015; Rodríguez-Hidalgo et al. 2017). Knapping activities are also commonly observed in these archaeological contexts (Serangeli and Böhner 2012; Serangeli and Conard 2015; Rodríguez-Hidalgo et al. 2017; Pope et al. 2020). In the case of Hadza populations, evidence of fire use has also been identified in some cases (O’Connell et al. 1992).

Following O’Connell et al. (1992), the archaeological characterization of a kill/butchering site is generally carried out using data related to the faunal/lithic assemblage, the internal structure of the site, and the spatial distribution of the archaeological remains. However, it must be noted that in some cases all these features are not taken into account at the same time when interpreting a locus as a kill/butchering site; in those cases, only the taxonomic, taphonomic, or/and lithic information is commonly used (Domínguez-Rodrigo 2008).

The spatial information extracted from Hadza populations has shown that the spatial analysis is especially important in differentiating kill/butchering sites, which are similar to those sites produced by carnivores during hunting. Kill/butchering sites are (in some cases) extremely large due to the repeated episodes of occupation/use of the space, especially when intercept hunting techniques are applied (O’Connell et al. 1992).

For the specific case of the Middle Pleistocene of the Iberian Peninsula, few kill/butchering sites have been excavated and most of them have been identified using only taphonomic/technological features, and the spatial information contained in them has not been sufficiently explored yet. Thus, the level TD10.2 of Gran Dolina (Burgos) (Rodríguez-Hidalgo et al. 2016, 2017) and the levels CB1, CB2, and CB3 of Cuesta de la Bajada (Domínguez-Rodrigo et al. 2015) can be mentioned as clear kill/butchering sites without a published characterization of their spatial patterning. Other elephant-dominated assemblages as Ambrona (Villa and d’Errico 2005; Villa et al. 2005a, b; Sánchez-Romero et al. 2016, 2022), Torralba (Villa and d’Errico 2005; Pineda and Saladié 2019, 2022), PRERESA (Yravedra et al. 2012, 2019; Panera et al. 2014), or Áridos 1 and Áridos 2 (Santonja et al. 1980; Villa 1990) can be mentioned as examples of butchering sites on the Iberian Peninsula. Nevertheless, at these sites, the agency involved in the death of the carcasses cannot be directly related to hominin hunting, and scavenging cannot be rejected, especially in the cases of Ambrona (Villa and d’Errico 2005; Villa et al. 2005a, b) and Torralba (Villa and d’Errico 2005; Pineda and Saladié 2019, 2022).

This lack of spatial research in the Middle Pleistocene kill/butchering sites is not generalized. Recently, an important effort has been carried out in order to investigate the spatial pattern provided by Schöningen 13II-4 archaeological remains (Hutson et al. 2020, 2021; García-Moreno et al. 2021). This site has been interpreted as a kill/butchering site for horses (Equus mosbachensis) which were hunted and accumulated by hominins (Voormolen 2008; Van Kolfschoten et al. 2015a, b; Hutson et al. 2020, 2021; Bonhof and van Kolfschoten 2021) in a main cluster of remains during different occupations. In addition, other taxa were also hunted by hominins, which produced a different spatial distribution of remains (e.g. cervid and bovine remains are segregated). The spatial pattern detected at Schöningen 13II-4 (in which horses appear to be clustered in a specific area) is especially interesting due to their differences/similarities with other types of specialized sites.

If ethnographic information related to kill/butchering sites is considered, it seems likely that Schöningen 13II-4 can relate to some sites created by Hadza when they generate accumulations of carcasses on the same locus during different episodes of carcass acquisition (O’Connell et al. 1992). However, this ‘clustering’ pattern identified with Hadza populations has not been properly (i.e. statistically) analysed, and it is possible that those clustered materials identified by O’Connell et al. (1992) may actually represent a different type of spatial pattern (i.e. random or regular patterns) (Bevan et al. 2013; Sánchez-Romero et al. 2022).

Thus, it is important to do additional spatial analyses of this type of sites in order to create a better referential framework about the kill/butchering spatial pattern variability. This aspect could help in the future to identify the presence of these sites in the archaeological record, especially if the assemblages are not well-preserved, or if the palimpsest nature of the site is very pronounced (Schiffer 1983; Vaquero 2008; Dusseldorp 2009).

In this work, we have investigated the spatial pattern provided by the Middle Pleistocene site of Cuesta de la Bajada, which is characterized as an Early Middle Palaeolithic kill/butchering site of mainly equids (Equus chosaricus), and secondly of cervids (Cervus elaphus) (Domínguez-Rodrigo et al. 2015; Moclán et al. 2022). Thus, this is a site that could be adequately used to be compared with the spatial pattern provided by the lacustrine site of Schöningen 13II-4 (Stahlschmidt et al. 2015; Hutson et al. 2020, 2021; García-Moreno et al. 2021). Thus, our main goal here is to improve the spatial framework of kill/butchering sites in archaeological contexts, especially in relation to the European Middle Pleistocene.

Additionally, the spatial analysis of Cuesta de la Bajada could be useful to identify features related to the biostratinomic and fossildiagenetic (Fernández-López and Fernández-Jalvo 2002) history of the site. The Cuesta de la Bajada archaeological levels are associated with a small deformation depression where a lacustrine environment had been identified with very low-energy and autochthonous sedimentary conditions (in the case of CB3 to CB1 levels [see below]) (Santonja et al. 2014). Weak water alterations have also been found at the site (Santonja et al. 2014; Domínguez-Rodrigo et al. 2015), which were able to modify the original spatial pattern. Thus, this paper has also the aim of spatially characterizing the archaeological features of Cuesta de la Bajada, and the possible degree of modification of the archaeological assemblage underwent during biostratinomy and sedimentation.

Cuesta de la Bajada

Cuesta de la Bajada is located in the Alfambra river valley within the Teruel Depression (Fig. 1a). The site is part of the T4 terrace (+ 50–53 m above the current Alfambra river level) which has been subjected to synsedimentary subsidence (Santonja et al. 2014). Cuesta de la Bajada is located 18 m above the known base of T4, with two fossiliferous areas (i.e. Western Sector and Eastern Sector), which contain deposits with different thicknesses and sedimentary characteristics (Fig. 1b) (Santonja et al. 2000a, b, 2014).

Fig. 1
figure 1

a Location of the study area; (1) Mesozoic and Paleozoic of the Iberian range; (2) Neogene depressions; (3) Faults FC (Concoud) and FT (Teruel); (4) rivers; (5) Cuesta de la Bajada Site; (6) Teruel City. b Stratigraphic section of terrace T4 of the Alfambra river and stratigraphic details of the sequences of the Western (lower right) and Eastern (upper right) Sectors. For a detailed view of the legend of the stratigraphic sections, see previously published profiles published by Santonja et al. (2000a, b, 2014). c Overview of the Eastern Sector during the excavation season of 2011. d Detailed picture of some faunal remains from level CB3 recovered during 2010 excavation

The first archaeological excavations at the site were conducted in the Western Sector between 1990 and 1994. In this area, the different fluvial deposits contain archaeological materials, with up to four fluvial gravel pavements along the sequence (E, G, H, and I) and a succession of fluvial channels and floodplain facies. Levels 12, 16, 17, 18, Pavement I, and 19 provided archaeological materials characterized by being non-autochthonous due to they have been accumulated under high-energy fluvial conditions (Santonja et al. 2000a, b; Sesé et al. 2016).

The second phase of excavations was conducted between 1999 and 2011 in the Eastern Sector in an area of 72 m2. This area shows a series of levels in which the archaeological material is found in an autochthonous position in relation to lacustrine low-energy levels. These levels, which from base to top have been named CB3, CB2, and CB1, have deposits of ascending fine sediments. At the top of this 1.5-m sequence, there is another meter of floodplain facies (level P) (Santonja et al. 2014). In the case of CB1, CB2, and CB3, it must be noted that although the levels have low-energy conditions, water flows could have slightly modified the original position and preservation of the faunal/lithic remains in different areas of the site (Santonja et al. 2014; Domínguez-Rodrigo et al. 2015).

A combination of ESR (Duval et al. 2017) and OSL (Arnold et al. 2016) methods has yielded a date for CB3 between MIS7 and MIS9; the analyses of the herpetofauna have shown that it is possible to date the site in a humid and cold period such as the MIS8 or the MIS9b (Blain et al. 2017). This interpretation is consistent with the dating methods applied to a tuff deposit located at the top of the T4 terrace of the Alfambra river (Simón et al. 2005; Gutiérrez et al. 2008).

The analysis of micromammals (Sesé et al. 2016) and herpetofauna (Blain et al. 2017) of the Eastern Sector has shown the presence of a high number of species, which cannot be related to the anthropic activity (i.e. micromammals seems to be mainly accumulated by raptors and herpetofauna by small carnivorans). From an environmental point of view, the paleoenvironmental reconstruction of the Eastern Sector based on the herpetofauna and pollen analyses suggests a poorly forested, patchy landscape with a large representation of dry meadows, and scrubland habitats, together with well-evidenced aquatic habitats (Blain et al. 2017).

Cuesta de la Bajada is a Middle Pleistocene site at which some of the earliest evidences of Middle Palaeolithic stone tool traditions (Fig. 2) have been documented (Santonja et al. 2000a, b, 2014; Santonja and Pérez-González 2001). Lithic production at the CB3 level (specially in the case of silicified limestones) represents a technology focused on debitage knapped at the site (all phases of the chaîne opératoire are identified), characterized by the presence of ramified production sequences, high percentages of retouched tools and pieces with cortical surfaces, the presence of Levallois reduction strategies, and by the presence of recycling of flakes via resharpening of tools and cores (Santonja et al. 2014). A detailed study of the taphonomic features of the lithic remains and the use-wear analysis are currently under development (P. Bello-Alonso, personal communication).

Fig. 2
figure 2

Examples of lithic remains from levels CB3 and CB2 of Cuesta de la Bajada. CB3: 1–4, scrapers; 5–6, denticulates. CB2: 7, Levallois core; 8–9, scrapers; 10–11, denticulates; 12, awl. Drawings, Raquel Rojas

The zooarchaeological and taphonomic analyses of the levels CB1, CB2, and CB3 have shown that the accumulation of macrofaunal remains is anthropogenic (Domínguez-Rodrigo et al. 2015; Moclán et al. 2022). These levels (mainly CB2 and CB3) have been interpreted as a kill/butchering site of equids (Equus chosaricus) and cervids (Cervus elaphus) which would have been hunted near the site. These carcasses were processed at the site and then butchered remains would probably have been transported to other locations (e.g. base camp).

Materials and methods

Analysed sample

The spatial analysis of Cuesta de la Bajada has been carried out using the archaeological material provided by the Eastern Sector due to its autochthonous nature (Santonja et al. 2014; Domínguez-Rodrigo et al. 2015). During the excavation of this sector, faunal remains ≥ 50 mm were mapped. However, the positions of all lithic industry elements were recorded (Fig. 3).

Fig. 3
figure 3

Spatial distribution of archaeological samples of the levels P, CB1, CB2, and CB3 of the Eastern Sector of Cuesta de la Bajada site. See Supplementary Files 2 and 3 in order to see a detailed map of the distribution of materials of the levels CB2 and CB3

Faunal remains were mapped when they were smaller than the abovementioned measurements if during fieldwork it was considered that they would be of palaeontological/taphonomic interest. However, we have excluded these remains from the analysis to overcome the possible bias produced by the excavators.

Levels P (lithic industry = 49; faunal remains = 45) and CB1 (lithic industry = 307; faunal remains = 202) have provided small, mapped samples, while CB2 (lithic industry = 849; faunal remains = 798) and CB3 (lithic industry = 1927; faunal remains = 1405) show clearly larger samples.

Before starting the spatial analyses, we have carried out a studentized permutation test employing a K inhomogeneous function and a locally scaled K function to test if the samples can be related, or not, with a correlation-stationary or with a locally scaled model (the samples show an inhomogeneous spatial point pattern, see below) (Baddeley et al. 2016). These analyses were conducted using first a tessellation map with 3 divisions (ny = 4) and then with 6 divisions (ny = 4).

These analyses have shown that the level P shows a small sample to test the type of inhomogeneity of the spatial point pattern and, in the case of the level CB1, both correlation-stationary and locally scaled models were rejected. Thus, we have excluded these levels from this study, and we have only analysed CB2 and CB3 (which have been related to a correlation-stationary model) (Table 1).

Table 1 Results provided by the studentized permutation analyses to evaluate the type of inhomogeneity of the samples. Note that only in the case of CB2 the sample was too small to use the first tessellation grid which was used with CB3 and CB1 (x = 3; ny = 4) and we needed to use the second grid (x = 6; ny = 3). None of them was useful to analyse the sample of the level P. C-S, correlation-stationary; L S, locally scaled

It must be noted that we have not included a refitting analysis in this work. While conducting the technological, zooarchaeological, and taphonomic analysis of the samples, we identified the existence of lithic and bone refits. However, we have not yet been able to conduct an exhaustive study of these materials, but we plan to do so soon.

Archaeostratigraphy and spatial modelling

An archaeostratigraphic analysis was conducted to identify the presence of temporal hiatuses in the vertical distribution of CB2 and CB3. These analyses were performed to detect the possible existence of different sublevels of occupation (Fruitet 1991; Canals 1993, 1996; Canals and Galobart 2003; Canals et al. 2003; Rosell et al. 2012; Sánchez-Romero et al. 2017). For this purpose, we used different archaeostratigraphic sections (hereinafter ‘transects’) some ~ 30 cm in width which crossed the site. We have analysed a total of eleven transects per level. These datasets were examined using ArcGIS Desktop v.10.5 software.

However, our archaeostratigraphic approximation was unable to identify different sublevels. As we found it impossible to divide CB2 and CB3 into different sublevels, we devised another method to check whether the material is homogeneously distributed along the vertical section. This aspect is key to understanding the spatial meaning of the accumulations because there could be important spatial differences along the vertical distribution.

First, we divided both levels into two parts (i.e. Upper and Lower), trying to include half of the archaeological material in each of the divisions. Modelling through polynomial cubic regression was performed to create 15 replications of each part of the site (i.e. complete sample and Upper and Lower part of each level) using first all archaeological materials (i.e. fauna and lithic industry) and then differentiating each type of material. The models obtained through spatial modelling were compared with a studentized permutation test (T and \(\overline{T }\) statistics were used). The divisions of CB2 and CB3 were analysed using the methods applied to unmarked patterns (see below) to evaluate through other techniques the similarities/differences between the Upper and Lower parts of both levels.

Spatial analysis

Methods

The spatial analysis of the site was performed using the ‘spatstat’ (Baddeley and Turner 2005) library in R (R Core Team 2022). RStudio (R Studio Team 2021) was used to write the code, which is included as a supplementary document (see Supplementary File [SF] 1). This library allows us to analyse spatial point patterns (SPP) represented in 2D to evaluate their homogeneity, the intensity of the patterns, if the samples are randomly distributed, or not, and the relationships between different types of remains, among other issues. For a complete explanation of the methods included in ‘spatstat’, see Baddeley et al. (2016).

Before explaining the applied methods, it is important to note the difference between a SPP with ‘marks’ and those without ‘marks’. An unmarked SPP is that defined as containing information about only one type of element (e.g. faunal remains, lithic industry, retouched tools) while SPPs with marks are those that consider different types of elements at the same time (e.g. comparison between the SPPs of faunal remains and lithic industry). This issue is important because the type of statistical analysis differs depending on the type of SPP (i.e. unmarked or marked/multitype).

In both cases, the first step was to test the homogeneity of the SPP through the χ2 test. In this case, we have performed the χ2 test using four different grids to do an accurate analysis of the samples (see SF 1 for more details). As most of the analysed SPPs showed an inhomogeneous distribution, we considered using the inhomogeneous versions of the functions (see below). Then, we performed a first map which shows the intensity of the SPP using the Kernel density test with likelihood cross-validation correction.

Then, in the case of unmarked SPPs, we performed a second map which shows the presence of hotspots by calculating the likelihood ratio test. This map was also calculated using only those areas where the presence of archaeological material was 99% significant. For marked patterns, spatially varying probabilities were determined for each mark. When the number of marks was > 2, a joint assessment was made of the probability of certain marks appearing within the same space as another.

For unmarked patterns, it is key to understand if complete spatial randomness (CSR) can be rejected or not. A CSR model is the null hypothesis of some analyses that evaluate what is the type of distribution of each SPP. If CSR is rejected, a clustering (materials share a positive relationship among them) or a regular (the spatial relation among materials is negative) pattern can be proposed depending on the results of the tests.

First, Hopkins-Skellam and Clark-Evans tests were performed; however, these tests are extremely sensitive to inhomogeneous samples and commonly return a false ‘cluster’ result. However, in the case of small samples (e.g. tooth marks) whose homogeneity could not be rejected through the χ2 test, these tests have been useful. When inhomogeneous distributions were identified, these analyses can be used to reaffirm the results provided by the χ2 test.

We have also used the pair correlation function (only with those samples of the Upper and Lower divisions of the levels) and the K, L, F, G, and J functions to test the null hypothesis of CSR with an inhomogeneous Poisson process. In addition, we have used another version of the K function (hereinafter ‘modified’ K function) (Cobo-Sánchez 2020); the reason for using this version of the K function is that the ‘normal’ version uses a standard inhomogeneous Poisson process as the null hypothesis, while the ‘modified’ version uses an inhomogeneous Poisson process with the same intensity than that provided by the analysed spatial pattern. With this version of the function, we can overcome the possible error of characterizing a high-intensity variation with positive correlation of the points (i.e. clustering).

Finally, we have also used the maximum absolute deviation (MAD) and Diggle-Cressie-Loosmore-Ford (DCLF) tests, which can identify whether a SPP corresponds to the CSR model (in this case H1 = CSR rejected, but cluster or regular patterns are not proposed). MAD and DCLF tests have been performed using the K, L, F, G, and J functions.

In the case of marked SPPs, the main objective is to identify the possible codependency between marks. In this sense, the rejection of CSR is not as key as identifying the possible relation between different archaeological elements (i.e. codependency). Thus, we have used the segregation test and the nearest-neighbour contingency table analysis to do a first exploratory analysis of the relationship between several types of materials. Nearest-neighbour contingency table analysis (NNCTA) was performed using the R library ‘dixon’ (de la Cruz 2008). However, these tests are often not very robust (Baddeley et al. 2016) and it is preferred to use other methods to test the codependence between marks. In this sense, we have performed both Kcross and Lcross functions for inhomogeneous samples which are capable to examine the presence of codependency between marks by pairs in a more accurate manner.

Note that CSR is used in the case of homogeneous distributions (Baddeley et al. 2016) while in the case of inhomogeneous distributions the null-hypothesis is referred as ‘inhomogeneous Poisson distribution’. However, in order to be clear during the exposition of the results, we will use the concept of CSR as synonym of inhomogeneous Poisson distributions.

Selected variables

In the case of the faunal remains, we have analysed as unmarked SPPs the distribution of all faunal remains, specimens identified as large- (mainly horses) and medium-sized animals (mainly cervids), the presence of manganese oxidations, biochemical marks, rounded and polished bones, trampling marks, green-fractured bones, cut and percussion marks, and carnivore tooth marks. We have analysed faunal remains as marked SPP with the following variables: maximum length, degree of manganese-oxidated surface, the degree of manganese accumulated on remains (i.e. isolated [low], clustered [low-medium], widespread [medium–high], and complete [all the surface is affected by manganese]) (Cáceres 2002), cortical preservation (i.e. using different degrees of preservation [0 = 0%; 1 = 1–25%; 2 = 25–50%; 3 = 50–75%; 4 = 75–99%; 5 = 100%] or dividing the sample between ‘good-preserved specimens’ [degrees 3 to 5] and ‘bad preserved specimens’ [degrees 0 to 2]) (Moclán et al. 2021), rounded and polished bones using different degrees of preservation (i.e. 0 = non-modified; 1 = low damage; 2 = medium damage; 3 = high damage) (Cáceres 2002).

For the lithic industry, we have analysed the complete sample, pieces over 9 mm and pieces over 19 mm. Then, we analysed other variables, but we have always used only those materials over 19 mm to overcome the possible bias produced by the excavation process (i.e. smaller remains are not always identified during the excavation process and they are found during subsequent dry-sieving of the sediment) (Santonja et al. 2014). We have analysed as unmarked SPPs the distribution of each raw material (i.e. silicified limestone, ‘normal’ limestones, quartzites, quartz, and flint) and then we have tried to analyse the distribution of flakes, retouched tools, and cores per raw material (note that some samples are too small and they could not be analysed, e.g. flint flakes or retouched tools). We have also analysed as marked SPPs the distribution of different types of tools with cortical surfaces (raw materials have also been taken into account in these analyses). In addition, we have compared as marked SPPs the distribution of lithic pieces using the information related to the maximum length of the items and the raw materials.

As in the case of the faunal remains, we have analysed the rounded artifacts. First, rounded tools have been analysed as an unmarked pattern and then we have used different degrees of rounding (R0 = non-modified; R1 = low damage; R2 = medium damage; R3 = high damage).

Finally, we have compared the distribution of faunal and lithic remains in order to test the codependency between them (in this case, we have used only those lithic pieces over 49 mm to overcome the bias produced during the mapping of faunal remains). This comparison has been carried out using the complete samples and using the maximum length of the samples (smaller remains = 50–69 mm; larger remains =  > 70 mm). We have also compared those water-affected (rounded and polished items) remains and those that were interpreted as percussion tools with those bones with percussion marks.

In Table 2, a summary of the variables analysed in this work can be visualized.

Table 2 Different archaeological variables statistically analysed through the ‘spatstat’ library of R software

Some of these variables have been analysed to identify patterns of anthropic activity at the site, such as the distribution of faunal and lithic tools, the distribution of each raw material or the distribution of anthropic marks, and the location of the large- and medium-sized animals. Other variables have been analysed to test the presence of bioestratinomic and diagenetic modifications of the assemblage after anthropic abandonment of the carcasses and lithic tools (e.g. tooth marks, biochemical marks, trampling, manganese oxidations, rounding). However, other variables can be used to evaluate both types of issues, such as the relation between faunal and lithic remains (water-affected or not), even when they are analysed considering the specimen length.

It must be noted that when the samples were especially small in order to carry out the statistical analyses, we have performed only those that could be done. In those cases, the results must be interpreted with caution (this situation is typical when marked SPPs are compared and one or more marks are unbalanced in comparison with others).

Results

Archaeostratigraphic analysis

The archaeostratigraphic analysis of the level CB2 is completely shown in Supplementary File 2. The visualization of the transects has shown that there are no clear hiatuses in the level that could be interpreted as different archaeological sublevels of the level CB2 (Fig. 4). Although there are some areas along the surface of the level where there seem to be hiatuses at the middle part of the section (e.g. xy view of transect 8, SF 2 p. 9), these empty spaces are not extended along the surface of the level and we have rejected the possibility of differentiating archaeological sub-levels.

Fig. 4
figure 4

Archaeostratigraphic transects 3 and 10 of CB2 and CB3 (view xz axes). Note that sublevels cannot be identified, and the slight slope is oriented to NW. Also note that the scale of each transect is not the same (axis values are shown in meters)

The most interesting aspect related to the archaeostratigraphic analysis is that all materials are clustered along 30 cm depth. It must be noted that this situation is clearer in the west/northwest part of the site, while in the east/southeast part, the accumulation of remains along the vertical section is clustered in a thinner area. This aspect is interesting because CB2 is 50 cm deep in some areas and its sedimentological characteristics are different (moderate-higher energy) at the top of the level. Clustering of materials along only 30 cm (or even less in some cases) can be considered as a sign of good preservation of several of the original spatial features.

Another characteristic of the assemblage that can be inferred through the visualization of the transects is that the slope is slight along the surface, although the slope of the level is oriented (i.e. northwest) to the current fluvial channel of the Alfambra river. This situation is expected for a floodplain like that represented by Cuesta de la Bajada.

The archaestratigraphic results of CB3 are shown in Supplementary File 3. This level has provided comparable results to those provided by CB2.

First, all materials are clustered around ~ 30–40 cm deep, especially in the case of the northwest part of the site. However, in this case, this observation cannot be considered meaningful because CB3 is vertically thinner than CB2. The observations about the slope of the level and the absence of clear hiatuses are also clear at CB3.

Spatial analysis

All spatial analyses carried out using ‘spatstat’ are shown as supplementary files (SF 4, SF 5 and SF 6). Below we have introduced some figures as general guidelines, but for a complete interpretation of the results please refer to the supplementary documents.

Testing vertical distribution

The spatial modelling of archaeological materials from CB2 and CB3 and the Upper and Lower half of the levels through the studentized permutation test have shown that there are clear similarities along the vertical section of both levels (Table 3).

Table 3 Results provided by the studentized permutation test applied to the simulations created with the spatial point pattern of the levels CB2 and CB3

If the complete archaeological sample of CB2 is tested (i.e. faunal remains and lithic industry at the same time), significant values (< 0.05) are obtained when all materials are compared with the Upper and Lower parts of the level and also when the Upper and the Lower parts are compared between them. This result is also obtained when faunal remains are considered. In the case of the lithic industry, significative values are obtained when the complete level is compared with the Upper part and when the Upper part and the Lower part are compared using the \(\overline{T }\) statistic. However, when all lithic remains are compared with the Lower part of CB2 and when the Upper and the Lower parts are compared using the T statistics, non-significative results are obtained.

In the case of CB3, similar results are obtained but few differences are shown. The comparison of the Upper and the Lower parts of the levels has always shown significant values (< 0.05). However, when all materials of the level are compared with the Upper or the Lower parts of the level, some values are non-significative (i.e. complete sample with T and \(\overline{T }\) statistics, lithic industry with T statistic, faunal remains with T and \(\overline{T }\) statistics).

The spatial analysis of the SPP of CB2 and its Upper and Lower divisions are shown at the SF 4 (pp. 1–6). The results obtained are extremely similar when all materials are compared with the distribution of the Upper and the Lower parts.

Hopkins-Skellam and Clark-Evans tests have always shown a cluster distribution of the samples. Pair correlation functions have shown a weak cluster model in short distances (except for the lithic industry of the Upper part of the level which cannot reject the CSR model). The application of G, F, and J functions has provided inconclusive results because they could not reject the CSR model; however, the K and L functions show a model with a weak cluster at short distances and a regular pattern at long distances (K function applied to all distribution of faunal remains and those lithic tools related to the Upper part of the levels shows only a regular pattern). The application of the ‘modified’ K function suggests that the pattern is regular for all analysed samples.

The analysis of CB3 has shown similar results (SF 4, pp. 7–12). Hopkins-Skellam and Clark-Evans tests showed a cluster pattern for all samples. The pair correlation function shows again a cluster pattern at short distances when both faunal and lithic remains are analysed at the same time and when the lithic industry is analysed (except for the Lower part of the level which cannot reject the null hypothesis). In the case of the faunal remains, the null-hypothesis of inhomogeneous CSR could not be rejected for the complete sample of CB3 and the Upper part of the level, while it was rejected for the Lower part. F, G, and J functions did not reject the CSR model while the L function showed a model characterized by a slight cluster distribution of faunal remains at short distances and regular at longer ones (the K function and the ‘modified’ version of the K function shows regular SPPs).

The combination of the spatial modelling analysis and the unmarked spatial analysis of both levels has shown that there are clear similarities between the Upper and the Lower parts of the levels, which cannot be differentiated from a statistical point of view. In this sense, we consider that the spatial analysis of both levels can be conducted using all mapped materials of each level without the need to split the levels into archaeological sub-levels.

Spatial analysis of CB2 level

All spatial results related to the level CB2 can be seen in Supplementary File 5.

The analysis of all faunal remains (Fig. 5a; SF 5, p. 1) and all lithic tools (complete sample, > 9 mm and > 19 mm) (Fig. 6a–c; SF 5, pp. 19–21) show a similar spatial distribution. The intensity and hotspot maps show that there are no clear clusters but rather the spatial pattern is characterized by significant accumulations of most of the materials on the central/central-western part of the excavated area. The application of the functions has rejected the CSR model when K, L, and ‘modified’ K functions are used (Fig. 7). All these functions identify the SPP as regular (K and L also identify a slight cluster pattern at short distances in the case of lithic tools; in the case of faunal remains only the L function identifies a slight cluster pattern at short distances).

Fig. 5
figure 5

Maps showing the kernel smoothed intensity and the hotspots (‘TRUE’ areas) with a 99% confidence interval for the different faunal remains examined from level CB2. a) Total sample of faunal remains; bc) faunal remains considering animal body size; dl) different taphonomic alterations of the faunal remains

Fig. 6
figure 6

Maps showing the kernel smoothed intensity and the hotspots (‘TRUE’ areas) with a 99% confidence interval for the different lithic remains examined from level CB2. ac) Total sample of lithic remains; dg) Silicified limestone; hk) ‘normal’ limestone; lo) quartzite; p) quartz; q) flint; r) rounded remains

Fig. 7
figure 7

Results obtained by the inhomogeneous versions of the K, L, and the ‘modified’ version of the K functions from level CB2

These results are also supported (cluster pattern identified by K and L function is not always identified) by the analysis of other faunal-related variables as the presence of large- (Fig. 5b; SF 5, p. 3) and medium-sized (Fig. 5c; SF 5, p. 4) animals, manganese oxidations (Fig. 5d; SF 5, p. 5), biochemical marks (Fig. 5e; SF 5, p. 9), rounded bones (Fig. 5f; SF 5, p. 10), and green-fractured bones (Fig. 5i; SF 5, p. 15). However, the null hypothesis of CSR cannot be rejected if the spatial distribution of polished bones (Fig. 5g; SF 5, p. 13), cut-marked bones (Fig. 5k; SF 5, p. 17), and percussion-marked bones (Fig. 5l; SF 5, p. 18) are analysed. Tooth marks (Fig. 5j; SF 5, p. 16) show a different pattern because they show a homogeneous distribution and Hopkins-Skellam and Clark-Evans test did not reject the CSR model; the ‘modified’ K function has identified a slight cluster pattern (although it should be noted that the sample is very small).

When considering the marked SPPs of faunal remains, relationships among materials are identified. If the analysis is performed using the length of the specimens, it can be noted that those remains of 60–69 mm of maximum length are clustered with all the other analysed groups if Kcross and Lcross functions are used (NNCTA shows overall segregation among groups). There is a positive correlation between bones with good and bad cortical preservation (SF 5, p. 8), although the analysis cannot be performed (i.e. Kcross and Lcross functions) with different degrees of cortical preservation due to the small sample size (SF 5, p. 7). No codependency (or very weak) is detected if different degrees of manganese oxidations or rounding are analysed.

If the distribution of lithic industry per length groups (SF 5, p. 22) is analysed, the sample size does not allow us to perform the Kcross and Lcross functions. However, the results provided by the segregation test and the NNCTA are conflicting due to the first one identifies segregation among groups and the second one could not identify segregation. Nevertheless, if the analysis is carried out using the different raw materials (SF 5, p. 23), cluster-type codependency is detected among all raw materials except for flint and quartz (quartz shows a positive relationship only with limestone).

Silicified limestone shows a regular pattern if ‘modified’ K, K, and L functions are used (Fig. 6d; SF 5, p. 24), although the latter also shows a cluster pattern at short distances (as in the case of all lithic remains). However, it must be noted that the hotspot analysis clearly shows the presence of three high-intensity areas (see Fig. 6d–g), which are related to the presence of flakes (Fig. 6e; SF 5, p. 25), retouched tools (Fig. 6f; SF 5, p. 26), and cores (Fig. 6g; SF 5, p. 27). These significative hotspots however are not considered as clusters if functions are used, but rather the ‘modified’ K function is also showing a regular pattern. The analysis of the cortical tools (SF 5, p. 28) has shown positive correlation between cores and chunks, cores and flakes, and cores and retouched tools when the Kcross function is applied.

The SPP of ‘typical’ limestones shows a completely different spatial distribution (Fig. 6h; SF 5, p. 29). If all remains are analysed at the same time, contradictory results are provided by K and L functions, which provide a regular pattern interpretation, and the ‘modified’ K function, which shows a cluster pattern (this cluster pattern seems to be consistent with the hotspots analysis which has shown a large high-intensity area in the central part of the site). However, the analysis of flakes (Fig. 6i; SF 5, p. 30), retouched tools (Fig. 6j; SF 5, p. 31), and cores (Fig. 6k; SF 5, p. 32) have not provided significant results (i.e. CSR cannot be rejected). Cortical surfaces cannot be analysed due to the small sample size with cross function and NNCTA (SF 5, p. 33).

Quartzite elements also show one high-intensity area at the northwest area of the site (this area is similar to that provided by ‘normal’ limestones) (Fig. 6l; SF 5, p. 34), which can be related to the presence of all types of tools (i.e. flakes [Fig. 6m; SF 5, p. 35], retouched tools [Fig. 6n; SF 5, p. 36], and cores [Fig. 6o; SF 5, p. 37]). However, a regular pattern is proposed for all quartzite elements, flakes, and cores and the CSR model cannot be rejected for retouched tools. No clustering or segregation of cortical quartzite elements is detected (SF 5, p. 38).

In the case of quartz (Fig. 6p; SF 5, p. 39) and flint (Fig. 6q; SF 5, p. 40) tools, CSR could not be rejected by any function.

Rounded lithic elements show a regular pattern if the ‘modified’ K function is considered, and clustering at short distances and regular at long distances if K and L functions are used (Fig. 6r; SF 5, p. 41). Codependency among different degrees of rounding is not identified by cross functions if non-rounded elements are excluded from the analysis (SF 5, p. 42), whereas if these latter are included strong positive correlation is identified among slight (R1) and medium-rounded (R2) lithic elements (SF 5, p. 43).

If a comparison between the distribution of faunal and lithic remains is performed, strong positive correlation is identified between both types of materials (Fig. 8a; SF 5, p. 44). This type of codependency is also identified if the length of the materials is considered (Fig. 8b; SF 5, p. 45). In this case, small lithic remains (50–59 mm) show positive correlation with both large and small faunal remains while large lithic remains seem to be not related to the distribution of faunal remains. Small and large faunal remains show weak regular codependency between them.

Fig. 8
figure 8

Results for the Kcross and Lcross function analyses of the faunal and lithic industry remains (a, b), percussion tools and faunal remains with percussion marks (c), and water-affected remains (d) from level the CB2

If percussion tools and percussion-marked bones are analysed, the cross functions do not identify correlation between lithics and faunal remains (Fig. 8c; SF 5, p. 46). This situation is also identified in the case of the comparison between water-affected bones and lithics (Fig. 8d; SF 5, p. 47). However, for both cases, NNCTA shows the possible segregation between materials (in the case of percussion/percussive remains the segregation test also shows the possible existence of segregation between materials).

Spatial analysis of CB3 level

All spatial results provided by the level CB3 can be seen in Supplementary File 6.

As in the case of CB2, at the CB3 level, the faunal (Fig. 9a; SF 6, p. 1) and lithic (all, > 9 mm and > 19 mm) (Fig. 10a–c; SF 6, pp. 21–23) remains are mainly distributed at the north-western half of the site, and a regular SPP can be proposed if the K and the ‘modified’ K function are used to identify the typology of the SPP. As in the case of CB2, if the L function is performed, a slight cluster pattern is detected at short distances and a regular pattern at long distances (Fig. 11).

Fig. 9
figure 9

Maps showing the kernel smoothed intensity and the hotspots (‘TRUE’ areas) with a 99% confidence interval for the different faunal remains examined from level CB3. a) Total sample of faunal remains; bc) faunal remains considering animal body size; dl) different taphonomic alterations of the faunal remains

Fig. 10
figure 10

Maps showing the kernel smoothed intensity and the hotspots (‘TRUE’ areas) with a 99% confidence interval for the different lithic remains examined from level CB3. ac) Total sample of lithic remains; dg) Silicified limestone; hk) ‘normal’ limestone; lo) quartzite; p) quartz; q) flint; r) rounded remains

Fig. 11
figure 11

Results obtained by the inhomogeneous versions of the K, L, and the ‘modified’ version of the K functions from level CB2

This pattern is also identified in the case of other faunal variables as the presence of large- (Fig. 9b; SF 6, p. 3) and medium-sized animals (Fig. 9c; SF 6, p. 4), manganese oxidations (Fig. 9d; SF 6, p. 5), biochemical marks (Fig. 9e; SF 6, p. 9), trampling marks (in this case K and L functions show only the existence of a regular pattern) (Fig. 9h; SF 6, p. 16), and green-fractured bones (Fig. 9i; SF 6, p. 17).

In the case of the rounded (Fig. 9f; SF 6, p. 10) and polished (Fig. 9g; SF 6, p. 11) bones, the ‘modified’ K, K, and L functions identify a regular pattern (for rounded bones the L function also identifies a cluster pattern at short distances). However, the visualization of intensity and hotspot maps show that in the north part of the western area of the site (near the north-western limit of the excavation area) high-intensity areas are clear with water-affected bones. In this sense, the water-affected faunal remains are differently distributed than non-affected specimens.

If rounded bones are analysed using the different degrees of this modification (SF 6, pp. 11–12), the absence of codependency between degrees is detected (NNCTA identify segregation between degrees when non-rounded bones are excluded [SF 6, p. 11] but it rejects this possibility when non-rounded specimens are included in the analysis [SF 6, p. 12]). The absence of codependency among degrees of polishing is detected too (SF 6, pp. 14–15).

Anthropogenic and carnivore marks appear to be distributed in a different way than most of the faunal remains. In the case of tooth marks (Fig. 9j; SF 6, p. 10), it is not clear if the sample shows a homogeneous or an inhomogeneous distribution (χ2 tests contradictory results depending on the analysed grid). CSR could not be rejected (except for the Hopkins-Skellam test). However, it must be noted that a high-intensity area is detected by hotspot analysis that is different to those detected by the analyses of anthropogenic marks. The analysis of cut marks (Fig. 9k; SF 6, p. 19) could not reject the presence of a CSR model, which can be rejected if percussion marks (Fig. 9l; SF 6, p. 20) are analysed with the K function and the DCLF and MAD tests with the L function. These analyses are probably weak due to the small samples, but what seems to be clear is that a different high-intensity area is identified is both types of marks (i.e. anthropogenic and carnivore marks) are compared.

If faunal remains are analysed using the maximum length of the specimens, slight regular codependency seems to be identified by the Lcross function between few length intervals. This result is supported by the segregation test and the NNCTA (SF 6, p. 2).

The analysis of different degrees of manganese oxidations and the cortical preservation also shows slight relations of codependency among materials which are not supported by all statistical analyses (SF 6, p. 6–7). In the case of the manganese oxidations (SF 6, p. 6), clustered and isolated-oxidated bones are not related between them, but both appear to be clustered with widespread-oxidated specimens if the Kcross function is applied (however this analysis is not supported by the segregation test and the Lcross function). In the case of the cortical surfaces, a slight clustering relation is proposed by Kcross and Lcross among degree 4 and the other degrees of preservation (SF 6, p. 7); however, if bad and good degrees of preservations are compared, slight regular codependency is identified by Kcross and Lcross functions (although segregation tests and NNCTA reject this relation between materials).

As previously noted, the lithic industry shows a regular pattern if ‘modified’ K, K, and L functions are used to characterize the SPP (or even a cluster pattern at short distances if K and L are considered) (Fig. 10a–c; SF 6, p. 21–23). However, if specimen length is considered, positive correlation is identified among 60–69-mm-length materials and other length groups (SF 6, p. 24). Nevertheless, in the case of CB3, the absence of codependency is detected among diverse types of raw materials (SF 6, p. 25) if cross functions are performed (although the segregation test and NNCTA seem to identify segregation among groups).

Cortical analyses have provided weak results due to the small sample sizes obtained for some types. In this sense, we have discarded these results which have only been done as exploratory. However, we will expose differences between types about their distribution when they show results clearly identified de visu.

The analyses conducted with all silicified limestone remains and with flakes show SPPs characterized as regular ones by ‘modified’ K, K, and L functions (the latter also identify a slight cluster pattern at short distances in the first case) (Fig. 10d–g; SF 6, pp. 26–29). Retouched tools are characterized by the presence of a regular SPP in the case of the K and L functions (‘modified’ K function has not rejected the CSR model) (Fig. 10f; SF 6, p. 28). All these materials seem to be distributed along the western area of the site without a clear high-intensity area; however, in the specific case of the cores, they actually seem to be ‘clustered’ (note that the pattern has been characterized as regular) in the north-eastern part of the western area of the site (Fig. 10g; SF 6, p. 29). In this sense, a difference appears at least in the case of cores and all other types. Although the sample size is small to consider cortical elements from a statistical point of view (SF 6, p. 30), it must be noted that high-intensity areas of cortical cores, chunks, and flakes are located in the same area, which can be an indicator of codependency among them.

The sample of ‘normal’ limestones shows a different pattern (Fig. 10h–k; SF 6, pp. 31–34) than which that we have been identified for silicified limestones. If all pieces are analysed at the same time (Fig. 10h; SF 6, p. 31), a regular pattern is identified by the ‘modified’ K, K, and L functions, which is characterized by three high-intensity areas identified by the hotspot analysis. The western high-intensity area is clearly related to the presence of simple flakes (Fig. 10i; SF 6, p. 32), which also appear to be distributed in a regular way. In relation to the retouched tools (Fig. 10j; SF 6, p. 33), the sample is too small (n = 9), and the results can be discarded. However, cores (Fig. 10k; SF 6, p. 34) seem to be related to the eastern high-intensity area; in this case, only the ‘modified’ K function has identified a regular pattern which probably could be rejected due to the results provided by Hopkins-Skellam, Clark-Evans, and χ2 tests. The distribution of cortical elements seems to be non-related in this case (SF 6, p. 35).

Quartzite pieces (Fig. 10l–o; SF 6, pp. 36–38) show a regular distribution (K, L, and ‘modified’ K functions) and SPP characterized by the presence of two high-intensity areas. These two high-intensity areas are mainly composed of simple flakes (Fig. 10m; SF 6, p. 37) (which also are characterized by the presence of a regular pattern); retouched flakes (Fig. 10n; SF 6, p. 38) are also included in these high-intensity areas, but the functions could not reject the CSR model. In the case of cortical elements, it seems that a relationship between cortical cores and cortical flakes (SF 6, p. 39) is detected due to their higher intensity areas being located in the same spaces of the site.

Both quartz and flint elements have small sample sizes but each one is showing a completely different model to those previously mentioned for other raw materials. In the case of quartz (Fig. 10p; SF 6, p. 40), the CSR model cannot be rejected but there is only one high-intensity area identified by hotspot analyses. However, in the case of flint elements (Fig. 10q; SF 6, p. 41), the spatial pattern does not show a clear high-intensity area and the materials are ‘homogeneously’ distributed for the western half of the excavation area (in this case only the K function have identified a slight regular pattern).

Rounded lithic tools (Fig. 10r; SF 6, p. 42) show a regular-distributed SPP if the ‘modified’ K function is performed, while clustering at short distances is detected by K and L functions (a regular SPP is also identified at long distances). The high-intensity areas are extremely similar to those identified by the complete lithic sample (see Fig. 10a–c). If different degrees of rounding are considered (SF 6, p. 42), the Lcross function and the NNCTA identify slight positive correlation between R1 and R2 degrees. Lcross function has also identified this relationship if non-rounded pieces are included in the analysis, a negative relationship between non-rounded pieces and R1 pieces is also detected (this last relationship is also detected by the Kcross function) (SF 6, p. 44).

When faunal and lithic remains are compared, a slight cluster relationship at short distances is detected if the Lcross function is performed (Fig. 12a; SF 6, p. 45). This relationship is confirmed by the comparison of lithic and faunal remains considering the length of the specimens (Fig. 12b; SF 6, p. 46). In this second case, a clear positive correlation is identified between small faunal and lithic remains, while a regular relationship is identified if large and small faunal remains are compared.

Fig. 12
figure 12

Results for the Kcross and Lcross function analyses of the faunal and lithic industry remains (a, b), percussion tools and faunal remains with percussion marks (c), and water-affected remains (d) from level the CB3

The same type of relationship (positive relationship [i.e. clustering]) is identified if percussion tools and bones with percussion marks are compared (Fig. 12c; SF 6, p. 47). However, if water-affected elements are compared (i.e. rounded and polished materials) (Fig. 12d; SF 6, p. 48), positive correlation is identified (only Lcross function) at short distances and a regular relationship is detected at long distances.

Discussion

The general SPP of Cuesta de la Bajada

Although the spatial characteristics of kill/butchering sites have been proposed as a key element to take into account in order to identify this type of site (O’Connell et al. 1992), little archaeological information is available to help us to understand the SPPs of these sites. Here we have shown a complete spatial analysis of two archaeological levels of Cuesta de la Bajada site that could be useful in order to understand the variability of the SPPs of kill/butchering sites.

Previous technological and taphonomic analyses of the levels CB2 and CB3 (Santonja et al. 2014; Domínguez-Rodrigo et al. 2015; Moclán et al. 2022) of Cuesta de la Bajada pointed out a clear relation of these levels with the presence of a kill/butchering site, with butchering and knapping activities carried out in situ. However, the same studies (Santonja et al. 2014; Domínguez-Rodrigo et al. 2015; Moclán et al. 2022) were cautious about the possible similarities between CB2 and CB3 due to the smaller sample size of the level CB2.

Our spatial analyses have shown that both levels show SPPs with clear regular distributions. This result allows us to propose that the spatial pattern produced by hominin activity at both levels was homologous from a general point of view.

This phenomenon is especially interesting due to the previous assumption of other authors about the possible existence of clustered patterns in kill/butchering sites identified in the ethnographical (O’Connell et al. 1992) and in the archaeological (Hutson et al. 2020; García-Moreno et al. 2021) records. In the Middle Pleistocene site of Schöningen 13II-4, equids appear to be clustered, while cervids and bovines show a cluster distribution only at short distances and regular at long ones (Hutson et al. 2020; García-Moreno et al. 2021).

In this sense, Cuesta de la Bajada is showing a spatial distribution of archaeological materials (both lithics and fauna) that could be defined as unexpected. Thus, two different kill/butchering sites of the same period (and with technological similarities) (Serangeli and Conard 2015) are clearly showing different spatial distributions that are mutually contradictory.

It must be noted that the application of K and L functions has shown in some cases a slight cluster model at short distances in both Cuesta de la Bajada levels as in the case of cervids and bovines of Schöningen 13II-4 (Hutson et al. 2020; García-Moreno et al. 2021). However, the ‘modified’ K function seems to be more accurate to determine the real SPP of the site. This statement is based on that both levels show strong differences of intensity in different parts of the studied area (i.e. southeast part of the site shows a clear lower intensity in comparison to the north-western part), and this fact can create weak or incorrect results when ‘normal’ functions are applied (Cobo-Sánchez 2020). In this sense, we think that it could be interesting to apply the ‘modified’ version of the K function with other previously analysed samples due to it improves the knowledge about the spatial pattering of those sites. In any case, as previously said, the SPPs of faunal remains and lithic tools of Cuesta de la Bajada could be considered as regular.

However, high-intensity areas have been identified in some cases, especially when some lithic samples are analysed. In this sense, a question can be done from an archaeological point of view: is the existence of high-intensity areas contradictory to the regular pattern of the materials? This question is important due to some of our results have shown that the areas occupied by different materials seem to be clearly distinguishable (see for example the spatial distribution of silicified limestone tool types on both levels). However, the existence of high-intensity areas is not contradictory with a general regular trend at the site. This type of area can be identified irrespective of whether the SPP is characterized as regular, cluster, or even CSR (Baddeley et al. 2016).

In any case, it seems to be key to know if the regular pattern is identifiable at the site. At the ‘Results’ section, we have used spatial methods related to the use of distances among specimens to do the spatial analyses. If we performed a new approximation to the SPP of the site using the neighbours of the points instead of their distances, we probably can identify if the regular pattern has been properly identified. Thus, we have carried out another test (nearest-neighbour cleaning) with the function ‘nnclean’ of ‘spatstat’ (Baddeley et al. 2016) library of R (R Core Team 2022) to identify a possible clustered area and those points that are not included in that possible cluster (the tests have been performed using faunal and lithic remains and k values of 10, 20, and 50). Here, this approximation is useful due to it assumes that a cluster area is present in the assemblage (Baddeley et al. 2016).

As can be seen in Fig. 13, most part of the assemblage has been identified as ‘clustered’. This situation has been interpreted by us as new evidence of the absence of a clear clustered SPP at Cuesta de la Bajada site.

Fig. 13
figure 13

Nearest-neighbour cleaning for CB2 and CB3 samples of faunal remains (≥ 50 mm) and lithic industry (≥ 20 mm). Note that 10, 20, and 50 values have been used when k was modified

This main result provided by Cuesta de la Bajadas is especially interesting because its distribution is clearly different not only from that of the equids of Schöningen 13II-4 (Hutson et al. 2020; García-Moreno et al. 2021) but also from other sites interpreted as butchering sites, such as ST4 site at Peninj (Domínguez-Rodrigo et al. 2007). Here, we can mention the case of the Sivatherium Floor of Thiongo Korongo site at Olduvai Gorge, where a regular pattern has been identified in the case of faunal remains and a cluster and regular combined pattern have been identified in the case of the lithics tools with an inhomogeneous L function (Panera et al. 2019) (this example is probably the most similar case to those results provided by Cuesta de la Bajada SPPs).

At the same time, it must be noted that Cuesta de la Bajada regular SPPs are also different from those SPPs provided by other typologies of sites, like referential camps created by modern (Domínguez-Rodrigo and Cobo-Sánchez 2017b) or by Lower (Domínguez-Rodrigo et al. 2017; Domínguez-Rodrigo and Cobo-Sánchez 2017a; Cobo-Sánchez 2020; Diez-Martín et al. 2021) and Upper (Vaquero et al. 2001; Vallverdú et al. 2010; Rosell et al. 2012; Fernández-Laso et al. 2020) Pleistocene hunter-gatherers or the pattern identified at a hunting camp dated during the European Middle Palaeolithic (Moclán et al. 2023), where clustered patterns are clearly identified. However, a similar regular pattern has been identified at another Middle Palaeolithic camp which has been characterized as a residential camp in the case of the access to cervids and as a hunting camp in the case of equids (i.e. Level P of Abric Romaní) (Marín et al. 2019).

Thus, it is clear that the SPP of Cuesta de la Bajada is especially interesting to understand the variability of spatial results provided by kill/butchering sites. In any case, the example of the Cuesta de la Bajada shows that kill/butchering sites do not necessarily have to show a cluster distribution, against other archaeological or ethnographical data (O’Connell et al. 1992; Hutson et al. 2020; García-Moreno et al. 2021), but may show some variability.

In this sense, it would be important to check in the future the reason for this variability, which may be influenced by different factors such as the type of taxon captured, the number of preys, or the intensity of the occupation, among others such as the seasonality of the occupation. At the same time, it could be interesting to check what is the reason for the existence of regular patterns in other types of sites and what is the possible relationship among them.

At this moment, we can only hypothesize the possible reason for the existence of this spatial distribution at Cuesta de la Bajada. Regular patterns are characterized by ‘segregation’ of points, and in most cases this aspect requires a causal mechanism that produces the segregation process (Baddeley et al. 2016). In this sense, from an archaeological point of view, it can be difficult to understand that anthropic activities that are repeatedly taking place at the same space produced segregation between materials instead of clustering. At the other end, postdepositional abiotic processes impacting any given assemblage produce CSR, not regular patterns (Domínguez-Rodrigo et al. 2017). However, if the formation process of the site is considered under a different perspective, the regular pattern documented at Cuesta de la Bajada can be understandable.

Cuesta de la Bajada was used as a kill/butchering site to which animal carcasses were transported during different events separated in time. These carcasses were subsequently butchered there. If this situation is considered, we think that it is logical that the regular pattern exists at long distances, because during the different occupation episodes separate occupation spots were probably produced in the same area. Thus, we hypothesize that the regular SPP of Cuesta de la Bajada is reflecting the palimpsestic nature (longitudinal sequence of depositional events) of the site. This situation is probably analogous to that identified at the level P of Abric Romaní or at the sites of TK and DS, where different occupation events produced debris in different occupation loci (Panera et al. 2019; Marín et al. 2019; Cobo-Sánchez 2020).

In addition, it must be noted that the segregation among materials has been identified only when a specific type of item is analysed (e.g. faunal remains), while if correlation between faunal and lithic remains is tested this is positive (i.e. both materials are clustered). In this sense, we think that the explanation for this situation is also related with the existence of a palimpsest with different anthropic occupations that resulted in non-concentrated remains in each of them. If a natural process had formed the site, the different types of materials would not have been ‘clustered’. This can only be the result of human activity during different occupation events. Thus, we think that the general pattern of the site is showing a SPP where different butchering events were produced. During these episodes, faunal and lithic remains were positive correlated since the lithic tools were used to butcher the carcasses, although the materials are segregated from those that were deposited during other butchering episodes. In the future, it will be of interest to carry out analyses of lithic and bone refittings that allow us to delve into these issues.

In addition, in the future, it is necessary to test this hypothesis with other kill/butchering sites, especially due to the differences between Cuesta de la Bajada and the horse assemblage of Schöningen 13-II 4 (Hutson et al. 2020; García-Moreno et al. 2021).

Specific aspects of CB2 and CB3

Conversely, if the spatial results obtained for CB2 and CB3 are carefully analysed, they can show that some differences are clear between both levels. This is especially interesting due to it could be useful in order to understand the variability of two kill/butchering sites that occupied the same space.

In general, the results provided by the analyses of the faunal remains of the level CB2 are showing that except for the medium-sized animals (mainly red deer), there are no clear high-intensity areas that could be interpreted archaeologically as activity areas or ‘clusters’ of anthropic activity.

In fact, if all faunal remains are considered at the same time (or if large-sized animals are analysed), the hot spot areas are similar to those provided by the presence of diagenetic alterations such as biochemical marks and manganese oxidations. This latter result can be interpreted as evidence of low alteration of the assemblage.

The analysis of anthropic (i.e. cut and percussion) and carnivore marks of this level have shown that the alterations appear to be more ‘clustered’ than other analysed variables. As in the case of Teixoneres (Zilio et al. 2021) or Amalda (Sánchez-Romero et al. 2020) caves, this situation could be used in order to define the sequence of occupation of both agents. However, the sample size of taphonomic marks at the level CB2 is too small and we prefer to be cautious about the obtained results. Nevertheless, it must be noted that the hot spot of carnivore marks appears to be located at the same area (i.e. north-western part of the site) and this could be indicative of carnivore ravaging of the carcasses left by hominins due to the assemblage have been produced in an anthropogenic way (Domínguez-Rodrigo et al. 2015).

On the other hand, if water-affected bones are considered (i.e. rounded and polished remains), it seems to be clear that the hot spots are more ‘clustered’ than in the previous analyses and this feature could be interpreted as a possible sign of water alteration of the assemblage. However, the application of the functions has shown that the SPP of rounded bones is regular and the SPP of the polished bones is random. Thus, we can reject the possible existence of a high degree of alteration of the faunal remains by water. These results are consistent with those related to the anatomical part representation provided by the previous taphonomic analysis (Domínguez-Rodrigo et al. 2015), which points out that the main anatomical parts represented at the site are those identified as the first that are extracted from the sites by water action (Voorhies 1969).

When lithic remains are considered analysing the different raw materials, high-intensity areas are clearer than in the case of faunal remains. In this case, the sample size can be used as an argument to note that maybe the hot spots cannot be considered as activity areas, but that is not the case when silicified limestone is analysed. This raw material is overrepresented in the assemblage (probably because it is the most suitable raw material for knapping in this environment) (Santonja et al. 2014) and the sample size can be considered sufficient for the interpretation process.

This raw material is showing a clear different distribution with three high-intensity areas where flakes, retouched tools, and cores can be found, and this characteristic cannot be a result of the alteration of water (Domínguez-Rodrigo et al. 2018). If silicified limestone tools with cortical surfaces are considered, there is a clear positive correlation between cortical cores, cortical flakes, and cortical retouched tools which could be interpreted as evidence of the preservation of the chaîne opératoire at knapping areas at the level CB2. In any case, the different high-intensity areas of flakes, retouched tools, and cores must be considered carefully. It is necessary to improve the analysis of these lithic remains by including a complete refitting analysis of the sample to confirm the possible presence of clear knapping areas.

As previously said, the other raw materials have smaller sample sizes, but when codependency among raw materials is tested most of the raw materials are clustered and this can be interpreted as a result of the anthropic activity and not as a result of natural processes.

Finally, the level CB2 has shown interesting results when lithic and faunal remains are considered at the same time. The existence of a cluster relation between them (although the general distribution is regular) is showing that the presence of both types of materials at the site has a common origin. Thus, hominins must be considered as the main accumulators of archaeological remains at the site, which have been preserved in a great manner. There is also a slight cluster relationship between water-affected materials if the Lcross function is used. This aspect can be interpreted as evidence of the existence of little areas affected by water that produced alterations in both faunal and lithic remains at the same time. However, we think that the alteration of the materials must have been slight because the cluster relationship is also identified when all faunal and lithic remains are compared.

Results provided by the CB3 level are similar to those provided by CB2 with some differences that must be commented on. In general, the distribution of faunal remains shows results extremely similar to the previous ones but it must be noted that here the hot spot maps of water-affected bones are showing clearer high-intensity areas which could be interpreted as evidence of further alteration of the level. However, again, the functions have shown that the SPPs are regular and not cluster and in this manner, we can reject the possibility of a high degree of alteration in those areas.

In this case, lithic tools are not showing a cluster relation among raw materials but here the presence of a cluster relation is present if raw lithic remains are compared taking into account the length of the tools. In this sense, we can consider a possibly higher degree of alteration of the sample which could modify the original position of some remains. However, rounded tools show again a regular distribution similar to those provided by the complete sample of lithic tools and we can reject the possible existence of a high degree of alteration of the site produced by water flows.

In this case, the identification of knapping areas is more complicated. The most common raw material (i.e. silicified limestone) seems to be distributed in the same way when different types of tools are considered and this allows us to be cautious about the results. However, in this case, we can note the case of ‘normal’ limestones, because flakes and cores seem to be differently distributed (nevertheless, the smaller sample size is important, and the results must be interpreted with care).

Again, the cluster relation between faunal and lithic remains has been identified (in this case it has been possible when the materials have been compared using the length of the remains). This situation shows again a low degree of alteration of the level. In the case of CB3, the existence of a positive correlation between those bones with percussion marks and those lithic tools used during percussion activities (anvils and hammerstones) (see Fig. 12) is especially important. This cluster relationship is showing that, as previously noted by the taphonomic analysis (Domínguez-Rodrigo et al. 2015; Moclán et al. 2022), there is a clear interest in the access to the bone marrow of the long bones to be consumed by hunters at the butchering site.

Finally, in the case of CB3, there are slight regular relationships between water-affected remains if the Kcross function is considered and cluster and regular if Lcross is used. These results are contrary to the result provided by the level CB2. Here we think that the degree of alteration of the sample is probably higher than in the case of CB2. As Santonja et al. (2014) showed, CB3 displays small floating gravels along the level, especially at the bottom. This evidence shows that in CB3 there would be more tractive than in CB2 and therefore, this could explain the higher degree of alteration. In any case, this degree of alteration is low in comparison with sites with high-energy conditions, such as the ancient levels of Cuesta de la Bajada (Santonja et al. 2000a, b; Santonja and Pérez-González 2001). It must be noted that here the anatomical representation of the site is mainly composed of elements which are easier to be altered by water flows (Voorhies 1969; Domínguez-Rodrigo et al. 2015).

Thus, we can be sure that the differences between both levels are particularly marked as far as the alteration of the remains by water flows. Both levels show SPPs that can be mainly attributed to the anthropic activity but in the case of CB3, the original SPP seems to be slightly more altered. However, this is not surprising. As previously noted (Santonja et al. 2014; Domínguez-Rodrigo et al. 2015), CB3 has a higher sample size and probably the archaeological assemblage was formed during a large period therefore bioestratinomic alterations such as resedimentation (Fernández-López 2000) could be more common at the assemblage than in the case of CB2.

Similarities between Cuesta de la Bajada and Schöningen

Above we insisted on stating that the SPPs of Cuesta de la Bajada and Schöningen 13II-4 are different, but if the specific results of Cuesta de la Bajada are considered it can be assumed that the SPPs of both sites show similarities. If we leave aside the horses’ results of Schöningen 13II-4, it is remarkable that the SPPs of cervids and bovines are quite similar to the SPPs of Cuesta de la Bajada (Hutson et al. 2020; García-Moreno et al. 2021). In this sense, the main question is if both assemblages have been produced similarly and if the equids of Schöningen 13II-4 show a different pattern, why do cervids and bovines of Schöningen show a similar pattern to those shown by Cuesta de la Bajada? Is there any subsistence answer to this question?

The geological features of both sites show that the sedimentation processes are different between them, although both of them are related with a lacustrine environment. In the case of Schöningen, the archaeological assemblages are related to an extremely low-energy sedimentary area characterized by the rise and fall of the lake level (Stahlschmidt et al. 2015) while in the case of CB2 and CB3 the sedimentary processes are related with the presence of a floodplain area where lacustrine areas appeared in relation with a small depression (Santonja et al. 2000a, b, 2014) and therefore the sedimentary processes should have higher energy due to the probably existence of small rills that could modify the archaeological assemblage in a higher degree than in the case of Schöningen.

Above we noted that although a slight water alteration has been identified at CB3 and especially at CB2, it is extremely weak and most of the assemblage still retains the original distribution pattern. This statement is also supported by the technological analysis (Santonja et al. 2014), which has pointed out the existence of all phases of the chaîne opératoire and the presence of debris elements (debris pieces were mainly collected during dry and wet sieving processes) (see Sesé et al. 2016 for a complete description of wet sieving processes). As previously said, the skeletal part profile of the levels is also indicative of a slightly modified assemblage (Voorhies 1969; Domínguez-Rodrigo et al. 2015; Moclán et al. 2022). At the same time, it must be noted that although Schöningen has almost ‘unbeatable’ archaeological conditions, this site has also suffered some alterations (see different areas defined by García-Moreno et al. 2021). Thus, if both sites have been only slightly modified by water, the spatial distribution of materials can only be related to different anthropic decisions.

Zooarchaeological analyses of the faunal assemblage of Schöningen 13II-4 have discussed if the accumulation of equids is the result of a single mass killing event or if the horses have been accumulated during multiple hunting events (Thieme 2005; Voormolen 2008; Van Kolfschoten et al. 2015a, b; Hutson et al. 2020, 2021; Bonhof and van Kolfschoten 2021). Recent approximations to this topic which included spatial and seasonal data (García-Moreno et al. 2021) stated that the accumulation of horses differed at Schöningen depending on the season, but in any case, the accumulation was related to multiple hunting events. The most interesting result about Schöningen is that the cluster pattern of equids was formed during summer-autumn periods and the occupation pattern was probably related to short and intense occupations which produced the cluster of equid carcasses (García-Moreno et al. 2021). However, the accumulation of bovines, cervids, and equids during winter-spring was more sporadic and produced a non-cluster/regular pattern at least in the case of cervids and bovines.

This situation opens new questions about the Cuesta de la Bajada assemblage due to in this site medium- (mainly cervids) and large-sized (mainly equids) animals (Domínguez-Rodrigo et al. 2015) show a regular distribution and high-intensity areas are difficult to be identified. If we use the model of Schöningen in order to interpret Cuesta de la Bajada, sporadic occupations can be proposed at the site (against short and intense ones) which were related to the access to different carcasses of equids and cervids during different episodes of occupation. In any case, it is necessary to carry out an analysis of the seasonality of Cuesta de la Bajada in the future in order to improve our knowledge about the origin of these assemblages which could be combined with an in-depth analysis of the processing and transport strategies of the carcasses of these animals.

Conclusions

The spatial analysis of the levels CB2 and CB3 of Cuesta de la Bajada has been very useful to clarify the possible variability of the spatial distribution of archaeological materials in a kill/butchering site. The SPPs of these levels have shown that clusters of archaeological materials are not present at the site although high-intensity areas can be identified especially in the case of lithic tools. Thus, the pattern has been identified as ‘regular’.

The spatial analysis has also shown that the alteration of the CB2 and CB3 assemblages by water flows or other bioestratinomic and diagenetic processes is weak. Thus, the spatial pattern identified during the excavation process can be considered as very slightly disturbed.

The analysis of the faunal remains has shown that, in general, the remains are mainly distributed in the north-western part of the site at both levels, regardless of taxa or the analysed taphonomic alteration. However, the analysis of the lithic industry has shown areas of higher intensity of materials depending on the type of tool, which could be indicative of the presence of knapping areas (especially in the case of the level CB2), although this aspect should be assessed in the future through a refitting analysis of the lithic assemblage. In addition, there is a clear positive correlation between fauna and the lithic industry, which demonstrates a clear relationship between the carcass acquisition by hominins during the occupation of the levels CB2 and CB3 and the knapping activities carried out at the site, which could be produced by hominins during the same episodes of occupation.

The general spatial pattern of this site is completely different from the main accumulation of the other European Middle Pleistocene kill/butchering site which has been analysed with spatial techniques (i.e. Schöningen 13II-4) (Hutson et al. 2021; García-Moreno et al. 2021). However, Cuesta de la Bajada and Schöningen are similar if the cervid and bovine assemblage of Schöningen is considered. This aspect could be indicative that Cuesta de la Bajada assemblages were formed during sporadic and short episodes of anthropic activity focused on the hunting of isolated animal carcasses which were subsequently transported to other sites (e.g. residential sites).

Nonetheless, it is still premature to evaluate the spatial significance of the disparities between Cuesta de la Bajada and Schöningen. It is crucial to explore these aspects more comprehensively by examining new sites to determine if the pattern of one of them is more typical of kill/butchering locations, while the spatial pattern of the other site may be the outcome of a specific factor related to the site’s conditions (such as topography or hydrography).