Using GIS and Geostatistical Techniques to Identify Neanderthal Campsites at archaeolevel Ob at Abric Romaní

Although intra-site spatial approaches are considered a key factor when interpreting archaeological assemblages, these are often based on descriptive, qualitative, and subjective observations. Currently, within the framework of research into spatial taphonomy and palimpsest dissection, several studies have begun to employ more quantitative and objective techniques, implementing tools such as geostatistics and geographic information system (GIS) methods. This is precisely the approach that the Abric Romaní team is following. In this work, we present GIS and geostatistics methods applied to the faunal and lithic assemblages from archaeolevel Ob, including an analysis of the spatial structure, the identification of clusters and sectors, size and fabric analyses, the projection of vertical profiles, and the reconstruction of a digital elevation model of the paleosurface. The results obtained indicate a clustered distribution, primarily concentrated into four dense accumulations. The predominance of remains < 3 cm in length and the absence of preferential orientations make it possible to rule out a generalised postdepositional movement affecting most of the site, although some local movement has been identified. The horizontal and vertical spatial analyses allow us to identify accumulations of a single material (lithic or faunal) in addition to mixed accumulations (lithic and faunal). Integrating all this data with the results of previous studies (zooarchaeological, refits, combustion structures, and partial lithic technological analyses), we evaluate and combine the interpretations proposed previously using different approaches, thereby improving the overall interpretation of the archaeolevel Ob. Finally, we also develop a preliminary comparison between Ob and some other levels at the same site (in particular M and P).


Introduction
Intra-site spatial approaches have long been recognised as a very useful, even necessary, tool when trying to understand and interpret Palaeolithic archaeological assemblages. Such approaches can provide key information on both anthropogenic activities and the actions of various natural agents, since all of these usually impact the spatial distribution of the remains. Consequently, intra-site spatial analyses provide valuable data on site formation (including both depositional and postdepositional processes) and past human behaviour.
Ethnoarchaeological studies have demonstrated the ability of humans to generate a recognisable spatial pattern and shown that the distribution of remains is related to the social behaviour of human groups (Yellen 1977a, b;Binford 1978a, b;O'Connell 1987;O'Connell et al. 1991; Bartram et al. 1991;Fisher and Strickland 1991;Kroll and Price 1991). This occurs because human activities are generally organised, that is, they occur in specific temporal and spatial sequences (Schiffer 1972(Schiffer , 1983Estévez et al. 1998;García-Piquer and Estévez-Escalera 2018;Domínguez-Rodrigo and Cobo-Sánchez 2017). However, certain intrinsic characteristics of the archaeological record hamper the use of ethnoarchaeological models for interpreting archaeological spatial patterns. For instance, most archaeological sites are palimpsests, understood as an accumulation and superposition of remains resulting from multiple different activities (both anthropogenic and natural), carried out over a variable period of time (Lucas 2005(Lucas , 2012Bailey 2007). As a consequence, the archaeological record has a much lower temporal resolution than the ethnoarchaeological work. To try to resolve this issue, some authors have opted for a spatiotemporal approach, introducing the concept of time into spatial analyses. The main goal for this is to dissect the palimpsest into archaeological associations with a higher temporal resolution, closer to "ethnographic time." To achieve this complex objective, a wide range of different disciplines, techniques, and methods have been developed and applied, including archaeostratigraphy, micromorphology, raw material unit (RMU) analysis, lithic and faunal refits, and computer-based tools (Chenorkian 1988;Canals 1993;de Loecker 1994de Loecker , 2004Audouze and Enloe 1997;Baxter et al. 1997;Baales 2001;Leesch et al. 2004;Vaquero et al. 2007Vaquero et al. , 2012Vaquero et al. , 2017Hovers et al. 2011;Enloe 2012;Carbonell 2012;Street et al. 2012;Turq et al. 2013;Machado et al. 2013Machado et al. , 2019Bisson et al. 2014;Mallol et al. 2013; Bargalló et al. 2016;Romagnoli and Vaquero 2016;Spagnolo et al. 2016Spagnolo et al. , 2019Spagnolo et al. , 2020aRomagnoli et al. 2018;Geiling et al. 2017;Gabucio et al. 2018a, b).
In terms of natural taphonomic modifications, various neotaphonomic studies and subsequent archaeological applications have shown the impact that certain natural processes, such as carnivore action or water flows, have on the location and disposition of previously deposited objects (e.g., Voorhies 1969;Behrensmeyer 1975;Frostick and Reid 1983;Frison and Todd 1986;Marean and Bertino 1994;Petraglia and Potts 1994;Bertran et al. 2006;Lenoble et al. 2008;Camarós et al. 2013;Arilla et al. 2014;Domínguez-Rodrigo et al. 2014Giusti and Arzarello 2016;Sánchez-Romero et al. 2016;Organista et al. 2017;de la Torre and Wehr 2018;Mendez-Quintas et al. 2019). Among these works, those that reproduce fluvial environments stand out, although the potential of other natural agents to mobilise pre-deposited items has also been studied. Recently, a new area of research in taphonomic studies, spatial taphonomy, has been developed, which underlines the importance of spatial proxies. Spatial taphonomy addresses the spatial analysis of taphonomic attributes by applying several multiscalar geostatistical methods (Domínguez-Rodrigo et al. 2018;Giusti et al. 2018). Its objective is to go beyond mere graphic representations, expanding the systemic vision of taphonomy and demonstrating that the spatial distribution of the attributes can help to identify the agents, as well as the superposition and interaction of different taphonomic processes (Domínguez-Rodrigo et al. 2018;Giusti et al. 2018).
Although intra-site spatial approaches are not new in archaeology, in most cases, these have been based mainly on visual, descriptive, and qualitative analyses, where subjective evaluations prevail (Bevan et al. 2013;Giusti and Arzarello 2016;Giusti et al. 2018;Domínguez-Rodrigo et al. 2018). This fact might limit the validity of the conclusions and, in addition, makes it extremely difficult to compare different assemblages. However, some authors have pushed the boundaries of these techniques, implementing statistical and, more recently, geographic information system (GIS) methods in their intra-site spatial studies (e.g., Whallon 1973Whallon , 1974Doran and Hodson 1975;Hivernel and Hodder 1984;Baxter et al. 1997;Bevan and Conolly 2006;Lenoble et al. 2008;Alperson-Afil 2008Crema et al. 2010;Keeler 2010;Moseler 2011;Gallotti et al. 2011;Nigst and Antl-Weiser 2012;Bosch et al. 2012;Bevan et al. 2013;Domínguez-Rodrigo et al. 2014Baxter 2015;Giusti and Arzarello 2016;Sánchez-Romero et al. 2016Spagnolo et al. 2016Spagnolo et al. , 2019Spagnolo et al. , 2020aDomínguez-Rodrigo and Cobo-Sánchez 2017;Organista et al. 2017;de la Torre and Wehr 2018;Giusti et al. 2018;Gonçalves et al. 2018;Discamps et al. 2019;Mendez-Quintas et al. 2019;Marín et al. 2020;Saladié et al. 2021). New areas of research, such as palimpsest dissection and spatial taphonomy, are employing these techniques in search of a more quantitative and objective approach. As a result, statistical tools previously implemented in other disciplines, including geology, ecology, and epidemiology, are being applied to intra-site spatial analyses. In addition, GIS methods, commonly applied to inter-site and regional studies, are increasingly being used in intra-site studies.
An approach of this type is being applied to the study of the Abric Romaní, one of the main key sites for the study of Neanderthal behaviour. The rock shelter, with a stratigraphy of approximately 30 m of mostly travertine sediments, has been dated to between 110 and 40 ka (Sharp et al. 2016) (Fig. 1). The site is characterised by a rapid sedimentation rate, estimated at about 0.46 mm/year (Vallverdú et al. 2012b;Vaquero et al. 2019), which is one of the features that make this site exceptional, conserving combustion structures, wood imprints, and minimising palimpsest formation. Almost the entire sequence corresponds to the Middle Palaeolithic, except for level A, which is Proto-Aurignacian.
The research team led by Prof. Eudald Carbonell has been working on this site for 40 years, extensively excavating and gathering information in great detail to document the maximum occupied area in each layer, over a surface area of approximately 300 m 2 (Carbonell 2012). The excavation system was designed with the final objective of reconstructing the behaviour of Neanderthals from a palaeoethnographic perspective. This goal conditioned the research methodology, both in the field and the subsequent studies, promoting a multi-and transdisciplinary approach.
These working methods, together with the rapid sedimentation rate, have enabled various intra-site spatiotemporal studies based on different archaeological levels from the site Vallverdú et al. 2005Vallverdú et al. , 2010Vaquero et al. 2007Vaquero et al. , 2012Vaquero et al. , 2015Vaquero et al. , 2017Vaquero et al. , 2019Fernández-Laso 2010;Rosell et al. 2012;Carbonell 2012;Bargalló 2014;Gabucio 2014; Bargalló et al. 2016Bargalló et al. , 2020bRomagnoli and Vaquero 2016;Modolo and Rosell 2017;Gabucio et al. 2018a, b;Marín et al. 2019;Fernández-Laso et al. 2020). Despite the rapid sedimentation, the palimpsests still have a much lower temporal resolution than that typical in ethnoarchaeological works, a fact that limits their direct comparison. Between the various levels, there are differences between the amount of remains and their distribution. Consequently, some show well-defined concentrations, produced in a reduced number of occupations, while others -such as level O -correspond to a more developed palimpsest. In order to dissect these palimpsests and bring them as close as possible to ethnographic time, a wide range of techniques, methods, and disciplines have been applied, including archaeostratigraphy, micromorphology, RMU analysis, lithic and faunal refits, computer-based tools, and mathematical methods.
Currently, the Abric Romaní research team is promoting a more quantitative approach to the intra-site spatial organisation of the archaeological record. The aim of this strategy, already initiated by Vaquero (1999), is to obtain more quantitative spatial data and thus favour a diachronic comparison between the different levels at the site. The approach has already been applied to some of the assemblages, including the lithic assemblage from level M (Romagnoli and Vaquero 2016) and the faunal assemblage from level P (Marín et al. 2019). In this work, we propose a similar proxy for archaeolevel Ob, including both lithic and faunal assemblages.
Abric Romaní level O, dated at around 54.6 ka BP (Vaquero et al. 2013), was excavated between 2004 and 2011, and provided more than 40,000 archaeological remains. The study of this level has led to numerous publications (Vallverdú et al. 2012a;Courty et al. 2012;Gabucio et al. 2012Gabucio et al. , 2014aGabucio et al. , b, 2018aLópez-García et al. 2014;Picin et al. 2014;Chacón et al. 2015;Bargalló et al. 2016Bargalló et al. , 2020aCarrancho et al. 2016;Allué et al. 2017;Fernández-García et al. 2018Gómez de Soler et al. 2020a, b;Eixea et al. 2020Eixea et al. , 2021. Some of these (Bargalló 2014;Gabucio 2014;Bargalló et al. 2016;Gabucio et al. 2018a), based on faunal and/or lithic remains and aimed at dissecting the original palimpsest from a spatiotemporal perspective, led to the identification of two main archaeolevels, Oa and Ob (Fig. 2). In addition, for both archaeolevels, the authors proposed a site function (mainly as a campsite), intra-site structure (different activity areas), connections between areas (long connections between the E and the W of the rock shelter in Ob), and settlement dynamics (short occupations in Oa; a higher number of individuals and possibly longer occupations in Ob) (Bargalló 2014;Gabucio 2014;Bargalló et al. 2016Bargalló et al. , 2020aGabucio et al. 2018a, b). The lithic and faunal results were obtained separately, but the authors shared information with each other throughout all the phases of the study. Nevertheless, comparisons between lithic and faunal results can be difficult since the methodology applied, although very similar, is not exactly the same. In addition, the analysis of the lithic technology from archaeolevel Ob is still in progress.
In this work, we present the first application of GIS and geostatistical methods to archaeolevel Ob. This includes a determination of the spatial structure (basically Besag's L (r) function), the identification of clusters and sectors (Kernel density and k-means analysis), the exploration of size and orientation distribution patterns, the visualisation of the vertical distribution of the remains, and a palaeotopographic reconstruction using a digital elevation model. In this way, we contribute to the quantitative approach that is currently being developed at the site and facilitate its diachronic study. Unlike previous works within this area, we selected a specific archaeolevel, Ob, identified as a result of a previous exhaustive archaeostratigraphic study. Another singularity of this work is the simultaneous application of a great variety of analyses to both lithic and faunal remains.
In summary, our main objectives are to: (1) characterise the point pattern of the lithic and faunal remains in Ob; (2) contribute to the dissection of the palimpsest; (3) improve our knowledge of the deposit's formation processes and postdepositional disturbance; (4) provide an optimal basis for the integration of lithic and faunal data; (5) produce quantitative data that facilitates the diachronic study of the site; and (6) take the first steps towards such a diachronic comparison.

Materials and methods
This study focuses on Ob, the most extensive and complex archaeolevel in level O. It varies in thickness, being thinner in the inner area (maximum 25 cm) and reaching 55 cm in the outer area of the rock shelter. A similar trend is observed in the slope, which is 12.92° towards the theoretical SW, although it is almost flat in the inner area. A total of 21,748 piece-plotted lithic remains and 8596 piece-plotted faunal remains have been recovered from Ob. Until 2010, all the lithic and faunal remains detected during the excavation were piece-plotted. Instead, as of 2011, it was decided to piece-plot lithic remains > 1 cm and faunal remains > 2 cm. Figure 3 shows the location of the travertine blocks, the wood imprints, and the 50 combustion structures assigned to the archaeolevel Ob.
Spatial analyses and geostatistical techniques were applied considering both lithic and faunal remains. Cartesian coordinates are available for all the lithic and faunal remains that were recorded during the excavation. Any remains for which the spatial position was not recorded (recovered from the sediment after water-sieving using nested meshes of 5 mm and 0.5 mm) have not been taken into account in this study. First, we analysed the two assemblages (faunal and lithic) together, to better characterise archaeolevel Ob, and then separately, to compare the two materials. This should enable a clearer understanding the formation processes and facilitate the future integration of the faunal and lithic results, once the latter is fully available.
The spatial and geostatistical analyses were performed using R and its environment (R Core Development Team 2011) and (QGIS 2016) software (version 2.18 Las Palmas and version 3.26 Buenos Aires). Occasionally, we used other open-source software, such as CrimeStat software (National Institute of Justice, Washington, DC, USA) and SAGA GIS. The following subsections detail the methods applied to the specific analyses.

Spatial structure
Besag's L (r) function was used to to characterise the spatial structure of the three assemblages (lithic and faunal remains/ lithic remains/faunal remains) by specifying whether the points had a uniform, random, or aggregated distribution (Besag 1977). Complete spatial randomness was used as the null hypothesis (Diggle 1983). This procedure is key in point pattern characterisation, and is now being used more ◂ systematically in intra-site spatial studies, although in most cases Ripley's original K (r) function (Ripley 1976(Ripley , 1977 is selected (e.g., Bevan and Conolly 2006;Crema et al. 2010;Giusti and Arzarello 2016;Romagnoli and Vaquero 2016;Spagnolo et al. 2016Spagnolo et al. , 2019Spagnolo et al. , 2020aThacher et al. 2017;Domínguez-Rodrigo and Cobo-Sánchez 2017;de la Torre & Wehr 2018;Domínguez-Rodrigo et al. 2018;Gonçalves et al. 2018;Giusti et al. 2018).
We prefer the modified function L (r) proposed by Besag (1977) to the original Ripley's K (r) function because the former is easier to interpret, especially when the function is represented with respect to 0. In the case of aggregated distributions where this type of representation is used, the highest point of the observed curve clearly indicates the approximate radius of the clusters. Moreover, Besag's L (r) function is normalised, so it allows the point pattern of assemblages with different number of points to be compared.
Although we prefer Besag's L (r) function, with the intention of broadening our understanding of the data and facilitating future comparisons, we also calculated Ripley's K (r) function and added the results into the Supplementary Information . Likewise, we included Monte Carlo simulations of the tests (Robert and Casella 2004), generating 50 sets of artificial points with the same density  Gabucio et al. 2018a). Bottom: Horizontal distribution of lithic and faunal remains according to the rock shelter limit as the observed points but with random spatial locations . In this way, we checked if the K and L values deviated with respect to the expected theoretical K and L values for random point patterns.
The "isotropic" edge correction was preferred as this was developed for irregular polygons (Ripley 1977), such as the excavated area of level O. However, in the Supplementary Information, we provide the calculations of Ripley's K (r) function and Besag's L (r) function with both isotropic and border corrections .
Before using Besag's L (r) function, a homogeneity test was carried out on the three assemblages to be able to determine the most appropriate method (homogeneous or inhomogeneous). This step is not always performed, since homogeneous formats (which tend to provide more powerful and clearer results) are sometimes selected directly. Thus, although inhomogeneous tests are considered to be better adapted to the features of archaeological sites, in practice they are little used. In this regard, Domínguez-Rodrigo and collaborators intentionally used the two methods with specific objectives: the homogeneous test to show spatial trends; and the inhomogeneous one to correct the former

Clustering pattern
As a first step, we used density maps generated using the kernel algorithm as an informal approach to spatial cluster analysis (Baxter et al. 1997). In this task, we considered the radius suggested from the results of Besag's L (r) functions. Kernel density estimation (KDE) uses point data to create a smoothed density map taking into account the proximity between the points (Baxter et al. 1997). We used QGIS software to create the density maps, generating both 2D and 3D maps. KDE has been applied to different archaeological contexts in several previous studies, including recent examples relating to Palaeolithic sites (e.g., Keeler 2010;Moseler 2011;Nigst and Antl-Weiser 2012;Blasco et al. 2016;Spagnolo et al. 2016Spagnolo et al. , 2019Spagnolo et al. , 2020aRomagnoli and Vaquero 2016;Sánchez-Romero et al. 2016de la Torre and Wehr 2018;Giusti et al. 2018;Marín et al. 2019).
As a second step, a k-means analysis was performed, after calculating the ideal number of clusters using the elbow method (Kintigh and Ammerman 1982;Baxter 2015;Kasambara 2017). In elbow graphs, the point where the curves abruptly change trend indicates the ideal number of clusters, which must be read on the x-axis. We used the R language and environment to calculate the ideal number of clusters (R Core Development Team 2011), and CrimeStat software (National Institute of Justice, Washington, DC, USA) to carry out the k-means analysis. This analysis allowed us to classify all the items in different sectors and identify the main clusters (one cluster within each sector). Then, we visualised the results in QGIS. K-means analyses have been widely used in archaeology since the 1970s, although their use in Palaeolithic contexts came somewhat later (Doran and Hodson 1975;Ammerman et al. 1983;Kintigh and Ammerman 1982;Simek and Larick 1983;Simek 1984;Rigaud and Simek 1991;Baxter et al. 1997;Vaquero 1999). Nowadays, the use of k-means analysis is declining, although some authors still use this method, complementing it with other techniques, or use similar methods (Baxter 2015;Sánchez-Romero et al. 2016;Domínguez-Rodrigo and Cobo-Sánchez 2017;Mendez-Quintas et al. 2019;Sánchez-Romero et al. 2021).
In order to gain a better understanding of the differences and similarities between the distinct materials, all these methods were applied first to the combination of lithic and faunal remains and then for the two assemblages separately. However, we then proceeded to examine the remains in greater detail, taking into account only the sectors and clusters calculated for the combined lithic and faunal remains. For instance, considering the area and the number of remains (lithic/faunal/ both), we calculated the average intensity of each of these sectors and clusters, as has already been performed for level M at the same site (Romagnoli and Vaquero 2016).

Size distribution pattern
Both faunal and lithic remains were classified into the following size classes: ≤ 1 cm; > 1-2 cm; > 2-3 cm; > 3-5 cm; and > 5 cm. These measurements correspond to the longest axis of each item, regardless of their technical or anatomical features. Despite the fact that many remains ≤ 1 cm were surely not noticed nor piece-plotted during the excavation (especially the bones recovered from the year 2011, due to the change in the cut-off size for piece-plotting), we decided to add the category ≤ 1 cm because the amount of small remains is huge in the assemblage. However, we will take this bias into account when interpreting the results.
Next, we explored the horizontal distribution of these categories using two methods. The first consisted of generating density maps (KDE, representations in 2D and 3D) for each size class, taking into account the two materials studied (lithic/faunal). The second method was based on quantifying the number of remains in each size category and material contained in the different sectors and clusters identified using the k-means analysis. This bimodal approach allowed us to address the size distribution pattern in both a visual and quantitative way.
The distribution of the remains according to size categories is considered very informative in terms of assessing the degree of postdepositional disturbance in archaeological deposits (Petraglia and Potts 1994;Bertran et al. 2006Bertran et al. , 2012Spagnolo et al. 2016Spagnolo et al. , 2020ade la Torre et al. 2018;Mendez-Quintas et al. 2019). For instance, water flows tend to displace the lightest pieces, which are often the smallest (although other factors are also involved, such as density and shape). Consequently, deposits in fluvial environments usually have scarce, or even no, small remains. In contrast, when the remains are dispersed by gravity, it is the heavier remains (usually the largest) that move downslope.

Fabric analysis
During the fieldwork, the general orientation data (NS, NESW, EW, NWSE, or SQUARE -in the case of items with the same length and width) of each remain was recorded prior to its removal. To assign an orientation to the remains, the theoretical north of the site, located on the wall at the bottom of the rock shelter, was taken as a reference. In line with previous work (Bertran and Texier 1995;Bertran and Lenoble 2002;McPherron 2005;Benito-Calvo and de la Torre 2011;Domínguez-Rodrigo et al. 2014;Sánchez-Romero et al. 2016;de la Torre and Wehr 2018;Giusti et al. 2018;Spagnolo et al. 2020a, c), we selected a subsample of remains appropriate for the orientation analysis, comprising those larger than 2 cm and with an elongation index (length/ width) > 1.6. As in the case of the size categories, we generated 2D and 3D density maps for each orientation (NS, NESW, EW, and NWSE) and quantified the number of items with each orientation in the different sectors and clusters. In both cases, we differentiated lithic and faunal remains.
Fabric analysis is particularly relevant for detecting reworking processes, since different processes lead to different orientation patterns (Voorhies 1969;Petraglia and Potts 1994;Bertran et al. 1997;Bertran and Lenoble 2002;Lenoble and Bertran 2004;McPherron 2005;Benito-Calvo and de la Torre 2011;Sánchez-Romero et al. 2016;de la Torre and Wehr 2018;Giusti et al. 2018;Mendez-Quintas et al. 2019;Spagnolo et al. 2020a, c). For instance, generalised water transport usually generates preferential orientations. Thus, in fluvial contexts, the longitudinal axis of most small remains will be aligned with the direction of the flow, whereas the larger remains tend to align transversally to this. In contrast, debris-flow deposits usually present massive and poorly bedded mixtures of unsorted sediments and random clast orientation (except at the flow margins). Finally, undisturbed archaeological assemblages tend to have a planar orientation pattern.

Vertical distribution
To observe the vertical distribution of the remains, a total of 13 profiles were drawn using the Z-values. These included both longitudinal and cross-sectional profiles, covering the complete surface area of the rock shelter. All profiles were plotted with a thickness of 20 cm. In each one, lithic and faunal remains were differentiated (different colours) in order to compare their vertical distributions. In addition, the location of the different sectors and clusters was indicated in all the profiles (considering the lithic and faunal remains together).
The generation of profile projections is considered a key factor in the study of site formation processes, although it is not applied systematically in all studies (Cacho et al. 2016;Giusti and Arzarello 2016;Giusti et al. 2018;Discamps et al. 2019). For instance, vertical profiles show the thickness and slope of the layer, the possible existence of different sublayers with a greater temporal resolution, and the morphology of the palaeosurface on which the layer accumulated. In addition, vertical distribution analyses help to identify massive processes and differentiate between autochthonous and allochthonous assemblages.

Digital elevation model (DEM)
We used the depth values of the remains as a proxy to reconstruct the morphology of the base of the archaeolevel Ob. From a database containing the Z-values of both lithic and faunal remains, a grid file was generated using SGA GIS software. During this process, we indicated that, in the case of multiple Z-values for the same pixel, the lower value would be selected. Next, the grid file was interpolated (Multilevel B-Spline Interpolation) using the same software. The result of the interpolation was represented in 3D by applying the QGIS tool "QGis 2threejs." Finally, using the same interpolation, a terrain analysis (Analytical Hillshading) was performed. As a result, we obtained detailed 2D and 3D images of the surface of archaeolevel Ob.
Similar methods have been used by other authors to produce palaeotopographic reconstructions (e.g., Giusti et al. 2018;Bargalló et al. 2020a;Spagnolo et al. 2020a, b, c;Sánchez-Romero et al. 2020). Determining the relief of the surface on which the archaeological remains are preserved is important when interpreting the spatial distribution of the materials, since the existence of steep slopes, depressions, high points, and so on can influence the location of the remains and some of their traits, such as orientation and slope.
Besag's L (r) function showed that the point patterns were clustered (Fig. 4), ruling out the null hypothesis of complete spatial randomness. The point pattern is clustered both when the lithic and faunal remains are considered at the same time and when they are processed separately. The presence of the observed line (L value) above the poisson line in all the three cases indicates this. However, the results suggest that the faunal remains are more strongly clustered.
In the graph for the two materials together, the maximum height of the observed curve suggests a radius of about 80 cm for the clusters. A similar result (between 60 and 80 cm) is observed in the graph of faunal remains. In contrast, the graph of the lithic remains suggests a radius of approximately 120 cm.
The Supplementary Information includes the results of the homogeneous and inhomogeneous tests (Ripley's K (r) function and Besag's L (r) function), using isotropic and border corrections, and their corresponding Monte Carlo simulations . In all the homogeneous tests, Monte Carlo simulations confirm the clustered pattern of the assemblages, whereas in the inhomogeneous tests, they only confirm the clustered pattern of the faunal assemblage. In no case do the results indicate a random distribution of the remains.

Clustering pattern
We plotted all the lithic and faunal remains horizontally (X/Y) and then compared the distribution of these remains with the location of the structural elements (combustion structures, travertine blocks, and wood imprints). The resulting scatterplots show that the archaeological remains tend to cluster around and within the hearths, especially those located in the inner half of the rock shelter (Figs. 3 and 5). However, there is also a very strong tendency of the remains to concentrate towards the theoretical NW (grid squares R-U/58-62, approximately) and N (grid squares U-W/51-54) of the rock shelter (Fig. 5).
Although the general distribution pattern is similar, there are some significant differences between the lithic and faunal . The x-axis represents the distance (cm), the continuous black line represents the observed Besag's L (r) function, and the discontinuous red line represents the poisson line ▸ remains (Fig. 5). For instance, lithic remains are more abundant in the NW corner, where they clearly predominate over faunal remains, while bones are more abundant in the N, where the proportions between both materials are more equal (even faunal remains predominating over the lithic remains in grid squares V/52-53). In addition, there is a well-defined cluster towards the theoretical W (approximately grid squares Q-P/57-60) composed almost solely of lithic remains. Finally, the outermost area, especially towards the E, contains proportionally more faunal than lithic remains.
We can better appreciate these differences using kernel density maps, especially when displayed in 3D (Fig. 5). Density maps were developed choosing a radius of 80 cm, as suggested by Besag's L (r) functions. If we look closely at the theoretical NW corner, we can identify three highdensity points towards the NW. This is more evident in the case of lithic remains, although they can also be appreciated from the distribution of faunal remains.
Once we had confirmed that the spatial structure is clustered and before performing the k-means analysis, we explored the ideal number of clusters by means of the elbow method. We applied this method to the faunal and lithic assemblages both together and separately, in order to compare the results (Fig. 6). Considering only the lithic remains, the elbow method suggested four or five clusters. When only the faunal remains were taken into account, this method suggested that the ideal number of clusters was three or four. Finally, when both materials were analysed together, the elbow method indicated four clusters. Consequently, we decided to conduct the k-means analysis selecting four clusters, classifying the assemblage according to four sectors. Figure 6 shows the distribution of these sectors considering lithic and faunal remains together and both assemblages separately. Light-coloured polygons correspond to sectors whereas dark ellipses correspond to the main clusters. The results obtained from all the remains are practically identical to those obtained from lithic remains alone (which are quantitatively more numerous than faunal remains). Here, three clusters are concentrated in the NW area, whereas the fourth cluster extends through the eastern half of the rock shelter, coinciding with the four high-density points visible on the density maps (Fig. 5). In contrast, the results obtained from the faunal remains alone are quite different: there are two clusters in the W, one in the N and another in the E.
In the following analyses, the classification into four sectors and the clusters obtained considering the lithic and faunal assemblages together were taken as a reference. We named the sectors KM1-4, starting from the W and going clockwise (Fig. 6). Likewise, the main clusters were named by adding a "c" to the beginning of the sector name to which they belong (cKM1-4).

Fig. 6
Results of the elbow method and k-means analysis calculated first for lithic and faunal remains separately, and then for the two assemblages combined. The blue and red lines in elbow method correspond to the inter-and intra-class curves. The labels on the k-means analysis considering the two assemblages together correspond to the names of the sectors (light-coloured polygons) and accumulations (dark ellipses) The first step in the characterisation of these sectors and clusters involved calculating their average intensity (Table 1). Considering the sectors, both the lithic and faunal remains were concentrated more intensely in KM2 (which is the smallest sector), followed by KM3, KM1, and, at some distance, by KM4 (which has a larger area). If we look at the clusters, we can see that the greatest intensity of lithic remains was accumulated in cKM2, followed by cKM1, cKM3, and cKM4 (which is the cluster with the largest area), while the accumulations of faunal remains had similar intensities in cKM3 and cKM1, followed by cKM2 and cKM4 (where the extent of the area blurs the large concentration of bones observed in the N of the rock shelter, coinciding with the NW area of cKM4). Figure 7 presents the density maps (in 2D and 3D) of the complete assemblage by length categories. In contrast, Figs. 8 and 9 contain the density maps related to lithic and faunal remains, respectively. These figures can be visually compared to Fig. 5, which contains the density maps of all the remains (lithic and faunal remains together and separately). In addition, Table 2 provides quantitative information on the distribution of both materials by length categories, sectors, and accumulations.

Size distribution patterns
The lithic remains ≤ 1 cm in length are concentrated in two points in the NW of the site, coinciding with clusters cKM1 and cKM3 (Fig. 8). If we compare this with the general distribution of all the lithic pieces ( Fig. 5), we can see that, although many lithic remains are concentrated in cKM2, the items smaller than 1 cm are scarce in that area. This fact is even more evident if we look at Table 2. The distribution by sectors and accumulations of the lithic remains that are between 1 and 2 cm long, which are the most abundant, is very similar to that of the total lithic remains: there are significant concentrations at three points in the NW (cKM1-3) and another to the N (several points inside cKM4). However, as the length category increases, the global number of remains decreases and their proportion in KM2 gradually increases. Thus, of the 734 lithic remains > 5 cm, 244 are in KM2 (33.24%). Likewise, Fig. 8 indicates that several remains > 3-5 cm and > 5 cm are located in the outermost fringe of the rock shelter, where the total number of remains is very low.
As far as faunal remains are concerned, Fig. 9 shows that the remains ≤ 1 cm in length are mostly concentrated to the N, mainly at two points (one much more pronounced than the other) within the cluster cKM4. There are also two lesser points in the NW zone, coinciding with the clusters cKM1 and cKM3. However, as in the case of the lithics, cKM2 presents a low density of small remains. Table 2 illustrates this: the proportion of remains ≤ 1 cm is very high in sector KM4 (26.71%, 33.24% in cKM4) and very low in sector KM2 (5.50%, 5.75% in cKM2). In fact, of the 1599 faunal remains ≤ 1 cm in length, 1037 (64.85%) are in KM4. Both the figure and table indicate similar dynamics with respect to faunal remains > 1-2 cm: there is a greater density of remains in the N zone, followed by the NW, but with an increase in the proportion of remains in the NW area (especially in clusters cKM1 and cKM3, but also cKM2). In the higher length categories, which have increasingly fewer items, most remains are concentrated in the NW area, mainly in the KM2 sector, followed by the KM3 sector (but not in the cKM3 accumulation).

Fabric analysis
Figures 10, 11, and 12 contain, respectively, the density maps by orientation categories for the complete assemblage, the lithic sub-assemblage, and the faunal sub-assemblage. Likewise, Table 3 shows the number and percentage of remains by orientation, material (lithic and faunal), sector, and cluster.
Both the figures and the table show a fairly similar proportion of the four orientations (NS, NESW, EW, and NWSE)  Nevertheless, there is a mild general tendency for the NS and EW orientations to be slightly better represented than the NESW and NWSE ones (Table 3). This trend is more evident among lithic remains than among faunal remains.
In Table 3, we can see some small differences between the lithic and faunal subsets. For example, for lithics, in all sectors and clusters, the lowest percentage corresponds to the NWSE orientation, while in the case of faunal remains, the lowest percentage sometimes corresponds to the NWSE orientation (sectors KM1-3 and cluster cKM3) while other times this is represented by the NS orientation (sector KM4 and clusters cKM1, cKM2, and cKM4).
There are also some discrepancies in the orientations with higher percentages. On the one hand, the lithic remains present a higher proportion of NS orientations in KM1 and KM4 and the cKM2 and cKM4 clusters, while KM2 and KM3 and the cKM1 and cKM3 accumulations have a higher proportion of EW orientations. On the other hand, the faunal remains show a higher percentage of NS orientations only in KM3 and cKM3, with EW orientations being more abundant in all other sectors and clusters.
The strongest contrast between the orientations of the lithic and faunal assemblages is found in KM4, especially in cluster cKM4, where NS orientations are the most frequent among the lithic remains (33.21%) but the rarest among the faunal remains (18.21%). However, it should be considered that, in general, the differences between the percentages of each orientation are low, with cKM4, the most extensive cluster, presenting the greatest difference. Finally, it should also be mentioned that, of the ten slope categories determined during the field work (N, NE, E, SE, S, SW, W, NW, vertical, and flat), the flat slope clearly predominates, this being the case for approximately 30% of the remains in each of the clusters (KM1-4).

Vertical distribution
A total of 13 profile projections were made, five oriented EW (A-E), and eight NS (F-M) (Figs. 13, 14, and 15). The vertical view of these profiles shows that the inner half of the rock shelter (approximately lines Q-W) is very flat. Thus, in the NS-oriented profiles, a change of slope is observed between the inner zone, flat, and the outer zone, more inclined. This change is seen towards the Q line in the western area of the site (profiles G, H, and I), towards the wall of the rock shelter in the central zone (in section K, the change occurs on line U). In the eastern area, this change is less noticeable (not observed in section L; in section M, it seems to occur in line N).
In both the EW-oriented profiles and, especially, the NSoriented ones, the Ob archaeolevel only shows internal stratification at some points. Thus, in the NW (extreme N end of sector KM1 and sectors KM2 and KM3) and N (W end of sector KM4) zones, two layers can be observed. These layers were also reported in previous publications Gabucio et al. 2018a, b), where they were named microlevel Ob1 (upper) and microlevel Ob2 (lower). Unlike archaeolevels Oa and Ob, which are separated by a continuous sterile layer, microlevels Ob1 and Ob2 are separated by a discontinuous sterile layer. While the N zone of Ob2 contains more remains than Ob1, in the NW zone, the number of remains from the two microlevels seems more equal. These microlevels have not been detected in any other area of archaeolevel Ob. Finally, at some points, a few  (in 2D and 3D) showing the distribution of remains (lithic and faunal) according to orientation categories remains were found at a level significantly lower than the others, but these do not present any type of lateral continuity (profiles B, G, and H).
The vertical distribution also shows differences between the faunal and lithic remains. Thus, concentrations of a single material were detected in some profiles, some of which had already been identified when analysing the horizontal distribution of the remains. For example, in profiles G, H, and I, an accumulation of lithics was observed in squares Q-P/57-60. Likewise, profile I also showed a cluster of faunal remains in squares M/58-59. In addition, the large concentration of bones identified in UW/51-54 (especially in V/52-53) Fig. 11 Density maps (in 2D and 3D) showing the distribution of lithic remains according to orientation categories appeared clearly in profiles A and K. Most of the remains from this last faunal accumulation belong to microlevel Ob2. The profile projections also revealed another difference between the distribution of lithic and faunal remains that could not be detected by analysing the horizontal distribution. Thus, profiles J and K showed an inverse proportion of the materials by microlevels, while in Ob1, the lithic remains apparently dominated over the faunal ones, and in Ob2, the faunal remains seemed to predominate over lithics.
In order to integrate the combustion structures into the study of the vertical distribution of the remains, we have classified the structures from the NW and N zones into the two  (Table 4). In turn, this has allowed us to relate these combustion structures to the archaeological material (lithic and/or faunal) associated with them. In both microlevels, there are medium-sized hearths separated from the wall by a couple of metres and associated with a large quantity of both lithic and faunal remains. This is the case of the combustion structures AR08-10-11/1 (Ob2) and AR07-08-11/1 (Ob1). Smaller hearths closer to the rock shelter wall have also been documented in both microlevels, such as AR11-1 (Ob2), AR11-8, AR11-9, AR11-11, and AR11-12 (Ob1). These combustion structures, unlike the previous ones, are usually associated with a small number of remains. Also noteworthy is the identification in some structures of different combustion phases, each phase coinciding with a different microlevel. This occurs in AR11-10, where both phases are similarly associated with both lithic and faunal remains. In contrast, in AR06-07-10-11/1, each of the two phases is associated with different archaeological material: with a predominance of lithic remains in Ob1 and an accumulation of bones (the one identified in grid squares V/52-53) in Ob2.
Finally, the visualisation of the profiles also made it possible to trace the morphology of the lower limit and the thickness of the more relevant clusters (cKM1-4). As for clusters cKM1-3, profiles B, C, D, and E showed that these accumulations are fairly thick, but when viewed more closely, for instance through profiles F, G, H, and J, it could be seen that this thickness is due more to the existence of two microlevels (with a discontinuous sterile layer in the middle) than to the thickness of each of the microlevels. Profile F, in addition, showed how in that area both microlevels describe a shape similar to that of a cuvette, joining cKM1 and cKM2. In fact, many remains from Ob1 were in one of the deepest areas of this microlevel, which coincides partially with cKM1. Likewise, most of the deepest remains from KM2 were accumulated in cKM2. Profiles C and D suggest that the western part of the cKM3 accumulation could be part of this cuvette. Inside the cuvette, the remains from sector KM2 had Z-values that were generally somewhat higher than those of sectors KM1 and KM3 (projections B, C, D, F, and G). Profiles I and J indicated that the central part of KM3 presented no concave morphology and it was in quite a high and flat area. However, the zone between KM3 and KM4 seemed somewhat more depressed than the zones around it towards the E and W (projections C and D). Regarding cluster cKM4, all the EW-oriented profiles (A, B, C, and D) showed a very flat surface, although in C, an accumulation of remains towards T51 coincides with the lower Z-value in the part of the profile corresponding to cKM4. In addition, in profile K, it could be seen how the huge accumulation of bones was in a fairly flat area that did not coincide with the most depressed zone of cKM4.

Digital elevation model
Both the 2D and 3D representations of the palaeosurface confirmed that the slope in the inner part of the rock shelter was slight (Fig. 16). Likewise, they confirmed that the change towards a more pronounced slope occurred towards the Q line in KM1, between the T and U lines in KM3, and the W half of KM4, while the E half of KM4 maintained a flatter surface and only began to descend significantly from the lines L-M. This change in slope became more pronounced towards the theoretical W of the rock shelter, highlighting the depressed area towards the S of cKM1. Figure 16 shows the elliptical morphology of the NW zone, where the cluster cKM2, and part of the cKM1 and cKM3 clusters were located. This morphology fits well with the distribution of the remains presented in Fig. 5: the distribution of both lithic and faunal remains (although it is more evident in the latter) describes an oval, with a much emptier central space.

General spatial pattern of the faunal and lithic remains
One of the main contributions of this work is the simultaneous study of the point patterns of lithic and macrofaunal remains, enabling their comparison and allowing common conclusions to be drawn from the two materials. The application of different GIS and statistical methods (scatterplots, KDE, Besag's L function, and k-means analysis) to the lithic and faunal remains separately as well as to the two assemblages together allowed us to distinguish similarities and differences between the distributions of the two materials. Both assemblages show a clustered spatial structure, although the faunal remains show a stronger trend in this sense.
Both faunal and lithic remains tend to be concentrated towards the inner zone of the rock shelter. This trend is especially strong when considering only the remains ≤ 2 cm in size. The densest clusters are in the NW zone (cKM1-3) and N zone (cKM4). However, while lithics clearly dominate over faunal remains in cKM1-3, cKM4 presents a more equal proportion between both materials (although lithic remains are still more abundant than faunal remains). In addition, a predominance of faunal remains is observed in the external area, especially in the theoretical SE of the rock shelter.
In terms of the distribution of faunal remains in the inner area, however, a previous study focused on the macrofauna from level O (Gabucio et al. 2014b) qualifies these results. That study (where contour plots were used to locate different taphonomic and zooarchaeological variables, including categories other than length) showed that, if the remains recovered from bags of general findings and wet sieving are taken into account (i.e., the remains recovered but not pieceplotted), the densest accumulation of faunal remains is in the NW corner (as we see here in the lithics), followed by the N zone. This mismatch could be related to the fact that the smallest remains are the most sensitive to factors such as differences between different field methodologies, such as the change in the cut-off size for piece-plotting. In addition, the smaller remains are more visible (and therefore more often piece-plotted) when grouped into well-defined accumulations, as is the case of the accumulation of bones identified Fig. 13 Situation of the profile projections made for the Ob archaeolevel in relation to the distribution of the remains in the N of the rockshelter (cKM4). These remains, moreover, were found inside a combustion structure, which was excavated very carefully, and they presented striking taphonomic features (they were burnt, mostly calcined) favouring their visibility (Gabucio et al. 2014b(Gabucio et al. , 2018a.
Despite the differences observed between the lithic and faunal spatial patterns, as lithic fragments are more numerous than bones, the sectors and clusters calculated from the combined materials almost mirror those calculated using only lithics while they are quite different from those calculated using the faunal remains. This fact may have distorted the study of the faunal assemblage distribution, which was not carried out according to its ideal subdivisions. However, we believe that it was worth doing, since in this way, all the remains were classified into the same categories, allowing the two materials to be directly compared in terms of their dimensions, orientations, and so on. Likewise, the extent of sectors KM1 and KM4 and the high density of the clusters recognised using KDE and k-means made it difficult to identify and study more modest accumulations, which have had to be explored using other techniques, such as vertical analyses.
Both the horizontal spatial analysis and the vertical projections allowed us to identify accumulations with mixed materials and others formed exclusively (or almost exclusively) of either faunal or lithic remains. The main clusters in the NW zone, cKM1-3 (except the W zone of cKM-3), present a mixture of the two materials in both microlevels, although lithics are generally more abundant. In the NE zone of the KM4 sector, a similar mixed distribution of materials is observed, this time with a lower density, in a single layer and with more similar proportions between the bone and lithic remains.
As for the single-material accumulations, in the inner half of the rock shelter, the most evident case is the accumulation of bones in the central-northern zone, which coincides with the NE end of the cKM4 cluster (Figs. 9 and 14 profile A, Fig. 15 profile K). Many of these bones are burnt, as evidenced in previous work (Gabucio et al. 2014b;Chacón et al. 2015). The other accumulations of one specific material were identified in the external part of the site, where the density of remains is lower. For instance, a cluster of lithic remains was identified in the central area of KM1 (Fig. 15, profiles G, H, and I). Somewhat more towards the theoretical S, also in KM1, is a faunal accumulation (Fig. 15, profile I). Likewise, in the SE quadrant of KM4, there is a very evident preponderance of faunal remains over lithics. In addition, in profiles H, I, J, and K (Fig. 15), it is possible to appreciate a difference in the proportions of materials according to microlevels in KM4 and part of KM3, while in Ob1, lithic remains dominate over faunal ones, in Ob2 faunal remains dominate over lithics. This fact had already been pointed out in the N zone (extreme NW of sector KM4) in Chacón et al. (2015), but this is the first time the phenomenon has been identified in the NW zone (sector KM3).
Before delving into the interpretation of these different accumulations (both the mixed ones and those that mainly present either lithic or faunal remains), the site formation processes must be addressed.

Site formation processes
Earlier publications recognised Neandeprthals as being the main accumulators and modifiers in archaeolevel Ob (Gabucio 2014; Gabucio et al. 2012Gabucio et al. , 2014aGabucio et al. , b, 2018a. Likewise, previous taphonomic, spatial, and refit analyses evidenced that, although the preservation of the remains and their spatial properties are generally very good, some postdepositional processes have altered the assemblage, occasionally causing some local movement (Gabucio 2014;Gabucio et al. 2012Gabucio et al. , 2018a. This study has allowed us to clarify and reinforce the interpretation of postdepositional processes and their effect on the spatial distribution. The predominance of remains < 3 cm in length, which represent more than 70% of the coordinated lithic and faunal remains, in all sectors and clusters, corroborates the absence of a generalised transport phenomenon that would have affected a large part of the assemblage. The fact that small-scale remains are proportionally more abundant in clusters cKM1, cKM3, and cKM4, located inside the shelter, when the slope of the rock shelter descends from the flat area of the interior to the exterior, reinforces this idea, ruling out general water transport. Furthermore, the fabric analysis shows similar proportions of the four orientations (NS, NESW, EW, and NWSE) and a fairly even distribution of these in the rock shelter, confirming this conclusion. The slight predominance of NS and EW orientations could be related to a natural tendency for excavators to indicate NS and EW orientations more frequently than NESW and NWSE, as already suggested by Romagnoli and Vaquero (2016) for the lithic assemblage from level M at the same site. All this evidence ruling out the generalised transport of archaeological materials is added to other data from previous analyses, such as the presence of most of the calcined fragments inside the combustion structures, the abundance of lithic and faunal refitting groups, and the predominance of very short connection lines (Bargalló 2014;Gabucio 2014;Bargalló et al. 2016;Gabucio et al. 2018a).
However, as we noted before, there is data that reveals the effects of local and short postdepositional movements. Previous publications analysing the taphonomy of the macrofaunal assemblage (Gabucio 2014;Gabucio et al. 2018a) have evidenced the existence of a few items (68, 0.79%)

Fig. 14 Profile projections oriented in the EW direction (A-E). Continuous lines indicate the location of the sectors (KM1-4), and discontinuous lines show the location of the clusters (cKM1-4)
◂ that show a high degree of rounding on all their surfaces (R3, following the methodology proposed by Cáceres 1995), suggesting active role in the abrasion process and therefore local removal (displacement of items accumulated over the substrate). In addition, the location of a few calcined remains a certain distance from the hearths (following the natural slope), and some refitted dry broken bones separated by medium distances have been interpreted as the result of local and isolated movements on the surface caused by occasional water flows or gravitational movements. Finally, the very few refitted bones connecting microlevels Ob1 and Ob2 only in the NW corner of the rock shelter indicate that not only horizontal but also vertical postdepositional movements took place in this area. In this regard, a local taphonomic phenomenon, interpreted as a pond that would have filled and dried cyclically, was identified in this area (Gabucio 2014;Gabucio et al. 2018a).
This study has added further evidence in this regard. The main change in slope detected in the vertical profiles (Fig. 15) and the DEM (Fig. 16) from the flat inner area to the lower outer area can serve as a reference to locate the cornice limit (differentiating an interior space from an exterior space). In turn, the cornice limit is useful for identifying other related phenomena, such as drip zones and scour surfaces, which may have played an important role in the postdepositional mobilisation of archaeological remains. For instance, the depressed area visible in the W of Fig. 16C and profile G (Fig. 15) could be interpreted in this way. Likewise, the fact (easier to appreciate in the lithic material) that in the outermost fringe of the rock shelter, located downslope, there are larger remains than in the inner area, located upslope, could be related to gravitational transport of these heavier objects (Figs. 7, 8, and 9). However, we cannot rule out other explanations for this pattern, including the lack of small remains in the outer area of the rock shelter.
In the NW corner, the horizontal distribution of the remains in the form of an ellipse with many remains on the perimeter and fewer remains inside has been complemented with vertical projections, which reveal a cuvette-shaped arrangement of the remains in this area. Furthermore, the projections show that, inside the cuvette, the KM2 sector is generally in a slightly higher position than the KM1 and KM3 sectors. The digital elevation model complements the delimitation of this area: its superimposition with the location of the clusters (Fig. 16B) shows that the dynamics of the pond fully affected cKM2 and partially affected cKM1 and cKM3. The combination of all this data and the scarcity of remains ≤ 1 cm in length in KM2 and cKM2 might suggest that the humidification and drying cycles of the pond would have caused the horizontal displacement of the smaller and lighter archaeological remains from the KM2 sector towards the nearby areas of sectors KM1 and KM3. However, as many remains ≤ 1 cm were not seen by archaeologists during fieldwork and, therefore, were not piece-plotted (especially in the case of the remains recovered from the year 2011, due to the change in the cut-off size for piece-plotting), to verify this possibility, it is necessary to review the remains stored in the bags of general findings and wet sieving. With this aim, we tested the content of these bags for a selection of squares. Specifically, we selected squares R60, R61 (both located entirely in KM1, practically all its surface area being included in cKM1), T58, T59 (both located in KM3, largely included in cKM3), T60 and T61 (both located entirely within KM2, practically all its surface area being included in cKM2). As the lithic assemblage is still under study, only the faunal remains were taken into account for this test. The results indicate that bags of general findings and wet sieving from T60 and T61 (in KM2) contain a lot of remains that are ≤ 1 cm in length, even proportionally more than the rest of the tested squares. In fact, if we add the uncoordinated faunal remains to the coordinated ones from these two grid squares and calculate the percentages of remains by length categories, we can see that the items ≤ 1 cm represent 61.57% of the remains recovered from these two squares, while in the rest of the grid squares (related to KM1 or KM3), the percentage of this category is around 50%. Consequently, the quantity and ratio of items ≤ 1 cm in the tested squares do not support the theory of the smallest debris from the KM2 sector being displaced to the KM1 and KM3 sectors.
When the analysis of the lithic material is complete (the flint remains in the NW zone, which are the most abundant, are currently being analysed), we will be able to address the taphonomic data obtained from the lithic and faunal analyses according to the approach developed in this work: using density maps, tables quantifying alterations by sector and cluster, and visualisating the vertical distribution of the remains according to their modifications. This proxy will allow us to further refine our understanding of the local postdepositional processes and their effects on the spatial distribution of archaeological materials. A first exploration of the results from the faunal analysis has allowed us to see, as an example, that the very rounded bones (R3 on all their surfaces) are more abundant in sector KM2 (1.86%), followed by sectors KM1 (1.03%), KM3 (0.84%), and KM4 (0.38%). In KM2 and KM1, most of these remains are located within the cuvette zone, although some of those in KM1 are found in the depressed zone at the W end or in the SW corner (related to the formation of a drip zone and scour surfaces), the two zones with the lowest Z-values. In KM3, some remains fall within the cuvette while others appear at the border with the KM4 sector, which presents quite low Z-values. Finally, the bones with R3 rounding from KM4 appear fairly uniformly distributed in the N zone of this sector, being totally absent from the E half (columns 39-45). All this suggests that the high degree of rounding is not related to long-distance transport, but rather to repetitive processes that involved the action of water, sediment, and archaeological remains in a reduced space, such as the pond, the drip zones, and short-distance scour surfaces. This evidence is in line with the results obtained from other levels of the site, such as the absence of size sorting between abraded bones observed in level J .

Neanderthal space management in archaeolevel Ob
Once it had been confirmed that the spatial distribution of the materials mainly reflects the activity of the Neanderthals (but keeping in mind the short-distance postdepositional movements detected), it was time to interpret the anthropic use of the rock shelter during the formation of archaeolevel Ob. To do this, we combined the results presented here with other previously published data: location and features of the combustion structures (Vallverdú et al. 2012a, b;Vallverdú 2018;Borrell 2018;Gabucio et al. 2018b); partial technological analysis (Bargalló, 2014); zooarchaeological analysis (Gabucio 2014;Gabucio et al. 2014aGabucio et al. , b, 2018a; and refits (Bargalló, 2014;Gabucio 2014;Gabucio et al. 2018a, b). Three figures in the Supplementary Information (SI Figs. 9-11) illustrate and summarise this subsection. The three clusters with the highest density of material (cKM1-3) are located in the NW corner of the rock shelter and coincide totally or partially with the combustion structures AR11-2, AR11-5, AR07-08-11/1 (cKM1), AR11-4, AR11-11, AR11-13 (cKM2), AR10-2, AR08-10-11/1, AR10-4, and AR10-3 (cKM3) (Fig. 3). Several of these hearths are small (Fig. 3), although two of them, AR07-08-11/1 (Ob1) and AR08-10-11/1 (Ob2), located a couple of metres from the wall, are medium-sized (Borrell 2018). The abundance of lithic debris in both microlevels, mainly around medium-sized hearths, indicates that knapping activities were frequent in this zone (Bargalló 2014). Likewise, in the same spaces, a high number of small faunal remains of different taxa show evidence of carcass processing and consumption (cut marks, low burning degrees, and percussion products) (Gabucio 2014;Gabucio et al. 2018a, b). Consequently, according to ethnoarchaeological data (e.g., Yellen 1977a, b;Binford 1978b;Hayden 1979;Kent 1987;O'Connell 1987;Gifford-Gonzalez 1989, Kroll andPrice 1991;Gamble and Boismier 1991;O'Connell et al. 1991), these accumulations have been interpreted primarily as the product of repeated domestic activities, although some peculiarities have been observed in KM2 (notably the percentage of lithic cores and faunal percussion products) (Bargalló 2014;Gabucio 2014;Gabucio et al. 2018a, b). This study, by evidencing the mixed nature of these accumulations around hearths in both the Ob1 and Ob2 microlevels (except at the E end of cKM3, where the microlevels show differential proportions of the materials), seems to support the interpretation of cKM1-3 as a recurrent domestic activity area, where technological activities and carcass processing and consumption might have occurred in an interspersed and continuous manner. The predominance of lithic remains over faunal ones suggests various possibilities: this domestic space was used more for technological than food purposes; this space oscillated between a domestic area and a technical area (which we cannot separate archaeostratigrahically, perhaps due to postdepositional processes); or there was more systematic cleaning of the bone refuse (for instance, accumulating and calcining this at the NW end of cKM4). However, the complexity of the postdepositional processes that occurred in this zone makes a more detailed interpretation difficult.
A similar mixed distribution of materials, although less dense and with more equal proportions between the different materials (lithic/faunal), is observed in the NE of KM4 (one part coinciding with the E half of cKM4 and another part extending towards the NE corner of the rock shelter). The remains are spatially related to several small and medium-sized combustion structures: AR07-1, AR06-1, AR07-8, XIIb-b (cKM4), AR07-3, AR09-2, AR07-6, AR07-2, AR08-2, AR07-5, and AR07-7 (Vallverdú et al. 2012a, b;Vallverdú 2018;Borrell 2018). As in cKM1-3, technological and zooarchaeological analyses point to a recurrent use of this zone as a domestic activity area where Neanderthals carried out their daily tasks (Bargalló 2014;Gabucio 2014;Gabucio et al. 2018a, b). Once again, this work reinforces this interpretation by demonstrating that the lithic and faunal remains present a joint distribution, this time in a single and less dense layer.
Between the mixed accumulations of material, interpreted mainly as domestic activity areas (cKM1-3 and E half of cKM4), and the wall of the rock shelter, there are some combustion structures that present several characteristic features of the hearths found in resting and sleeping areas: they are small, suitable for providing light and heat, separated from one another by approximately 1 m, located near the wall, and associated with scant archaeological material (Binford 1983;Vallverdú et al. 2010Vallverdú et al. , 2012aVallverdú 2018;Borrell 2018;Gabucio et al. 2018b;Spagnolo et al. 2019).
The clearest examples of these combustion structures are AR11-1 (KM1, Ob2), AR11-6 (KM1, Ob1), AR11-8, AR11-9 (KM2, Ob1), AR11-7 (KM3), AR09-1, and AR08-1 (KM4) (Fig. 3), although a geospatial study of the combustion structures using GIS methods did evidence other less visible sleeping areas (Borrell 2018). In microlevel Ob1, the aligned distribution and location close to the wall of these hearths are reminiscent of the combustion structures in the resting area identified in level N at the same site (Vallverdú et al. 2010). However, while level N is a highly visible archaeological assemblage resulting from short-term occupations, the palimpsest of level O is much more developed than that of level N. Possible resting and sleeping areas have also been proposed for level M and archaeolevel Oa (Gabucio et al. 2018b;Bargalló et al. 2020a, b). Other European and Near Eastern Middle Palaeolithic sites, mainly dated as MIS3, present similar evidence of possible resting and sleeping areas, for example, Oscurusciuto (Spagnolo et al. 2019(Spagnolo et al. , 2020b and Tor Faraj (Hayden 2012;Henry 2012).

3
The E half of cluster cKM4, occupying the inner central zone of the rock shelter, is more difficult to interpret when taking into account both the faunal and lithic remains. Considering only the technological analysis, this zone could also be considered a domestic activity area (Bargalló 2014). However, the faunal analysis indicates a more specialised use of this space. Firstly, the huge cluster of small bone fragments identified in cKM4 (Figs. 9 and 14 profile A, Fig. 15 profile K) actually corresponds to an accumulation of calcined bones recovered from inside a combustion structure (AR06-07-10-11/1) located primarily in microlevel Ob2 (although there are some calcined bones in microlevel Ob1). This accumulation has been interpreted in previous works as a possible post-hoc zone where skeletal remains would have been systematically burned to reduce the volume of waste and, possibly, take advantage of skeletal fat as a complementary fuel (Gabucio 2014;Gabucio et al. 2014bGabucio et al. , 2018a. Secondly, closer to the wall, the abundance of small bone and teeth fragments and the high number of items showing diagnostic elements of anthropogenic breakage (especially in Ob2) suggests that this zone was specifically used for fracturing bones to obtain the marrow (Gabucio et al. 2018a). The presence of an anvil and some hammerstones in the vicinity seems to support this interpretation (Bargalló 2014).
Nevertheless, the different lithic/fauna ratio identified during the vertical analysis of the two microlevels in part of cKM3 and cKM4 (Fig. 15, profiles J and K) explains the disparity when interpreting the functionality of this central inner area, since the main accumulations of lithics (in Ob1) and faunal remains (in Ob2) occurred at two different times during which this space would have been used for different purposes. Thus, during the formation of Ob2, this area could have been used for specialised tasks related to faunal resource management, while during the formation of Ob1, it could have been used mainly for lithic technological purposes or as a domestic area. Last but not least, the microstratigraphic analysis of the largest combustion structure in this zone revealed the existence of several fire-use areas separated by a little more than a metre and covered by heterogeneous carbonaceous beds, leading to its interpretation as a low-visibility sleeping and resting area (Vallverdú 2018;Borrell 2018). In short, this central inner zone would have been repeatedly occupied by Neanderthal groups but for diversified purposes, giving rise to a complex palimpsest of anthropogenic activities. Although the evidence of diachrony between these activities is numerous (existence of two microlevels, evidence of reuse of hearths, etc.), some synchronous relationships have also been identified thanks to refits  and to the archaeomagnetism of three combustion foci (Carrancho et al. 2016).
In the external part of the site, no relevant mixed accumulations have been found. In contrast, several accumulations clearly dominated either by lithics or faunal remains have been identified. An example is a cluster of lithic artefacts located in the central area of KM1 (Fig. 15, profiles G, H, and I), towards the NW side of combustion structure XVI. Several pieces of burned limestone were recovered, many of which are connected by refits. All these pieces have a high iron content, and are burned, fragmented, and incomplete, their interior part being missing. For all these reasons, it has been previously suggested that these remains could be associated with the production of a colouring product (Bargalló 2014). Also, in KM1, towards the S of combustion structure XVI, the faunal accumulation (Fig. 15, profile I) corresponds to an almost complete wildcat skeleton. This wildcat was processed and consumed by the Neanderthals at the same location (Gabucio 2014;Gabucio et al. 2014aGabucio et al. , 2018a. Likewise, in the SE quadrant of KM4, near combustion structures I, II, XVIII, and XX, some deer remains were refitted and identified as a single individual. Since the cranial skeleton of this individual is almost complete and only right elements of the postcranial skeleton have been recovered, it has been proposed that this animal could have been quartered in this area of the rock shelter and then some parts transported to other areas (Gabucio 2014;Gabucio et al. 2018a, b). However, the presence in this zone of KM4 of the only cuvette (I) and à event hearths (XVII and XX) (Vallverdú et al. 2012a) and the high percentage of bones burned to low degrees have also led to the suggestion that this area was used to carry out a particular cooking or food preservation technique (Gabucio et al. 2018a). Finally, another bone cluster in KM4, at the central E end (squares O-Q/40-41), was interpreted as a natural accumulation, comprising an abundance of rabbit bones and tooth marks, and a very low ratio of burned remains (Gabucio 2014;Bargalló et al. 2016;Gabucio et al. 2018a).
If we combine the sectors and clusters proposed in this study with the analysis of faunal refits (Gabucio 2014;Gabucio et al. 2018a, b) and the partial analysis (without the pending flint remains from the NW zone) of the lithic refits (Bargalló 2014; Bargalló et al. 2016), we obtain interesting, albeit provisional, information about the relationships between the different areas. Most refitted pieces are connected by short distances, although some long connections (up to 18 m) have been identified. The most noteworthy feature is the main orientation of the longest refits, which connect the theoretical E and W ends, horizontally crossing the inner half of the rock shelter, leaving the outer half without inter-sector connections (Bargalló 2014;Gabucio et al. 2018b). Both lithic and faunal inter-sector refits connect cKM1 with all the other sectors. The link between cKM1 and the NE corner of KM4 seems particularly significant, as both have been interpreted as domestic activity areas and they are related by bidirectional lithic refits and one faunal refit (mechanical green breakage). Bidirectional lithic refits also connect KM4 (NE corner) with all the other sectors. No faunal refits connect KM4 to either KM2 or KM3, but other criteria (preponderance of right or left aurochs elements) suggest possible food sharing between KM4 (NE corner) and KM2, and tooth microwear analysis suggests that horses from KM3 and KM2 could have died during a short period, for example, the same season (Gabucio et al. 2018a). The pattern of directionality observed in the lithic refits suggests a possible temporal order for the sectors: some concentrations of remains in KM4 would thus be older than others in KM2 and KM3, and these would in turn be older than others in KM1. However, faunal refits and other zooarchaeological criteria suggest food sharing over time between these sectors (Gabucio et al. 2018a, b). In any case, as the NW area was affected by a local postdepositional process, we should be cautious when interpreting the refits from this area, especially with the refits that connect the cKM1-3 clusters.

The Ob archaeolevel in the framework of the Abric Romaní site
Our long-term objective is to produce quantitative spatial data from several levels at Abric Romaní in order to be able to undertake a reliable diachronic comparison of the site. This task cannot yet be fully addressed since, prior to this work, studies of this type (including methods such as Ripley's K function, KDE, k-means, and similar) have only been implemented on the lithic assemblage of levels I, J, and M (Vaquero 1999;Vallverdú et al. 2005;Romagnoli and Vaquero 2016) and the faunal assemblage of level P (Marín et al. 2019). For level J, Sañudo et al. (2012) used density maps and size sorting to address the spatial patterns in archaeolevel Ja. In addition, Bargalló et al. (2020a, b) applied tools including palaeosurface reconstruction, density maps, and the k-means neighbour distance to archaeolevel Oa. Lastly, in all these levels, a vertical study of the distribution of the remains has been carried out. Thus, to date, according to the general features of the assemblages and the methods applied, levels M, P, and, to a lesser extent, archaeolevels Ja and Oa provide the best data to be compared with the results of this work.
However, before comparing the spatial patterns observed in the Ob archaeolevel and other layers at the site, it should be noted that during the formation of level O, there was a significant change in the configuration of the rock shelter. In the upper levels (from archaeolevel Oa upwards), the wall of the rock shelter closed towards the W side at the height of the stratigraphic testimonial, most of the remains located between this and the "Pou Romaní" pit (approximately, columns 45-54). In contrast, during the formation of archaeolevel Ob (and the lower units), the morphology of the rock shelter allowed the use of a further area towards the NW of the stratigraphic testimonial (squares R-W/54-62), not accessible in the upper levels. This meant that in the lower levels, there were larger flat areas in the interior than in the upper ones. Thus, at levels lower than Ob, such as P, Q, and R, the Neanderthal occupations were preferentially located on the flat surface towards the NW of the rock shelter. In turn, this fact conditioned the location of the combustion structures and activity areas, these being organised differently and displaced with respect to the patterns observed in the upper levels. It is also likely to have conditioned the main E-W orientation of long-connected refits in archaeolevel Ob.
Despite the morphological change in the rock shelter during the formation of level O, there are numerous similarities between the assemblages created before and after this restructuring. In levels M and P (horse remains), as in Ob, the remains show clustered spatial structures, tending to be concentrated in association with the combustion structures (Romagnoli and Vaquero 2016;Marín et al. 2019). Both in Oa and Ja, the distribution of the material is also linked to the location of the combustion structures (Sañudo et al. 2012;Bargalló 2014;Gabucio 2014;Bargalló et al. 2020a, b). In level M and archaeolevel Ob, the densest clusters are located on the flat inner areas (quite close to the wall, but not attached to it), although there are significant accumulations of material in the central area of level M. If we compare the average intensity of the inner accumulations of these two assemblages in terms of lithic remains (only the currently comparable data), Ob presents much higher values (up to 0.15018 in cKM2) than M (up to 0.01101 in M4), but we must remember that the main clusters in Ob are vertically subdivided into discontinuous microlevels, something that has not been observed in M (Romagnoli and Vaquero 2016). In addition, the study of Romagnoli and Vaquero (2016) does not consider any remains < 1 cm. Similarly, the remains concentrate in the inner zone in Oa (Bargalló 2014;Gabucio 2014;Bargalló et al. 2020a, b), although in Ja, the materials appear more frequently in the central area (Sañudo et al. 2012).
In both M and Ob, remains < 3 cm in length predominate in all accumulations (Romagnoli and Vaquero 2016;Vaquero et al. 2019). Small remains are also preponderant in Ja, where it has been observed, as in Ob, that the farther outward one proceeds, the fewer small remains are found. This characteristic of level J was initially interpreted as a possible divergence from the drop (interior) and toss (exterior) zones, but it was subsequently also linked to postdepositional phenomena related to the slope and the greater exposure of the outside area (Carbonell 2012). Nevertheless, several analyses rule out mass transport of remains in any of the studied levels at the site. As an example, the fabric analysis carried out in M and Ob, according to different clusters, rules out the presence of preferential orientations, and in the size analysis of several levels (including Ja, M, and Ob), no size sorting was identified (Sañudo et al. 2012;Romagnoli and Vaquero 2016;Vaquero et al. 2019).
In levels M and P, as in Ob, the results obtained have ruled out massive postdepositional movements of archaeological material, pointing to Neanderthals being the main cause of the spatial patterns. In level P (both archaeolevels, Pa and Pb), the different patterns of Cervus and Equus have led to the identification of two different uses of the rock shelter: as a hunting camp for horse exploitation; and as a shortterm transitory campsite related to the exploitation of deer (Marín et al. 2019). A differential pattern between different taxa was also noted in Ob, where horses appear to have been obtained over a short period of time, while aurochs, in contrast, evidence a long period of time (Gabucio et al. 2018a). In this sense, it seems that a similar trend is repeated at the two levels: horse remains would have been accumulated at the rock shelter over short periods of time during which the site would have been used as a hunting camp, while deer and/or aurochs would have been accumulated over longer periods of time while the site was being used for residential purposes (either in different separate short events, as appears to be the case for level P, or including some longer continuous periods of occupation, as might have occurred at Ob). As Marín et al. (2019) have already pointed out, this way of using of the same space for different activities is typical of groups of collectors with logistical mobility, who generally move radially within a territory (Mortensen 1972;Binford 1980;Kelly 1992).
In level M, differences in spatial patterns were observed between the inner and the outer zones of the rock shelter (Romagnoli and Vaquero 2016). Similar differences between the inside and outside were indicated from the study of hearth-related wood (Solé et al. 2013) and faunal remains (Fernández-Laso 2010; Gabucio et al. 2018b;Fernández-Laso et al. 2020) from the same level. However, when comparing the use of the inner and the outer areas between Ob and the upper levels at Abric Romaní, we must take into account the change in the morphology of the rock shelter, which probably altered the boundaries of the interior and exterior areas of the site. Thus, in Ob the innermost, most enclosed space is the NW area (cKM1-3), which is undoubtedly why it was occupied in a more recurrent way, generating a complex palimpsest with a very high density of remains (especially lithics). In addition, the rest of the N fringe of the surface area (cKM4) was protected by the cornice and can also be considered an inner area. Despite the change in the rock shelter configuration, the similarity of some parameters allows a preliminary comparison to be made between the inner/outer spatial patterns of archaeolevel Ob and level M.
On the one hand, the inner areas are characterised by an abundance of remains (concentrated in high-density clusters), a high number of interconnections (inter-accumulation refits) and spatial congruence (the same repeated use over time). In level M, these features were interpreted as the product of successive occupations of small groups that would have caused the development of a horizontal palimpsest (Romagnoli and Vaquero 2016). In archaeolevel Ob, the main clusters (cKM1-4) meet all these criteria. As already explained, accumulations cKM1-3 and the E half of cKM4 show a mix of faunal and lithic remains and have been interpreted as recurrent multifunctional domestic areas, although other occasional uses are not ruled out (Bargalló 2014;Gabucio 2014;Gabucio et al. 2018a, b). The W half of cKM4 could also have been used as a domestic space at times, but it was also the scene of more specialised repeated activities (Bargalló 2014;Gabucio 2014;Gabucio et al. 2018a, b). Finally, the fringe between the denser accumulations in Ob and the wall of the rock shelter have been interpreted as possible sleeping and resting areas (Vallverdú et al. 2012a, b;Gabucio et al. 2018b). Another possible resting area was proposed for level M, in the NE corner (Gabucio et al. 2018b), although its identification seems less clear. This bimodal spatial pattern (a sleeping area close to the wall and, next to it, a multipurpose area) has also been recognised in units SU13 and SU11 from the Oscurusciuto rock shelter (Spagnolo et al. 2019(Spagnolo et al. , 2020b and floors I-II of Tor Faraj (Hayden 2012;Henry 2012).
On the other hand, the outer areas present a low density of remains and a high visibility of specific events. In level M, the well-preserved isolated technical events and the recycling of several artefacts have led the researchers to propose short events, which probably occurred during the final phases of occupation Romagnoli & Vaquero 2016). In Ob, several individual episodes (following the nomenclature established by Gabucio et al. 2018a) were identified in external areas: possible colourant production from limestone; wildcat (towards the W); and red deer processing (towards the E, possibly related to the first phases of faunal chaîne opératoire or involving specialised cooking or food preservation techniques). These episodes have not been reliably linked to the inner areas of the rock shelter, so it is not possible to establish whether they occurred at times when the inner area was occupied or not. However, the vertical sections do not seem to support the idea that these accumulations were discarded during the last phases in which the site was occupied.
The differences between the inner and outer areas might be explained merely by a variability in intensity (related, for instance, to the number of occupations in each area, or to variations in the time span of the different events), or possibly by a differential use of space. Individual episodes are easily identifiable in the outer areas, where the intensity is low, because the preserved accumulations of materials are more isolated from each other. In contrast, the palimpsest nature of the inner areas, which were more intensely occupied, makes it difficult to identify specific events, which overlap with other associations accumulated at different times. In the case of Ob, this is especially evident in the NW zone, which, in addition to being the most intensely occupied area, suffered local horizontal and vertical postdepositional movements. This leads us to wonder if the mixed accumulations (with lithic and faunal remains), previously interpreted as domestic areas, were actually the result of several repeated mixed uses or, in contrast, of the overlapping of several accumulations of one or other material that, due to superimposition and/or postdepositional processes, we are now unable to delimit. Nevertheless, in the NE zone (E of the cKM4 cluster), which presents a lower density of remains and better preservation, a mixed accumulation was indeed identified, suggesting that the formation of this type of mixed accumulation could also have occurred in the NW zone. Similarly, in other levels at Abric Romaní with a higher temporal resolution, such as H, I, K, and L, the hearth-related clusters almost always include both lithic and faunal remains Vallverdú et al. 2005). In addition, despite the postdepositional movement, three high-density accumulations were observed in cKM1-3, all showing mixed lithic and faunal materials, seemingly reinforcing this same line of reasoning. Finally, the content of these clusters (small knapping products, bone percussion products, small burned bones, etc.) fits with the repeated use of these zones as domestic areas (Bargalló 2014;Gabucio 2014;Gabucio et al. 2018a, b). Undoubtedly, superimposition processes and postdepositional movements have interfered to some extent in the composition of the clusters, especially in cKM1-3, but they do not appear to be the underlying cause of the mixed accumulations of archaeological material.
In contrast, no mixed accumulations have been identified in the outer area of the rock shelter. In general terms, faunal remains are proportionally more abundant than lithics, contrary to that seen inside the rock shelter (especially in cKM1-3). Here, several accumulations of either one or the other material have been related to specific activities that are both spatially and temporally well-delimited, such as the processing and consumption of a wildcat, the production of a colouring material, and the separate processing of half a deer. Although it is evident that the lower density of remains in the outer area facilitates the identification of these events, we believe that their uniqueness (the only remains of Felis silvestris in the assemblage, the refitting between several burned limestone fragments, the refitting between deer remains, and the use of the only cuvette and à event hearths identified in the level) would have made them recognisable, to some degree, if they had also been present in the inner area of the rock shelter (as was the case with the accumulation of calcined bones, for example). In addition, it must be taken into account that the external area, despite the advantage of less intense occupation, presents other conservation problems, including the fact that it is less protected by the cornice of the rock shelter and, therefore, more at the mercy of inclement weather. Finally, the separation into different areas of domestic activities (preferential zones) from specific ones (marginal zones) has been ethnoarchaeologically documented (Yellen 1977b;O'Connell 1987). For all these reasons, we consider that, at least in the case of archaeolevel Ob, the data supports the idea that Neanderthals differentiated the use of space in the rock shelter between the inner (bimodal use, containing domestic and sleeping areas) and outer areas (more occasional and specialised used, particularly linked to faunal resources).

Conclusions
This work presents the first application of GIS and geostatistical methods to the study of the intra-site spatial distribution of lithic and faunal remains in archaeolevel Ob at the Abric Romaní site. As such, it fits into a proxy recently initiated by the Abric Romaní team, which, developing research lines such as palimpsest dissection and spatial taphonomy, is committed to a more quantitative and objective approach to the study of spatial patterns. This study is also novel in the fact that it is the first work of this type, in the context of Abric Romaní, to include both faunal and lithic remains. The analysis of the data, first combining the lithic and faunal remains and then separating them into two distinct assemblages has allowed us to characterise the Ob archaeolevel more accurately and facilitate the integration of the two materials in the interpretation. In addition, by including the results on the archaeolevel from previous work (technological, zooarchaeological, taphonomic, refitting, and combustion structure analyses), we have been able to employ a more comprehensive approach.
After confirming that the spatial structure of the assemblages is clustered, we explored the intensity of the accumulations using KDE and k-means analyses. In this case, despite its limitations, k-means analysis has been very useful as, in addition to determining the main accumulations (cKM1-4), it has allowed us to classify all the remains in different sectors (KM1--4). Thus, by using a parallel classification of the remains in sectors KM1-4 (including all the remains) and the accumulations in cKM1-4 (only the remains grouped into the main concentrations), we have been able to make a quantitative and in-depth study of the features of each subset, obtaining valuable information on the distribution of different variables, including dimensions, orientations, and materials in the rock shelter.
All these contributions expand the volume of information that we have for Ob, providing quantitative and objective data and favouring a more precise understanding of the deposit, allowing us to reach relevant conclusions about the formation and postdepositional processes. The predominance of remains < 3 cm in length and the absence of preferential orientations in all sectors and clusters (both lithics and bones) make it possible to rule out generalised postdepositional movement, leading to the conclusion that the observed spatial pattern mainly reflects the behaviour of the Neanderthal groups. However, some local postdepositional movement did take place. This work provides data that allow us to locate and interpret this with greater precision. For instance, the edge of the ledge has been located and associated with drip zones and scour surfaces. In addition, the area that would have been affected by the formation of an intermittent pond has also been more precisely delimited (cKM2 and the neighbouring areas of cKM1 and cKM3).
The horizontal and especially the vertical analyses have allowed us to identify similarities and differences between the distribution of the two types of archaeological materials. The densest accumulations of both materials occur in the inner area, especially in the NW (cKM1-3) and N (cKM4) zones. In cKM1, cKM2, and the E area of cKM4, the accumulations present a fairly homogeneous mixture of lithic and faunal remains. In contrast, in cKM3 and the W area of cKM4, there is vertical differentiation in the distribution of materials: in the upper microlevel (Ob1), lithic remains predominate over faunal ones while in the lower one (Ob2), faunal remains predominate over lithics. In the outer zone, which contains fewer remains, clusters mainly of lithic or faunal remains have been identified, although, in general, there is a higher proportion of faunal remains (especially in the SE area of KM4).
The combination of this data with the results of previous studies (technological, zooarchaeological, refit, and combustion structure analyses) has made it possible to interpret the accumulations with an anthropogenic slant, evaluating and uniting previous interpretations suggested by different approaches. In this case, it was necessary to go beyond the limitations of density and k-means analyses, differentiating discrete areas within the largest sectors (KM1 and KM4) and using other hierarchical elements of space, such as hearths, as a guide. Thus, the mixed accumulations identified in cKM1-3 and the NW corner of KM4, rich in small knapping products, bone percussion products, and small burned bones, have been interpreted as the result of repeated domestic activity, both in Ob1 and Ob2. Likewise, it is proposed that the fringe between these mixed accumulations and the wall, presenting few remains and small hearths, could have been used as sleeping and resting areas. The medial central area (W of cKM4) had different uses during the formation of the two microlevels, with a fauna-related specialised use in Ob2, inside the AR06-07-10-11/1 combustion structure (a possible posthoc area). In the outer half of the rock shelter, it is possible to identify different individual episodes: the processing of a block of limestone probably to produce a colouring product (central area of KM1), the processing and consumption of a wildcat (S area of KM1), and the processing of a deer (SE corner of KM4). Finally, a natural accumulation of rabbit bones was identified at the central E end of KM4.
The long-term objective of this approach is to optimise the diachronic study of the site and facilitate its comparison with other sites. In this regard, it has been possible to identify some similarities between archaeolevel Ob and other levels, especially levels M and P. The data obtained during the study of these three levels enabled us to rule out massive postdepositional displacement of the archaeological remains, the spatial patterns being mainly related to anthropogenic activities. Based on this premise, different uses of the rock shelter have been identified: as a hunting camp (horse exploitation) in level P and possibly Ob; as a short-term campsite in M and P (deer exploitation); and possibly as a long-term campsite in Ob (auroch exploitation, perhaps also deer). Likewise, in the M and Ob levels, different occupational patterns have been observed between the inner zones and the exterior areas. In the case of Ob, the presence of mixed accumulations only in the inner area (both in dense clusters and in other less dense and better preserved ones, in all cases containing large quantities of small knapping products, bone percussion products, and small burned bones) and other criteria (for instance, an increase in the proportion of faunal remains with respect to lithics and the identification of specific and temporally well-delimited activities in the outer zone) allows us to propose, for archaeolevel Ob, a differential use of the inner (more domestic) and outer (more marginal) areas.
Finally, we want to emphasise the convenience, as far as is possible, of approaching the study of archaeological assemblages from a high-resolution spatiotemporal perspective. This approach must be transdisciplinary, with a special focus on taphonomy. It should also be multiscalar and include the use of a wide range of tools (GIS, archaeostratigraphy, different geostatistical methods, refits, etc.), especially those that can provide objective and quantitative data to optimise the interpretation of the studied assemblages and facilitate their later comparison. writing. A.B.: conceptualization, lithic analysis, spatial analysis, investigation and data curation, revision of the manuscript. P.S.: geostatistics, investigation and data curation, revision of the manuscript. F.R.: investigation and data curation, revision of the manuscript. M.G.C.: investigation and data curation, revision of the manuscript. J.V. combustion structures analysis, investigation and data curation, revision of the manuscript. M.V. investigation and data curation, revision of the manuscript.
Funding Open access funding provided by Universitat Rovira i Virgili. Funding for fieldwork at Abric Romaní site is provided by the Ajuntament de Capellades and Romanyà-Valls. The Catalan Government (Generalitat de Catalunya) supported this research with the Quadrennial Project CLT009/18/00054, the project "Territoris prehistòrics de la Conca de l'Anoia (2022-2025)" (exp. ARQ001SOL-201-2022)  Data availability All data generated or analysed during this study are included in this published article (and its supplementary information files).
Code availability Not applicable.

Declarations
Ethics approval Not applicable.

Consent to participate
All authors consent to the participation.

Consent for publication
All authors consent to the publication.

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