The Spatial Structure of Galician Megalithic Landscapes (NW Iberia): A Case Study from the Monte Penide Region Authors

Abstract It is well known that Neolithic megalithic landscapes are the result of complex locational logics governing where communities chose to site their funerary monuments. These logics in turn respond to broader environmental and cultural affordances, and the relationship between these has been a major topic in the megalithic archaeological literature for the last few decades. Thanks to new approaches in spatial statistical modelling, there is now considerable opportunity to revisit traditional megalithic locational concepts from a more systematic point of view, not least in Galician studies (NW Iberian Peninsula). In the paper that follows, we apply such a modelling approach to a large set of megalithic monuments located in the south of Galicia (Monte Penide and surroundings) with a view to exploring locational choices, spatial hierarchy and territoriality in these funerary landscapes. The results indicate that the distribution of megalithic mounds in this region reflects a preference for locations with particular environmental properties, while at a more local scale the spacing of these mounds seems to reflect some kind of social partitioning of the landscape. Via spatial cluster analysis and a further novel method for testing site hierarchy, we conclude that the mound sizes within nine different mound clusters exhibits a non-random hierarchical structure, with a larger mound per group and smaller ones around that, and with what appears to be a preference for the large monument to be at or near the meeting point of several watersheds and upland ridge-routes.


Introduction
Megalithic tombs are a common social and funerary feature of Atlantic European landscapes during the Neolithic, involving the construction of large-scale stone monuments, typically but not exclusively collective graves covered by a mound (Laporte, Scarre 2011). Many explanations for the rise of megalithism and for the meaning of particular sets of megaliths have been offered by archaeologists in the past (Schultz Paulsson 2017), with an emphasis generally placed on their multiple rather than singular roles in society, and their likely symbolic marking of communally-shared resource landscapes belonging to segmentary groups (e.g. Saxe 1970;Fleming 1973;Chapman 1981;Tilley 1994;Sherratt 1990Sherratt , 1995. Such investigations therefore highlight the importance of megalithic monuments not only as a funerary monuments but also as visible, public landmarks for living communities (Renfrew 1976;Renfrew 1984;Delibes de Castro 1991). Their potential territorial significance has been approached from different perspectives, including theoretical landscape archaeology (e.g. Criado Boado, Villoch Vázquez 2000;, historical comparanda (e.g. Martinón-Torres 2001;Díaz Guardamino et al. 2015) and quantitative spatial analysis (e.g. Gillings 2009;Murrieta-Flores 2012, 2014Llobera 2015).
The development of the Neolithic Period and the megalithic complex of northwest Iberia probably reflects funerary ideas that originally came from the centre of Portugal (Rodríguez Casal 1990;Prieto Martínez et al. 2012;Schultz Paulsson 2017). Although as yet constituting only a small number of finds, and in contrast to the pattern suggested elsewhere in Iberia, the existing archaeological evidence in the north-west suggests a degree of continuity between the latest Epipalaeolithic hunter-gatherers and the first food-producing communities, leading to the development a mixed economy during the first half of the 5th millennium BC. The earliest radiocarbon date associated with archaeological structures from this initial period is 4720-4530 BC (on an organic-rich sediment; Ua-3267, 5780 ± 40 uncal bp), and comes from a storage pit at Monte dos Remedios (Fábregas Valcarce et al. 2007;Prieto Martínez et al. 2012: 220). Additional available radiocarbon dates suggest, in general dates, the appearance of simpler singlechambered Galician megaliths over the period 4500-3500 BC, a subsequent development of corridor dolmens (3500-2500 BC) and then a late phase of megalithism (2500-1800 BC) characterised by small cists or burials without above-ground structures.
Despite a strong tradition of scientific study of these monuments by Galician and other researchers (see e.g. Criado Boado 1989;Bradley 1991aBradley , 1991bBradley , 1998, spatial approaches have until recently rarely extended beyond the creation of distribution maps (Leisner 1938;Rodríguez Casal 1997), accompanied by, often sensible but informallyexpressed, comments about locational factors such as those outlined in Table 1 In the last few years, however, researchers have begun to adopt a more quantitative approach to the same observations via GIS and spatial modelling (Wheatley et al. 2010), although the emphasis has so far been on the relationship between sites and potential travel networks across the landscape (see e.g. Llobera 2015; Rodríguez Rellán, Fábregas Valcarce 2015; Carrero-Pazos 2018a; Carrero-Pazos, Rodríguez Casal in press). This growing body of work has identified a series of potentially relevant environmental Final version: Carrero-Pazos, M., Bevan, A., Lake, M. 2019. The spatial structure of Galician megalithic landscapes (NW iberia): A case study from the Monte Penide region. Journal of Archaeological Science 108. https://doi.org/10. 1016/j.jas.2019.05.004 variables (such as topographic prominence or natural transit routes) that exhibit considerable statistical interdependence with one another (i.e. are highly correlated).
Given this interdependence, many of these variables could, interpretatively, be seen as different ways of expressing the same general observation. In light of this, and after much wider exploration of the possibilities (Carrero-Pazos 2017) we have opted to consider just two of the simplest but most important variables -elevation and distance to major watershed boundaries (both variables were discussed by Bradley 1991a: 78;1991b in the case of south-west England burial mounds; by major watersheds here we generally mean those that have maritime outlets). This simpler starting point then allows us to conduct a more rigorous analysis and investigate both first and second-order patterns and site sizes (for the definition of these terms, see below and Bevan et al. 2013;Baddeley et al. 2016). We conclude that these Galician monuments made quite explicit use of the natural features of the landscape such as the ridgelines high up in major Galician upland watersheds which were both important transit routes for travel and major viewpoints. Moreover, we demonstrate that these monuments were organised into spatial clusters whose internal size hierarchies exhibited non-random structure suggesting some kind of organisational equivalence between clusters. These results can be independently reproduced with the data and analytical workflow provided in the accompanying repository (see Supplementary Data and Methods).

Study area
The study area is located in southern Galicia (NW Iberia, Figure 1  In geographical terms, this is a comparatively flat region rising from the coast to moderate elevations (400-500 mASL). From a hydrological point of view, the watersheds of Verdugo and Oitavén rivers stand out, as well as the Miñor watershed which bounds the region to the south. Pico San Vicente (432 m) in Redondela, is the highest point in the north, and has one of the largest concentrations of mounds, such as that at Monte Penide. The main mountain range of the area is further to the East, the elongated (N-S direction) Serra do Galiñeiro (Coto de Cales, 742 m).
Nevertheless, despite a long history of research, the Galician megalithic phenomenon has only very recently been studied from a spatial perspective (Carrero-Pazos 2017, B: Field plans of the excavations carried out by C. Mergelina (1936) in some of the monuments.

Data and methods
The archaeological data used in the analysis reported here is part of an archaeological database produced by the megalithic studies group at the University of Santiago de Compostela (GEPN-AAT). This group has catalogued 121 sites via both the accumulation of previous fieldwork (Carrero-Pazos 2017).
Site location analysis is one of the most prolific applications of GIS in archaeology over the last 35 years (early examples include: Kvamme 1983Kvamme , 1984Judge and Sebastian, 1988). Although these approaches have been criticised as encouraging undue Final version: Carrero-Pazos, M., Bevan, A., Lake, M. 2019. The spatial structure of Galician megalithic landscapes (NW iberia): A case study from the Monte Penide region. Journal of Archaeological Science 108. https://doi.org/10.1016/j.jas.2019.05.004 environmental determinism (Wheatley 1993(Wheatley , 1996(Wheatley , 2004Gaffney, van Leusen 1995), they have nevertheless remained popular, especially when allied with modern spatial statistics which, via their ability to discriminate between first and second order effects (Bevan et al. 2013;Baddeley et al. 2016), potentially offer a more nuanced account of the balance between environmental and social causes of past human behaviour (see e. g.  Table 2, which shows which covariates appeared to be better or worse predictors of the location of megalithic tombs in the study area. [Caption] Although previous work (Carrero Pazos 2017, 2018b) identified four covariates as important variables, an elevation cut-off is perhaps the most obvious one mentioned in previous discussion of megalithic monuments (Bradley 1991a: 78;1991b). As noted above, here we deliberately retain a certain simplicity by investigating the role that just two major variables played in this megalithic area: elevation and distance to watershed edges (Figure 3: A, B), not least because these two are also likely to be useful proxies for other meaningful locational priorities such as ridge-routes and high visibility locations. We first conduct a univariate analysis of how the intensity of megalithic sites varies with respect to these two key variables and then build a multivariate logistic model which can account for the first-order spatial inhomogeneity in the megalithic sites, in other words, the overall spatial trend explained by properties of the environment. We then use this model as a control which allows us to identify the second-order properties of the site pattern, that is, those which may be more explicable in terms of social processes, such as what appears to be a tendency for a pre-existing mound to encourage the construction of nearby ones. In the final part of the paper, we use the second-order results and a mark correlation function to suggest a meaningful clustering threshold with which to identify sub-groups of mounds using a dbscan method (Ester et al. 1996; as implemented in GRASS GIS v.cluster). Having identified these groups, we then use a novel rank permutation method to conclude that the tombs are distributed across them in a way that is that is broadly (even if not perfectly) hierarchical to an extent that is unlikely to occur by chance alone.

First-Order Location Model
The general distribution map of megalithic monuments documents a non-random distribution of sites over the study area (see figure 1) and we can use statistical methods to estimate the first-order trend behind much of this variation (Baddeley et al. 2016).
First, it is useful to consider a non-parametric summary of the univariate relationships between the dependent variable (presence/absence of megalithic sites) and each of the two covariates (elevation and watershed edges, Figure 3C-D). These suggest that sites are more likely to occur at elevations ranging from 400-500 m ASL and at distances closer to watersheds than expected by chance alone. Comparison with previous work confirms that these two variables explain most of the overall variation even when further possible covariates are included (see Carrero-Pazos 2018a). We therefore build a first order logistic model using both covariates and produce a prediction surface of megalithic site intensity across the landscape. This model is then used in a second stage to study the clustering trend of the sites by fitting known statistical models to their spatial distribution ( Figure 3E). We have shown that elevation and watershed distance go some way towards explaining spatial patterning in the distribution of sites, but we can also investigate the possibility that there are additional -perhaps more overtly social -causes. We do this by using the logistic model of first order spatial inhomogeneity to detrend for the unevenness of megalithic sites due to landscape preferences, which then allows us to detect whether there is additional second-order clustering or regularity in the site point pattern across multiple scales (Baddeley et al. 2000: 330;Palmisano 2012;Baddeley et al. 2016). The inferential logic of this enquiry involves two main steps: first, we use a pair correlation function -which summarises the typical point intensity found in a series of noncumulative buffers around each point in the dataset -to identify second order clustering and construct a 95% critical envelope for this function via random simulations conditioned on the first-order regression model (Baddeley et al. 2016: 225-230). We then consider whether the residual clustering can be explained with reference to one of the well-defined second order interaction models in the literature.

Second-Order Clustering of Mounds
A basic pair correlation function confirms that the points are spatially clustered at distances up to 1km (the black line remaining above the grey envelope in figure 4A-B).
More precisely, the observed function (in black) falls above the Monte Carlo critical envelope at these short distances, before dropping back down within the expected range after this. A second version, in which the envelope is conditioned on the first order model ( figure 4B), demonstrates that although first order trends are important, they do not account for the overall clustering (unless a major 'lurking' first-order variable has been ignored in both this and previous work, which we view as unlikely), and that the clustering is more likely to be endogenous and second-order in character (with the presence of one site making a nearby site more likely; see also Baddeley et al. 2016: 487). It is also interesting to note that megaliths possibly exhibit a slight dispersed pattern at distances of ca. 2000 meters (where the black line falls just beyond or at the lower limit of the grey envelope in figure 4A-B).
[Caption] Figure 4. Approaches to site distance and confidence. A: pair correlation function with a 95% envelope from wholly random Poisson process. B: pair correlation function of the observed sites with a 95% envelope conditioned on the first-order covariates model. C: pair correlation function with a 95% envelope also conditioned on both the first-order covariates and a second-order, area-interaction model (r=1500 m). Having identified clustering, we then fit a known point interaction model to the observed pattern. There are several Gibbs-type point processes that can be fitted: examples include a so-called hard core process, in which points strictly avoid each other up to a certain threshold, and a Strauss process in which points have constant influence within a certain distance threshold (Nakoinz, Knitter 2016: 131). However, in our case, the Widom-Rowlinson penetrable sphere model or area-interaction process (Baddeley, Lieshout 1995;Baddeley et al. 2016: 519) offers a better fit and more interpretative salience. This model generates inhibition and clustering patterns with reference to a buffer created for all the points of the distribution, which may be interpreted as modelling a scenario in which monuments have an area of influence within which they either attract or repel other monuments. Figure 4C shows that, after adding such a second-order interaction to our first-order model, the observed pair correlation function (black line) now falls within the critical envelope of the simulation at all interaction distances. This suggests that the jointly fitted first and second order factors are now able to account for the observed point pattern. In other words, the distribution of megalithic tombs can be accounted for by broad locational preferences for elevation and watershed boundaries, as described by the regression model in tandem with local tendencies towards mound clustering.

Mound Size and Shape
Given the observed second-order clustering of megalithic sites, it is interesting to further explore the nature of this clustering and what it might imply in terms of social organisation. If, for example, we use a mark correlation function to consider the spatial correlation of mound volumes (referred to hereafter as 'sizes'), there is evidence for significant auto-correlation of these sizes, not for mounds in very close proximity but rather for mounds spaced about 4.5km apart (Figure 5: A). Put very crudely, it would appear that whatever process dictates tomb size repeats at approximately 4.5km intervals, which further implies that it might be possible to detect meaningful groups by clustering tombs using a distance to neighbour threshold of approximately half that size (2000 m). Figure 5B shows the result of using a dbscan method (Ester et al. 1996) to generate groups of sites using this threshold. The resulting 9 clustered groups of mounds make visual sense, with groups concentrated by proximity. Thinking about possible social interpretations of these groups an obvious further question, given the mark correlation function results, is whether the mound sizes found within each group exhibit a non-random hierarchical distribution, in which each group contains at least one of the larger tombs, followed by a number of medium size tombs and finally a number of small tombs. For the sake of clarity, in this paper 'hierarchy' is understood as ranking, in the sense that our argument is about the distribution of tomb sizes characterised by their rank size. To our knowledge there is no existing technique available to assess the hierarchical nature distribution of tomb sizes given the uneven number of tombs with known sizes per group, so we introduce here a novel test based on permuting the rank order sizes of the tombs (for a complementary approach see Hennig and Lucianu 2000). The approach (see also the Supplementary Data and Methods) operates as follows: first, rank all tombs in descending order of size, so that the largest is ranked 1 and the smallest is ranked n, where n is the number of tombs. Next, create a hierarchy of tomb-size ranks for each group, so that the biggest tomb in each group is placed in hierarchical level 1, the second biggest tomb in each group is placed in level 2, and so on. The result is a set of tomb-size ranks for each hierarchical level, which then allows us to compare the mean and/or sum of ranks at each hierarchical level with what we would get if we simply allocated tombs to size-levels with no reference to their group membership. For example, if the distribution of tombs across groups was perfectly hierarchical in the sense suggested above and there were, say, four groups of tombs, then the tomb-size ranks in level 1 should be 1, 2, 3 and 4, the ranks in level 2 should be 4, 5, 6 and 7, and so on, which is exactly what we would expect if we ignored group membership. Since, in reality, not all groups contain equal numbers of tombs, we adjust the expectation accordingly, i.e. the hierarchical levels will not have equal numbers of members. If we compare the observed ranks per level (i.e. those expected if a perfect size hierarchy pertained), we see ( Table 3: IdealMeans) that the ideal mean ranks increase from one level to the next whereas that is not always true of the observed mean ranks. We also see that the level 1 and 2 ideal mean ranks are substantially numerically smaller (i.e. higher ranking in descending order of size) than their observed equivalents. The conclusion that can be drawn from this is that the observed allocation of tomb sizes across groups is not perfectly hierarchical, but it may nevertheless be closer to such an idealised hierarchical distribution than we might expect by chance. We can examine this possibility by conducting a Monte Carlo significance test (Hope 1968) to establish how frequently the observed mean rank at a given level in the hierarchy is numerically lower (i.e. higher ranking) than the equivalent mean ranks obtained in a large number of simulations in which tomb sizes are allocated independently of group membership.
This entails repeatedly permuting (randomly shuffling the tomb sizes) and re-running the allocation of ranks to levels by group. We conduct 999 such simulations and then in accordance with standard practice (Manly 1991) derive p-values by comparing the extremeness of the observed mean ranks at a given hierarchical level to the simulated mean ranks at that level ( Table 3: Pval). Note that we have not attempted to compute a p-value for the distribution of ranks across all hierarchical levels as it is a moot point whether the latter would constitute the testing of multiple hypotheses and so require correction of the p-values to control the family-wise error rate (Holm 1979). The result of this Monte Carlo simulation suggests that the observed distribution of tomb sizes in each of the first two hierarchical levels is closer to that expected of a hierarchical distribution (each group having a tomb ranked highly by size) than would be expected by chance (i.e. p <= 0.05). In fact, this is true for three out of the first four hierarchical levels. Consequently, our test supports the argument that the largest tombs are distributed across the groups in a way that is broadly (albeit not perfectly) hierarchical, to an extent that is unlikely to occur by chance alone.
[Caption] Table 3. Results of the permutation tests on megalithic site sizes by cluster group.

Discussion
Modern spatial statistical methods allow us to move beyond a now rather stale debate about environmental determinism, because they facilitate empirical investigation of the interplay of different causes, as opposed to the a priori assertion of primacy according to theoretical preference. In this case, we have been able to demonstrate that distribution of megalithic mounds in the Monte Penide region reflects a preference for locations with particular environmental properties, but at a more local scale the spacing of these mounds seems to reflect some kind of social partitioning of the landscape.
The results described above allow us to conclude that megalithic sites in the Monte Penide region concentrate at specific elevations (300-500 mASL) close to the ridgelines that define the main watersheds draining to the sea. Once these major locational trends are accounted for, sites still exhibit clustering within 1km of each other, probably implying that once some megalithic landmarks were established, these encouraged the construction of new ones nearby. These patterns elicit interesting archaeological interpretations. First, they suggest that different groups in the regions gradually built up mound groups on the upland side of possible local territories (see also Bakker 1976Bakker , 1991Bradley 1991a;Criado Boado, Villoch Vázquez 2000;Murrieta-Flores 2012), on ridge tops draining to the sea. These locations also placed such mounds in highly visible locations on major routes of upland travel across the landscape (see also Carrero-Pazos 2018b), but we make no attempt here to tease apart which of these motivations, if any, might have been the most important one for the mound-builders.
Beyond the broader locational preferences, the distribution of mound volumes/sizes in the nine identified clusters exhibit a non-random hierarchical pattern, with a larger mound per group and then smaller ones around that, with what appears to be a preference for the large example to be at or near the meeting point of several watersheds (and upland ridge-routes). At risk of straying too far beyond the evidence, this observation suggests that the Monte Penide region was probably occupied by several groups of roughly equivalent status. Furthermore, without wishing to delineate strict territories without better dating and wider evidence, it is still a useful thought experiment to imagine rough zones of activity that often gave each of these identified groups a strip of land from mountain to sea of ca.70 sq.km, and larger still if such areas also extended over into the inland valleys.
Such conclusions fairly match with the spatial structure of megalithic landscapes in other Iberian and European regions, such as for example in Falbygden, Sweden (Sjögren 2010) where mounds tend to be found in marginal zones and in clusters, with a large tomb in the centre of the area. In that case, the interpretation was that such patterning related to social factors, such as the organisation of settlement or social Furthermore, differences in tomb orientation may in certain circumstances indicate specific cultural or funerary trends (Hoskin 2001;Silva 2010;Prendergast 2016), but this is still under initial investigation in Galicia (see e.g. Corujo-Tilve, Domínguez Márquez 2014;González García et al. 2017;2018).
Future work in Galicia would clearly benefit from better radiocarbon dating of mound use-lives and mound groups, and some efforts in this direction are already underway in the Barbanza peninsula, west of Galicia (Rodríguez Rellán, Fábregas Valcarce in press) in order to elucidate a better 'biography' of megalithic sequences (Rodríguez Casal 1989;Alonso Mathías, Bello Diéguez 1995;Mañana Borrazás 2005;Bóveda Fernández, Vilaseco Vázquez 2015). Likewise, there is much still to infer about what the possible territories really mean in terms of anticipated community sizes and social structure, something more visible in later periods such as the Bronze Age. However the above exploration already suggests new theory and method that might be applied to a much wider set of Galician and European megalithic evidence.

Acknowledgements
The archaeological dataset used in this work was kindly provided by Prof. Antón A.
Rodríguez Casal from the University of Santiago de Compostela, as a part of the project "Archaeology and Ecology of the Megalithic Complex in the South of Galicia" (1998Galicia" ( -2000, funded by the regional government of Galicia (Xunta de Galicia). The digital elevation model (LiDAR based) was obtained from the Spanish National Cartographic

Supplementary Data and Methods
To enable re-use of our materials and improve reproducibility and transparency, we include the raw data and R code used for all the analysis and visualisations contained in this paper in our supplemental online material at (http://dx.doi.org/10.17632/3sb4hwrrw9.1). The data used in this paper is released under the CC-BY4.0 license and our code with the MIT license to encourage maximum reuse.