Bronze Age communication and trade networks from Mediterranean Ivory records

doi:10.21203/rs.3.rs-3410688/v1

Download PDF

Research Article

Bronze Age communication and trade networks from Mediterranean Ivory records

https://doi.org/10.21203/rs.3.rs-3410688/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 20 Jul, 2024

Read the published version in Archaeological and Anthropological Sciences →

You are reading this latest preprint version

The transportation of Bronze Age ivory raw material and artefacts across the Mediterranean has been in the focus of archaeological research for over a century now. However, tracing the flow of ivory has mostly been restricted to traditional theoretical models of raw material distributions deriving from socio-culturally centred considerations. Environmental conditions, potential transportation networks and dissemination routes have not yet been considered decisive for the spread of ivory raw material from the African shores and the Eastern Mediterranean towards the Iberian Peninsula. Implementing computational environmental and archaeological modelling, we present a fully reproducible quantitative approach to estimate potential communication and transportation networks based on environmental covariates. We deploy a Network Analysis model and a predictive model based on Least Cost Path density to propose a potential land- and sea-based movement corridor for the western Mediterranean Basin that could have enabled the cultural spread of ivory during the 3rd and 2nd millennia BC.

Mediterranean Basin

Chalcolithic

Bronze Age

Least Cost Path analysis

Computational archaeology

Network analysis

Mediterranean ivory trade networks have captured the attention of researchers for well over a century. Louis Siret's pioneering work in 1913 marked the beginning of discussions regarding East-West connections (Siret 1913). Siret's hypothesis suggested an Eastern origin for certain aspects of Western societies. This continued to be the main narrative until the 1950s, with V. Gordon Childe being one of the most prominent proposers of the “ex oriente lux” model (Childe 1925). However, the advent of radiocarbon dating in archaeology during the 20th century revealed that Iberian megalithic monuments predated their Eastern counterparts (Renfrew 1973). In our specific case, one of the main consequences of the radiocarbon revolution was the dismissal of the idea of trans-Mediterranean relations, favouring an African origin for most Western ivory productions (Chapman 1990). However, the discovery of an ivory workshop in Valencina de la Concepción (Sevilla, Spain) between 2005 and 2006 (Nocete et al 2013) and findings by T. X. Schuhmacher (2011) in collaboration with the International Center of Ivory Studies (INCENTIVS) reopened the question of ivory supplying routes. Their work led to the identification of substantial amounts of ivory from Asian elephants (Elephas maximus), even in its raw form, on the Iberian Peninsula. Subsequently, Asian elephant ivory has also been scientifically traced in other parts of the Central-Western Mediterranean, particularly in Sardinia. These findings have strengthened the hypothesis of an ivory trade flow originating from the East (Morillo León et al 2018). However, the character of these contacts remains obscure, and there is ongoing debate between the traditional approaches of step-by-step exchange and long distance trade. With a focus on long-distance exchange and communication routes in prehistory, this paper presents the first model of ivory trade networks across the Mediterranean. Our model aims to contribute to a more comprehensive understanding of transregional contacts between 3000 BC and 1500 BC building on environmental preconditions. Our approach combines archaeological data with quantitative models based on movement affordances, which are usually described as the energy expenditure associated with travelling through raster cells where values represent cost (Conolly and Lake 2006). We employ computational modelling and Least Cost Path (LCP) analysis to understand supraregional dissemination patterns of ivory across the Bronze Age Mediterranean Basin. Furthermore, we use the recently developed network v.net.model algorithm by Ducke and Suchowska (2022) that is implemented in GRASS GIS 7 (Geographic Resources Analysis Support System, GRASS Development Team 2017) to simulate land-based and sea-based distribution graphs. Together with the site database, we present a fully replicable computational model (Marwick 2017) that adds environmental conditions into the network of material distribution to explain potential movement and transportation corridors as well as predict potential site locations across Northern Africa, the Iberian Peninsula, and Italy.

The archaeological literature has extensively explored mobility and Least Cost Path (LCP) analyses (see, for example, Whitley and Hicks 2003; Howey 2011; White and Barber 2012; Herzog 2020; Lewis 2021). While this has led to a strong theoretical and methodological foundation, it primarily focuses on land-based movement. Seafaring is often overlooked due to its added modelling complexity arising from fundamental differences in transportation modes, as extensively discussed in the literature. Drawing from this existing body of literature, and integrating water networks (Ducke and Suchowska 2022), our objective is to construct an environmentally grounded model. This model not only seeks to explain ivory distribution and dissemination patterns but also aims to contribute to the ongoing methodological discourse surrounding LCP analyses.

The premise of this work was the creation of an up-to-date database including all the available archaeological data on ivory finds from the Central-Western Mediterranean basin. The collected information includes the typology of finds, the type of ivory used to create the object, its location, a detailed description of the object, its current storage location, and the available literature. However, based on this database, simple spatial distribution maps of archaeological features are rather descriptive if not connected to explanatory variables from socio-cultural or environmental background. Hence, we integrate topographic and hydrologic variables to the analytical workflow in order to estimate the impact of terrain roughness and stream-flow character on potential trade network development. In this section, we first introduce into the site database, following a description of the geographical settings and variables used and a detailed methodical outline of the environmental data used as well as the model set-up.

2.1 Archaeological background and site database

The study of ivory production, among other exotic raw materials, in Chalcolithic and Early Bronze Age contexts has been a focal point since the beginning of Iberian Prehistoric studies (Siret 1913; Guilaine 2018). Chapman (1990) suggests an African origin for the majority of the Iberian ivories, with a relevant part of fossil tusks sourced from the sandy terraces of Iberian rivers. However, since 1999, an interdisciplinary team led by T. X. Schuhmacher, conducted a study on Iberian ivory with the surprising result of a significant amount of Asian elephant (Elephas maximus) ivory, a species that during this time frame was inhabiting some ecological niches in the Near East (Pfälzner 2016). This discovery reignited the long-standing debate concerning East-West relations during the third and second millennia BC, and the potential existence of extensive contacts and exchange networks, challenging the paradigm of a Mediterranean divided in two geo-cultural regions (Renfrew 1967).

Building on this, several researchers focused on ivory as a socioculturally relevant commodity used to display the status of Iberian elites during Chalcolithic and Bronze Age (Garcia Sanjuan et al 2013). Ultimately, the objective was to trace the supply routes and understand their significance in the context of cultural interaction processes across the Mediterranean. To achieve this goal, several attempts were made to identify the actual provenance of the ivories, using diverse approaches based on archaeometry or typology (Lafrenz 2004; Pascual Benito 2012; Morillo León et al 2018). By now, it is accepted that ivory played a central role in the social discourse of Prehistoric elites, becoming a strategic element to reinforce their social prominence in different societies and groups across the Mediterranean Basin, with higher or lower relevance depending on the area (Luciañez Triviño 2018; Garcia Sanjuán et al 2013).

Due to its intrinsic physical characteristics, ivory is especially suitable for carving, producing smooth, shiny surfaces that were appreciated by different social groups across the globe in different periods. Ivory furthermore encompasses affordances beyond mere material properties, partly due to its distant origin and cost and network-intensive acquisition (Van de Noort 2012; Needham 2009). The connection between transport of the material and its individual value is pivotal for comprehending the role played by exotica in traditional societies, and a very relevant factor in the cultural interaction developed at least from the third millennium BC onwards across the Mediterranean Basin, the Atlantic Façade and continental Europe. Particularly, in the Iberian Peninsula the appearance of ivory in archaeological context is linked to the rise of social complexity during the Chalcolithic (García Sanjuán et al 2013). Different types of items have been identified, following local typologies and even workshops where the raw material was processed, supporting the idea of long distance exchange routes (Braciela et al 2022; Nocete et al 2013; Pascual Benito 2012; Liesau and Schumacher 2012). Surprisingly, there are no clear precedents for this working technology in the bone industry, suggesting again rather a possible influence of external elements within the local society - either in the form of foreign individuals or local people imitating the techniques (Morillo León et al 2018).

To acknowledge these archaeologically and socio-culturally complex relationships, we developed a site database that covers 64 sites across the Mediterranean Basin, covering a chronological frame spanning from approximately 3000 to 800 BC, roughly in coincidence with the Chalcolithic and the Early Iron Age in most Mediterranean areas. A typological study has been intentionally avoided, focusing instead on the provenance of raw materials, as during this chronology the vast majority of items correspond to local types that are also known in other raw materials, as it is the case of the italian ossi a globuli (small bone plaques decorated with a row of bosses or globules), the V-perforated buttons or the Early Minoan seals (Morillo León et al 2018; Morillo León 2019). During the acquisition of the archaeological data, the Iberian Peninsula became a central focus region due to intensive research over the past years. On the other hand, and despite the abundance of ivory production in the Levant during the Bronze Age, a more traditional approach that doesn’t particularly focus on raw material supply led us to the choice of limiting our study area to the western Mediterranean Basin.

Published research was screened to compile a database containing information about ivory finds from sites within our spatial and temporal window. The literature review included scientific papers as well as several museum collections (MET, Metropolitan Museum, New York; MAN, Museo Arqueológico Nacional, Madrid, MAS, Archaeological Museum of Seville). The final database contains 411 ivory artefacts from 64 sites and geographically spans from Spain (n = 20 sites), Italy (n = 9), Portugal (n = 8) to Greece (n = 9) and North Africa (n = 24). The following site distribution characteristics were included in the data collection: id, coordinates (WGS84, EPSG:4326, lat/lon), faunal species, site location, and chronology (Table 1).

Table 1

Site location and samples used in this paper. Coordinates are in WGS84, EPSG:4326. Coordinates are derived from the literature and refer to the archaeological site. A list of references for each site mentioned is provided as Supplementary Information.
ID	Latitude	Longitude	Genus	Site name	Chronology
1	37.409	-6.057	Elephant	Valencina de la Concepcion-Montelirio	Chalcolithic
2	37.416	-6.077	N/A	Valencina de la Concepcion-Matarrubilla	Chalcolithic
3	37.409	-6.059	Elephant	Valencina de la Concepcion-Sector Urbanizacion Senorio de Guzman-Divina pastora	Chalcolithic
4	40.076	-3.855	N/A	Numacia de la Sagra-Las Mayores	Chalcolithic
5	38.238	-6.015	N/A	Huerta de Dios-Llerena	Chalcolithic
6	37.683	-1.700	N/A	Los Cipreses-Lorca	Early Bronze Age-El Argar
7	37.835	-3.904	N/A	Torre del Campo-Jaen	Late Chalcolitic-Early Bronze
8	37.770	-3.789	N/A	Marroquies Altos	Chalcolithic
9	37.723	-3.966	N/A	Las Penuelas	Chalcolithic
10	39.043	-3.498	N/A	Motilla del Azuer	Early Bronze Age-Las Motillas
11	37.137	-3.548	N/A	Cerro de la Encina-Monachil	Early Bronze Age-El Argar
12	37.321	-4.011	N/A	Montefrio	Late Chalcolithic-Early Bronze Age
13	37.727	-2.514	Elephant	Cerro de la Virgen	Late Chalcolithic-Early Bronze Age
14	37.803	-3.745	N/A	Puente Tablas	Late Bronze Age
15	39.024	-2.028	Elephant	El Acequion	Chalcolithic-Early Bronze Age
16	36.963	-2.522	N/A	Los Millares	Chalcolithic
17	37.297	-1.880	N/A	La Encantada	N/A
18	37.265	-1.788	N/A	Almizaraque	Chalcolithic
19	37.994	-1.144	Elephant	Murcia	N/A
20	37.558	-2.423	N/A	El Malagon	Late Chalcolithic
21	39.505	-8.615	Sperm whale	Galeria da Cisterna	Late Chalcolithic
22	38.728	-9.282	N/A	Leceia	Chalcolithic
23	37.157	-7.546	N/A	Nora	Late Chalcolithic
24	38.889	-9.074	Sperm whale	Verdelha dos Ruivos	Late Chalcolithic
25	38.585	-8.873	N/A	Palmela	Late Chalcolithic
26	38.881	-9.169	N/A	Dolmen de Casainhos	Late Chalcolithic
27	38.899	-9.435	N/A	Samarra cave	Chalcolithic
28	40.791	8.449	Sperm whale	Monte dAcodi	Late Neolithic
29	40.633	8.327	Elephant	Anghelu Ruju	Late Chalcolithic
30	37.680	12.756	N/A	Marcita	Late Chalcolithic
31	37.143	15.030	Elephant	Pantalica	Late Bronze Age
32	37.217	14.633	Elephant	Madonna del Piano	Late-Final Bronze Age
33	40.288	18.426	Hippopotamus	Roca Vecchia	Final Bronze Age
34	45.033	11.650	Elephant	Necropoli delle Narde	Final Bronze Age
35	45.157	11.354	Elephant	Lovara di Villa Bartolomea	Late-Final Bronze Age
36	39.561	8.902	Elephant	Padru Jossu	Chalcolithic
37	37.528	22.874	Hippopotamus	Asine	Middle-Late Bronze Age
38	37.731	22.756	Hippopotamus	Mycenae	Middle-Late Bronze Age
39	35.299	25.161	Hippopotamus	Knossos	Middle Bronze Age
40	35.063	24.795	Hippopotamus	Ayia Triada	Early Bronze Age
41	35.024	24.924	Hippopotamus	Platanos	Early Bronze Age
42	35.235	25.160	Hippopotamus	Archanes	Early Bronze Age
43	37.599	22.800	Elephant	Tyrins	Late Bronze Age
44	37.059	10.063	Hippopotamus	Utica	Early Iron Age
45	33.497	-7.104	Elephant	Khef el Baroud	Late Chalcolithic-Early Bronze Age
46	31.410	-9.839	Elephant	Sidi Harraz Cap Sim	Neolithic-Early Chalcolithic
47	31.496	-9.786	Elephant	Morgador	Iron Age
48	36.353	6.637	Elephant	Khanghet Si Mohammed Tahar	Neolithic
49	33.857	-7.061	Elephant	Rouazi Skhirat	Chalcolithic
50	36.038	6.570	N/A	Bou Zabaouine	Chalcolithic
51	32.738	0.936	N/A	Brezina	Chalcolithic
52	35.368	1.320	N/A	Columnata Tiaret	Chalcolithic
53	36.353	6.638	N/A	Damous el Ahmar	Chalcolithic
54	36.468	4.598	N/A	Djebel Gueldaman	Neolithic
55	36.191	5.409	N/A	Mezloug	Neolithic
56	34.165	3.497	N/A	Messaad	Neolithic
57	36.770	5.096	N/A	Pic des Singes	Neolithic
58	33.979	-6.900	N/A	Dar es Soltan	Chalcolithic
59	33.522	-7.826	Elephant	El Kiffen	Chalcolithic
60	31.176	-8.333	N/A	Grotte de Toulkine	Neolithic
61	35.764	-5.937	N/A	Grotte des Idoles	Neolithic
62	34.660	-3.418	N/A	Hassi Ouenzga	Neolithic
63	35.567	-5.367	N/A	Kahf Taht el Ghar	Chalcolithic
64	35.359	8.397	N/A	Abri Clariond	Neolithic

2.2 Environmental datasets

The study area stretches from Portugal and Spain to Greece and from Morocco to the Alps, covering the Central and Western Mediterranean Basin (Fig. 1). Such a broad and supraregional environmental coverage contains multiple climatic, topographic, and hydrologic systems with a strong variation in terrain roughness, precipitation and temperature gradients, as well as streamflow characters of the multiple drainage systems. River run-off behaviour depends strongly on the seasonal cycle and varies according to the snow-melt regime during early summer and winter precipitation patterns.

Fundamental basis to model movement corridors and landscape permeability estimates is a digital elevation model (DEM). Here, we employed the SRTM 90 v.4.1 tiles (https://cgiarcsi.community/data/srtm-90m-digital-elevation-database-v4-1/, last accessed on 8th of May 2023) obtained from the CGIAR Consortium for Spatial Information (CGIAR-CSI) website. Developed by NASA, this dataset offers comprehensive and up-to-date open data coverage on a nearly global scale. Compared to the original version provided by the USGS, the CGIAR-SRTM data have been processed to fill undesirable voids in the model (Jarvis et al 2008). The dataset can be downloaded in tiles measuring 5 by 5 degrees in size (lon-lat) and has a resolution of 3 arc of a second, which translates to 90 m nominal resolution at the equator. A further possibility for download is using the elevation function from the geodata package in R (Hijmans et al 2023).

While higher resolution data, such as the 25 m resolution EU-DEM (Copernicus: Earth Observation component of the European Union’s space programme), are available for Europe, we decided to use the SRTM data for several reasons. Firstly, the vast size of the study area, coupled with the exponential increase in computational requirements, made the SRTM data more practical and efficient. Additionally, sourcing comparable high-resolution data for North Africa would have been time-consuming, offering marginal gains for our specific model. Moreover, to address hardware limitations, we rescaled the SRTM raster to a factor of 2, resulting in an approximate 200 m resolution elevation map, which would have made the choice of a 25 m resolution DEM irrelevant. As already mentioned, the 90 m resolution refers to measurements at the equator. Because of the projection from geographical coordinates (EPSG:3857), it derives that the resolution decreases at increasing distance from the equator. For the study area, the average resolution is slightly less than 100 m, resulting in the aforementioned approximate 200 m resolution. The scripts, including the functions and the workflow adopted to mosaic, crop and re-scale the elevation model are included in the repository. To facilitate movement modelling across areas covered by water, these were filled with values = 0, making them "operable" for the LCP model. However, uniformly assigning 0 values would result in a slope value of 0 across all the Mediterranean. Thus, a friction map was required to assign specific movement costs for crossing different terrain types. A friction map is a raster map in which each cell represents the friction cost required to cross it, irrespective of its slope value (Novaline Jaga et al 1993; Howey 2007). These maps are employed to represent variations in the ease or difficulty of movement across different areas, such as land-cover types, rivers, or other features in the landscape. Nevertheless, assigning friction costs is often a subjective process, and their use can potentially introduce bias, leading to results that reflect the initial parameter choice. For this reason we decided to proceed with caution, keeping the friction value changes to a minimum. A uniform friction value of 1 was given to all land areas, while applying friction values of 0, 2, and 5 to sea areas at distances of less than 10 km, 10–50 km, and beyond 50 km from the coastline, respectively. Rivers were given a friction value of 1, with an additional buffer of higher friction (3) around them to represent ease of travel along the riverbank but increased difficulty in crossing the river. Hydrologic streamflow data derives from the Copernicus European Union’s Space programme (https://land.copernicus.eu/imagery-in-situ/eu-hydro/eu-hydro-river-network-database, last accessed 21st of May 2023). The data comes in spatial lines format, which were extracted and filtered according to their stream flow hierarchy (Strahler order number). All river sectors with Strahler order number > 7 were selected from the catchments of the rivers (in alphabetical order) Danube, Duero, Ebro, Garonne, Guadalquivir, Mesima, Pinios, Po, Rhone, Tajo and Tevere. In some cases, secondary or minor river basins are provided aggregated to larger ones (for example the Po river basin includes the Adige too) and therefore not listed here. For North Africa, rivers were downloaded using the rnaturalearth R package (Massicotte and South 2023), and in particular the ne_download function which allows to download rivers and lakes data at the desired scale (using the argument “scale = 10”, downloading the highest resolution available). The terra (Hijmans 2023), sf (Pebesma 2018; Pebesma and Bivand 2023) and dplyr (Wickham et al 2023) R packages are used to process the river data. The script used to download the North African data and to process them with the European one is provided in the repository (01_rivers.R script).

The friction map has a resolution corresponding to a rescaling of 4 times of the original SRTM (i.e. ca. 400 m). The reason why it was not kept to an aggregation of 2 is that it is largely repetitive and uniform and therefore it is not necessary to achieve a greater resolution.

The choice of our friction parameters is based on the assumption that sea movement was faster and easier than land movement and, if within sight from the coastline (< 10 km) preferred over land movement. Rivers were not assigned lower friction values compared to land, as large alluvial basins already exhibit flatter terrain, making them more efficient to travel on. It derives that introducing a different friction for rivers would have resulted in redundancy or even a duplication of this pattern. However, we acknowledge that crossing major rivers may have presented challenges or increased effort, hence the higher friction buffer around rivers. The size of this buffer is of ca. 800 m (2 raster cells) for the only reason of avoiding the 0-cost crossing at the edge of a rasterised line. To acknowledge the fact that this choice is arbitrary and may not be shared by everyone, we provided a fully reproducible workflow, that includes the scripts necessary to create the friction map, allowing readers to modify friction values based on different theoretical assumptions and compare their results with ours. The friction map is created using map algebra in R, using the terra (Hijmans 2023) and sf (Pebesma 2018; Pebesma and Bivand 2023) packages. The script used to build the friction map is the 02_friction400.R in the repository.

2.3 Model set-up

We distinguish two different potential movement patterns: a land-based and a sea-based circuit or route network. As of today, it is more than clear that a very important component of movement in Prehistory was water-based, as testified by many archaeological finds and studies (Broodbank 2013; 2016; Blake 2014; Iacono 2019). The literature is rich in papers about computational approaches to mobility and LCP analysis in archaeological research, which resulted in a very well developed theoretical and methodological background. However, this is generally limited to land-based movement and tends to exclude seafaring due to increased complexity in the resulting model (e.g. Howey 2011; Llobera et al 2011; Verhagen and Jeneson 2012; White and Barber 2012; Whitley and Hicks 2003). The reasons for this are well-known and mostly related to the difficulties of combining two very different methods of transportation-movement that has been often pointed out in the literature. However, there are some studies trying to combine or to consider in some ways both types of movements simultaneously (Ducke and Suchowska 2022; Leidwanger and Knappett 2018). The leitmotif is to keep the models very simple, based on traditional movement land-based costs based on isotropic movement, coupled with friction costs that include water areas (rivers and sea). The more complexity is added to a model, the more variables need to be integrated and the more issues arise from the certainly significant differences in land and water movement costs. The advantages of conceptually simple models are multiple. Universality (read: reusability) is more at hand and computational work becomes less intense and more accessible.

Least Cost Paths (LCPs) are models that calculate the optimal or most efficient route between two (or more) locations in the geographical space. They find application in many fields, from urban planning to ecological models and archaeology (Fabrício Machado and Miranda 2022; Fjellström et al 2022; Balbi et al 2021). These models need two source data in order to be computed: spatial location (coordinates) of the points from and to which we want to calculate the travel costs and a surface with resistance or cost characteristics on which the route is calculated. This surface can depend on several factors such as elevation, slope or land cover. The output of a LCP is a single line connecting the input points. The most evident caveat of this method in archaeology is that it relies on the assumption that the connected points were directly linked during the period under study. This can be straightforward in documented historical contexts, as demonstrated by several studies (e.g. van Lanen et al 2015; Lewis 2021). Unfortunately, this is not the case for prehistoric case studies. Moreover, even if we are confident that two points were in contact, the presence of additional sites may have influenced the choice of an optimal path. The fact that the archaeological record is by definition incomplete and that we may not know some of the sites belonging to this network make the reconstructions even less reliable. Site size, catchment and relative importance of a node in the network add further layers of uncertainty. Although these issues are widely discussed, sometimes assuming contemporaneity and completeness of the dataset (or at leasts that new inputs would not substantially change the picture) seems to be the only possibility in order to reconstruct past mobility networks (Ducke & Suchowska 2022).

As already pointed out (Herzog 2014, 2022; Lewis 2021) there are also a number of issues concerning the computational side of LCP analysis. In this paper, it is not possible to address a discussion over these topics but it is important to consider them and to elaborate our models accordingly. The most important issue concerns the choice of the cost models and of the LCP algorithm (Herzog 2022). Another relevant issue is the choice of an appropriate DEM. It has been demonstrated that LCPs models based on high resolution DEMs (e.g. 5 m) are highly affected by modern infrastructures such as roads or highways (Verhagen and Jeneson 2012). In addition, it has been proved that the performance of LCPs does not improve considerably using higher resolution DEMs (Herzog 2022). However, when smoothing an elevation model or using a coarser one it is important to remember that also natural features such as small creek valleys are smoothed as the slope values generally used to compute the travel costs.

Part of these problems involve computational and storage requirements. In fact, increased complexity in the model has a strong impact on the size of the output and time required (two times higher resolution results in a quadratic increase of these parameters; Herzog 2022). Nonetheless, the increase in accessibility to computational power allowed the exploration of many different alternatives, improving the methodological and theoretical background of cost-based analysis. One of these approaches is the so-called “From Everywhere to Everywhere” (FETE) approach (White and Barber 2012). Compared to traditional LCPs or many other workaround solutions, such as the “flow” of movement using hydrological modelling (Llobera et al 2011), this method does not require the existence of an origin or a destination. What it does is to sample a regular grid of points within the study area and to compute every possible connection between these points, resulting in a network of paths that indicate the presence of high traffic corridors. Here, we propose three additional implementations to the original model: (i) the inclusion of water-based movements using a friction map that takes water into account, following in part what was suggested by Ducke and Suchowska (2022), (ii) instead of a visual assessment of the structure and location of the corridors compared to the archaeological sites, we study them in a formal way using statistical analysis in order to find out whether sites are located on high traffic routes more than expected by a random distribution or whether sites are located closer to the corridors than a random distribution. In this case a distance surface from the densest corridors (calculated using a KDE) is used as covariate; (iii) we noticed that this approach is affected by edge effects, due to the fact that points along or closer to the edges of the study area are less connected (i.e., they miss the connection to further points extended out of the study area). This would result in a higher density in the central areas, regardless of the topographical characteristics of the landscape. For this reason, we created a sensibly large area (500 km in every direction) that will then be cropped out and not discussed in the results.

These improvements greatly increase the confidence of our conclusions and also increase the formality of the analyses, not only compared to the “traditional” FETE but also when compared to a more classic LCP approach. This in turn may help the LCPs to regain centrality in the archaeological debate. The main downturn of this approach is the high computational costs required, especially because the area we want to analyse is large and encompass ca. half of the Mediterranean Basin. Only considering the land mass within the study area our region is about 50 times larger than the Oaxaca region analysed by White and Barber (2012). Therefore, we opted to downsample the SRTM-DEM by aggregating the original cells with a factor of 4, resulting in an approximate resolution of 200 meters.

The model used in this paper has been implemented in R programming language (R Core Team 2022), calling functions from GRASS GIS via the package rgrass (GRASS Development Team 2022; Bivand 2023). The script used to calculate the LCP density works as follow:

a grid of points is sampled within the study area at regular space intervals within landmass. Ideally, the spacing between each point in the grid is set to 125 pixels (ca. 25 km). For computational reasons we decided to break up this grid of points in 500 smaller groups and run the analysis on each subset separately and summing it at the end.

For each iteration, the script calculates walking paths between a source point and all other points using the r.walk function. This function computes the anisotropic cumulative cost of moving from the origin point to the input set of destinations. The costs are based on the elevation and friction maps provided to the function. The function returns two outputs: the cumulative cost of moving between each raster unit (200 m grid) and the movement direction taken. This function is a great improvement compared to the r.cost function also available in GRASS because travel times are anisotropic rather than isotropic (walking downhill or uphill matters).

Once the previous step is completed, the r.path function reads the results and computes the least cost paths between the origin and the target points. The input maps are used in a way that the path stops when the direction is 0 or negative.

The resulting paths are then used to calculate the kernel density, calculated using the density function from the spatstat package, summed together and scaled to the number of layers and iterations.

Evidently, our model is merely based on environmental and topographical constraints in order to compute the optimal path. This is because we are not interested in the route taken by prehistoric people but to study the relations of the characteristics of the landscape in terms of mobility (a sort of affordance) and the distribution of archaeological sites. Statistical correlation or lack of it is important to get a better grasp on human-environmental relations and the possible impact of socio-cultural factors.

2.4 Network and prediction model of Ivory sites

Using the recently developed network model algorithm v.net.models by Ducke and Suchowska (2022), we deploy the GRASS GIS interface in QGIS. Input data is a DEM with 200 m resolution and a projection in EPSG:3857. NA values were replaced with the value 1, which essentially fills the sea surface to prevent a cost value of 0. All negative values were equally set to 1 to remove artefacts and sinks from the DEM. A cost surface was produced using the r.cost algorithm in GRASS GIS as outlined above. For the v.net.models function, the input variables were set to default, the starting points were selected to be the 64 sites and the input raster is the cost surface. This set-up produces a “complete connectivity model” where all sites are interconnected (for more details, see: https://github.com/benducke/v.net.models, last accessed 21st of June 2023). A second cost surface was produced using a composite raster of the land mass and an additionally generated sea-surface raster that represents heterogeneous altitude from random point distribution. For this reason, a sample of 500 '000 points was generated across the buffered spatial extent (100 km) of the previous network model. The sample was assigned a consecutive number of values, buffered with 10 km and rasterized to the same spatial extent and resolution as the initial DEM. The raster was divided by 100 to range from 0 to 6082 m with peak values outside the land mass. From the extent of the DEM and the study area rectangle, a difference mask was created and the produced raster was cropped to the difference polygon. Finally, both rasters were merged, which resulted in a continuous raster representing the original DEM as well as heterogeneous raster values in areas covered by sea. From this raster, a cost surface was produced that includes fictive terrain roughness in all areas covered by water and hence do not automatically favour water-based transport but weighs more heavily land-based transportation systems.

Using the prediction function from the spatstat package, we leverage the calculated least cost path density as a covariate to make predictions about high probability locations (Baddeley et al 2016). These predictions are based on the point process and its relationship with the covariate. In essence, one is using the information derived from the least cost path density as a factor to estimate where high probability points are likely to be within the context of the point process. This approach permits to identify areas that are more favourable or significant according to the relationship between the point process and the least cost path density covariate. This involves several steps. Firstly, the data must be converted to a metre-based projection to be incorporated into the spatstat workflow. Next, we convert our raster data into a format recognized by spatstat, a real-valued pixel image using the as.im function. Simultaneously, we format our point pattern into an object of class "ppp" which stands for a planar point pattern, using the ppp function (Baddeley et al 2016).

With these preparations in place, we can employ the rhohat function from the spatstat package (Badelley et al 2012). This function provides a nonparametric estimate of intensity based on the covariate. In practical terms, let's imagine we have a collection of points labelled as "y", which represents spatial locations. Alongside these points, we have associated values, denoted as "X(u)", representing the spatial covariate at each respective location "u". We treat these points "y" as realisations of a spatial point process, denoted as "Y", where the intensity function "λ(u)" depends on the covariate through the equation λ(u) = ρ(X(u)) (Badelley et al 2012).

Our goal is to determine the function "ρ", also known as the resource selection function in the context of ecological site preference modelling (Badelley et al 2012, 2016). To achieve this, we create a plot where we compare the ratio estimator (referred to as "sites'') against the actual sites, incorporating a 95% confidence interval. We make an assumption that our point process follows a Poisson process, a straightforward model for events occurring randomly in space or time (Grandell 1976, 1; Orton 1982) (see code to this article for additional and supplementary plots of the point process).

On this plot, we observe "tickmarks" along the horizontal axis. Each tickmark corresponds to the parameter value associated with observed sites. Clustering of these tickmarks indicates areas where there is a higher intensity of sites, suggesting a preference for specific covariate value ranges. Peaks in the curve highlight locations with a denser concentration of points compared to their immediate surroundings. However, the significance of this intensity is determined by the confidence interval. A wider confidence interval implies lower confidence and might indicate the presence of isolated values or outlier clustering behaviour. Conversely, when we observe high site intensity along with small confidence intervals, this indicates high significance and, in such cases, a strong preference for a particular range of covariate values (Kempf 2021; Kempf & Günther 2023). In summary, this approach involves preparing and analysing spatial data to understand how the point process relates to a covariate, making use of nonparametric estimation and statistical visualisation techniques.

The subsequent step involves using the obtained results and the covariate as input variables for a prediction model. This model aims to generate a probability map, which indicates areas where potential sites could be distributed based on the raster data and the empirically observed point process. The predict.ppm function available within the spatstat package is employed to create this probability map (Badeley & Turner 2000). This map is rescaled to a range from 0 to 1, where a value of 1 signifies a high probability of finding sites. In essence, it provides a spatial representation of the likelihood of site presence based on the relationship between the covariate (least cost path density) and the point process.

Using the point process model (ppm) to represent our site data in ppp format and the LCP density in im format obtained from a stars object (sites_ppp ~ im_stars), we conducted 1000 iterations to simulate site locations based on the provided input variables. Subsequently, we assigned geographical coordinates to the ppm using the EPSG:3857 projection and calculated the intersection between the simulation results and the landmass window. To optimise computational efficiency and smooth the vector, we simplified the mask's outline.

This intersection process generated 33,458 points located on land, from which we performed a Kernel Density Estimation (KDE) with a sigma of 40 km. After trying different sigma in the range of 10–50 km we opted for 40 km due to a better representation of the local vs the regional density. In order to visualise areas with a high probability of containing sites, we established a statistically derived threshold for the KDE values. This threshold was defined as the mean plus the standard deviation (mean + SD), and all KDE values below this threshold were designated as NA. Conversely, all values surpassing the threshold were set to 1, creating a high-probability vector representing the potential distribution of sites. This vector was then superimposed onto the rescaled LCP density corridors (see Fig. 2). The scripts used to compute the rhohats and the predictive models are available in the “05_analyses.R” script, alongside with the code necessary to generate plots to visualise intermediate steps not shown in the paper.

In this section, the results from the two models are presented and discussed in light of Bronze Age ivory distribution and transportation across the Mediterranean Basin. Where plots and figures of the analyses are not provided in the actual paper, we refer to the repository and the code underlying this article.

3.1 Prehistoric ivory network model

The v.net.models results allow to estimate two different potential networks for transportation and raw material dissemination across the Mediterranean Basin (Fig. 2). From the first model, represented as blue lines in Fig. 2, a combined sea/land-based communication system can be derived that strongly depends on coastal transport and integrates sea-faring as the central part of connecting sites in the network. In essence, sea-based movement is enabled by assigning a very low travel cost to areas covered by water. Certainly, this assumption can only be used in the model environment and does not meet realistic standards of environmental variability. On the other side, it allows us to estimate the linkages of sites located close to the shoreline and how these empirically observed sites are connected to other sites within the network. Particular transportation hotspots can be identified that connect, for example, Sardinia with the north African coast or the Iberian Peninsula. Other important hubs in the network are, as expected, Sicily that links the east-west transportation system. The strong north-south network that is assumed by the model to connect north-western Africa and the Iberian Peninsula across the Strait of Gibraltar and the Alboran Sea to the east, however, seems to be less likely. Due to the predominant currents that characterise this region, sea-based transport during prehistory and antiquity is considered to be problematic. More likely, there have been particular routes that facilitated the crossing of the sea and that have been used more frequently depending on weather conditions and seasons. The second model output describes a combined sea/land-based network in which travelling across water is assigned to high energy expenditure. Eventually, this model outputs a network that strongly depends on land-based transportation systems (black lines in Fig. 2). Remote sites located more towards the inland show equal connection to neighbouring sites due to the same terrestrial values of the coast surface. Meeting the coastline, however, the network differs considerably from the rather sea-based transportation system. As expected, the connection of important nodes, such as Sardinia and Sicily are now connected via the shortest passages to travel across the Mediterranean Sea. In this model, the Strait of Gibraltar is the preferred axis that connects north-west Africa and the Iberian Peninsula. Still, coastal transport prevails in most areas, due to the topography. However, now mirrors rather inland transportation compared to the potential near-coast sea-based network of the first model output.

3.2 Least Cost Path density

The results of the LCP density can be also useful before further analyses to get a better understanding of the global landscape affordance in terms of travel and to assess the suitability of the proposed model. It is evident that part of the results are strongly influenced by the decision we made on the parameters of the model. The most evident result is the clear preference for travelling by sea at short distances from the coastline, which is the area with the lowest friction cost. Secondly, the way points were sampled cause peripheral regions to have a lower LCP density. This is due to the fact that the points located in these areas, even though they are destinations and starting points, are crossed only rarely (they have less neighbours). Edge effects have been accounted for cropping the study area to a smaller extent but there is still the possibility of underrepresenting the border zones. Moreover, the sampling method tends to underrepresent some islands that are not situated along major movement corridors. For example, the Balearic Islands have fewer simulated points within their boundaries, leading to a reduced number of simulated sources and destinations for the model. This is, in fact, a direct consequence of the islands' unique geographical and morphological attributes. It is important to note that this phenomenon is not a concern within our model; rather, it may contribute to detect possible deviations from environmental constraints in the archaeological records.

Additionally, the model attributes only a marginal role to rivers as corridors. It is essential to note that further reducing the friction of rivers, such as setting it to zero as in coastal areas, would lead to potentially redundant predictions. Rivers are naturally conducive to travel due to their modification of the landscape and the creation of alluvial plains and river valleys, which are already easier to cross compared to the surrounding terrain, due to the relatively low slopes. Consequently, making rivers even more suitable for travel would result in a self-sustaining prediction, potentially reducing the informativeness of the model, as already discussed.

A last remark due to the model is the lack of clear preferences for specific corridors in North Africa. This is due to the fact that South and East of the Atlas mountain range the morphology of the terrain is rather similar. In this case it is likely that the most important parameters were ivory availability (whether it was elephant or hippopotamus) and the presence of easy crossing areas in arid or semi-arid areas (water availability, oasis, or others). However, reconstructing the exact location of these is a challenge that would certainly require more time and field campaigns that could not be carried out at this stage.

In addition to these more technical remarks, there are several interesting findings from a visual assessment of the LCP density. Notably, the central-eastern Mediterranean reveals a significant European-African connection in western Crete, which is the closest point on the island to both Africa (Cyrenaica, Libya) and mainland Greece (SE Peloponnesus). The proximity of intermediate islands such as Kythera and Antikythera to the NW and Gavdos to the S, facilitated sea crossings. Contacts and exchange between Minoan Crete and Egypt are far known, including the depiction of Minoans carrying elephant’s ivory tusks in the Tomb of Rekhmire (Krzyszkowska 1990).

Similarly, the Otranto channel is the narrowest part of the Adriatic and the model shows it as a possible movement corridor. However, our model suggests two potential alternative routes, one northern, cutting southern Apulia by land and one southern, following the sea route into the Ionian sea. Northern Sicily and northern Calabria also stand out as a high-density corridor, connecting the Ionian and the Tyrrhenian seas between modern-day Paola and the Gulf of Corigliano and between the Gulf of Policastro and the Ionian sea. Furthermore, our model identifies a high-density movement corridor between these two seas between Policastro and the Gulf of Manfredonia, cutting through the southern Apennines. The sea is then crossed between the Gargano Peninsula, the Tremiti islands and the Dalmatian coast. A final significant connection between the Adriatic and Tyrrhenian is through a land corridor that cuts the Friulian plain and crosses the Apennines, roughly following the Secchia valley and reaching the sea at the Gulf of La Spezia following the Magra river.

The most suitable connecting point between Africa and Europe in the central Mediterranean is Sicily. From its western shores it is possible to reach Cap Bon (Tunisia) from two alternative routes, one to the north of the cape and the other to the south. The second shortest path is then the corridor connecting Sardinia (from Capo Spartivento) to North Africa, approaching from The Galite archipelago. This corridor could represent the easiest option to travel between North Africa and the French Midi and northern Italy, with Elba island and Corsica being part of this larger high traffic area.

In general, the areas exhibiting the highest LCP density include SE Spain ( Cartagena-Alicante-Valencia), NE Spain and S France (L’Escala-Perpignan-Gulf of Fos-Toulon), NW Italy (Liguria), the region between north-central Sicily (around Cefalù) and northwestern Calabria (modern-day Paola), and areas around Algiers and Annaba in Algeria, Bizerte, and the region between Hammamet and Monastir in Tunisia. Notably, many well-known classical and modern cities and archaeological sites are found across all these regions, suggesting a strong connection between their location and travel routes. Interestingly, the Gibraltar Strait, while serving as a node in the travel network and LCP corridors, does not rank among the highest-density areas, possibly due to the relatively lower number of points simulated on the Atlantic coast compared to the rest of the study area. Additionally, three other crossing points between southern Iberia and Africa are identified: Cartagena (Spain) and the area north of Mostaganem (Algeria); Cabo de Gata (SE of Almeria, Spain) and Oran (Algeria); Punta Entinas (SW of Almeria, Spain) and Cap des Trois Fourches (north of the Spanish enclave of Melilla), via the island of Alboran (Isla de Alboran, Spain).

3.3 Predictive modelling results

The results of the prediction of the point process and the covariate are shown in Fig. 3. Compared to the LCP density, this model is not meant to identify movement corridors but identify the best areas for sites using the LCP density as the predictive covariate. For this reason, water areas have been excluded from the model. It is worth noting that while the LCP density model relies solely on environmental variables, the predictive model considers the distribution of Bronze Age sites with ivory as input points, incorporating the information about point process.

The model shows a very high predictive power for the coastal areas of the western Mediterranean basin, from Ceuta to Melilla with some interruptions, maintaining a seamlessly high probability along all the North African coast up to Hammamet (Tunisia) and then extending to the rest of the Tunisian coast up to the island of Djerba, with minor interruptions. On the northern coastline of the Mediterranean, we observe a similar pattern of nearly uninterrupted predicted area from Gibraltar to the island of Elba (Italy). Additionally, south and eastern Sardinia and Corsica, exhibit high predictivity, reaffirming the previously mentioned patterns of movement between Tyrrhenian Italy, the French Midi, and northern Africa. Finally, the areas between Gibraltar and Cadiz (Spain) and the eastern fringes of the Po plain emerge as prominent areas in the predictive model. There are, however, some areas with patchier predicted density, notably along the Moroccan Atlantic coast, southwestern Portugal, southwestern Campania, northwestern Calabria, and the southernmost part of Apulia (Italy).

In Southeastern Europe, predictability is generally low, likely due to the scarcity of ivories in the region compared to other areas. Nevertheless, the model predicts some smaller coastal areas in this region, including Corfu, western Crete (Greece), and a few Dalmatian islands. In contrast to the LCP density, the prediction model identifies more restricted areas of higher predictability in the region south of the Atlas mountains. There are primarily two clusters of high predictability: one in central-northern Algeria (modern El Bayadh Province) and another in eastern Algeria and southwestern Tunisia.

3.4 Archaeological contextualization of the model

The predictive model and the observed LCP density presented in this study, exhibit consistency with both archaeological and historical evidence. Notably, the Iberian Peninsula is overrepresented due to extensive research in recent decades, in contrast to regions such as Italy, Greece, Southern France, and the Balkans, where ivory materials from the Bronze Age exist but have only been studied stylistically and typologically if at all (especially in the Balkans).

A first assessment of the model shows broader overlapping with the suggested directions of the cultural expansion of the Bell Beaker phenomenon across the Mediterranean Basin (Morillo León 2021). Moreover, the regions with increased connectivity shown in the LCP density, aligning with the most dynamic areas during the Chalcolithic and Early Bronze Age, underscoring the significance of mobility in shaping the Bronze Age as a socio-cultural phenomenon (Kristiansen and Larsson 2005). More generally, the movement corridors appear to reflect the presence of hypogeic burials, as described by Guilaine (2015).

A closer look at the LCP density shows further similarities between the easiest paths and the archaeological literature. To the west, the estuary of the Tagus River was a core region of the Megalithic phenomenon and a potential origin of the Bell Beakers horizon (Müller and Van Willigen 2001). Archaeologically, it exhibits close connections with the Guadiana mouth and the Guadalquivir river, along which Valencina de la Concepción (Seville, Spain) is located. These connections are established through river valleys and the Atlantic coasts, forming a Southwestern Iberian region. While this Southwestern focus has long been recognized by scholars based on archaeological evidence, the model's contribution lies in outlining the path along the Atlantic shores, both Iberian and Moroccan. This axis is supported not only by the presence of the same raw materials (L. africana) but also by a range of typologies supporting contact with Southwestern Iberia, including palmela arrowheads and Bell Beaker pottery f the same Maritime style (Schuhmacher 2014) The existence of local megalithism, though scarcely studied, suggests an earlier phase of contact during the Early Chalcolithic, potentially contributing to the origins of the Maritime/International Bell Beakers decorative style. Based on the fact that the comb-impressed decoration, while absent in the Iberian Neolithic, is well known in the Maghrebian Neolithic. This interaction axis likely played a role in the emergence of the Bell Beaker Horizon in the Atlantic and Mediterranean regions (Morillo León et al 2018). However, the predictive model shows this connection only partially, due to the lower density of LCP caused by the aforementioned reasons.

Southeast Iberia is another dynamic area during this period, marked by cultures such as Millares (Chalcolithic), El Argar (Early Bronze Age), Bronce Valenciano, or Culture of las Motillas, representing the transition from Bell Beaker Chalcolithic to the Bronze Age. Again, a network of sea routes connects this region to the North African shores (Fig. 2). Interestingly, during the Chalcolithic, most of the elephant ivory present in this area is of Asian elephant origin (E. maximus), with African elephant ivory appearing only in the Early Bronze Age (Schuhmacher 2011), along with hippopotamus ivory (Pau et al 2018). This evolution in ivory supply aligns with findings in the Sardinian hypogeum of Padru Jossu (Sanluri, Sardinia) (Morillo León et al 2018) and the well-documented historical connections between the Iberian Peninsula and Sardinia in Late Prehistory. However, the results of our model suggest that this connection may not have occurred directly or through the Balearic Islands, which is in line with the known challenging sailing conditions around these islands until at least the Bronze Age (Guerrero Ayuso et al 2007) (Figg. 2–3). Instead, the corridor passes through the Mediterranean Maghreb. Hotspots in Eastern Algeria and Cap Bon in Tunisia are noted, where regional megalithism has developed but remains largely unexplored (Camps 1995).

Cap Bon is the closest connection to Western Sicily, likely through Pantelleria Island, where an ivory bracelet of unknown origin has been documented (Ardesia et al 2006). Sicily acts as a crucial link between the two Mediterranean spheres as initially proposed by Renfrew (1967). The island exhibits connections with other regions, particularly related to the Lipari and Pantelleria obsidian trade network since the Neolithic (Nicoletti 2001). However, during the Bell Beaker period, Western Sicily becomes fully integrated into the Western Mediterranean network, as demonstrated by Guilaine (2009). Notably, a local variety of painted beakers emerges, representing a typological hybridization of Western impressed pottery with Aegean painted productions, belonging to the oriental tradition (Guilaine 2004).

Moving to the Adriatic region, the model reveals lower levels of connectivity, except for the traditional route from southeastern Apulia to Corfu (Fig. 3). In contrast, the v.net model shows a more intense interaction between southeastern Sicily, southern Calabria with the Peloponessos and southern Crete, where contacts are archaeologically observed and the tholos-type burials are prevalent (Heyd 2007; Morillo León 2019) (Fig. 2). These burials have also been connected to other Mediterranean regions, as well as other typologies in material culture, making them a possible indicator of long-distance contacts (Morillo León 2019).

We have chosen to exclude the chronology associated with the international ivory style in this region due to the extensive volume of material involved, which exceeds the scope of this study. It is reasonable to assume that this axis likely extends to include Cyprus and the Levant, which were suppliers of Asian elephant (E. maximus) ivory. However, the routes to Anatolia remain unclear, despite notable typological similarities with the western Mediterranean (Gonzalez Prats et al 1994; Heyd 2007; Morillo León 2021; Schuhmacher 2013, 2014). However, these areas are outside the study area of this study and are not further discussed.

Regarding the supply of Asian elephant ivory, it is noteworthy that it is absent throughout the entire study chronology in the Aegean and adjacent areas, where hippopotamus ivory serves as the primary raw material. The lack of characterization analyses in Sicily prevents the assignment of a specific role to this region, but its significance, especially in relation to the Maghrebian corridor, should not be overlooked. This underscores the importance of further on-field research in the underexplored North African region, which must be integrated into the broader Mediterranean cultural context.

Lastly, it is worth highlighting some similarities between the predictive and the ivory network models to the classic literature describing the African routes. In this regard, the two better-known reports are Hanno's Periplus (Hair 1987) and Balbo's expedition (Santana and Arcos 2009). In the former, the Phoenician navigator Hanno sailed to the Atlantic shores of Africa, possibly as far as Cap Sim, where abundant elephant herds were found, coinciding with the ivory supply to Western Iberia, as discussed previously. The latter expedition, led by Lucius Cornelius Balbus the Minor, traversed the Sahara Desert up to the Niger River, encountering hippopotamus in the area. Notably, Balbo's journey began in the Western Libyan shores and then inland, an area with high probability in our predictive model (Fig. 3). Certainly, these reports alone cannot be taken as a validation for our models, but they at least seem to support the observed patterns. Moreover, Balbo's Phoenician heritage from the Roman Gades (Cádiz) could reveal the persistence of traditional knowledge that could date back as far as to the Bronze Age, when hippopotamus and African elephant ivory replaced Asian ivory in Southeast Iberia and Sardinia. Following this argument, the integration of African communities during the Bronze Age seems evident, as testified by the depiction of Iberic types of weapons (hallbards, daggers and arrowheads) in North African rock art (Souville 1986) or the presence in the area of these types, either imported or locally produced (Chenorkian 2000).

In conclusion, this study presents a model of Bronze Age ivory distribution, travel affordance and probability across the Western Mediterranean. We combined two different models both based on environmental variables and site location, allowing for a more complete assessment of site distribution and testing the recently published v.net. models algorithm implemented in GRASS GIS. The network model identifies potential routes based on travel costs between existing sites, possibly indicating pathways for raw material and artefact dissemination. The results highlight the presence of significant hubs, such as Sardinia, Sicily, and the Strait of Gibraltar. However, the lack of sites in some areas, if caused by a research bias, would result in their underrepresentation. In order to overcome this issue, we employed a predictive model based on a FETE LCP density and the ivory point process. This is due to the fact that the covariate used relies solely on environmental variables, such as elevation, distance from the coastline, and the river network. This model highlights potentially favourable areas for transportation hubs within the study region.

A comparison between the proposed model, the distribution of ivories, and the potential source locations of raw materials (mostly Asian and African elephants and hippopotamuses) demonstrates consistency with historical and archaeological evidence, aligning with the state of the art in Late Prehistoric Mediterranean studies. Key findings correspond with the cultural expansions associated with the Bell Beaker phenomenon and dynamic regions during the Chalcolithic and Early Bronze Age, highlighting the role of mobility in shaping socio-cultural developments. The models shed light on connectivity between regions, such as the Southwestern Iberian-Moroccan axis, interactions between Southeast Iberia and North African shores, and connections between southeastern Sicily, southern Calabria, the Peloponnesos, and southern Crete. The models also identify potential areas for future investigation, particularly in North Africa, where the research has not been intense but is assumed to hold high significance in the broader Mediterranean socio-cultural context.

Both models incorporate water-based movement, often overlooked in mobility studies. Furthermore, as we provide the coordinates of the sites, the scripts to produce the models and the environmental variables are free and open source, it is possible to adjust and test different initial settings for the model. The size of the sample (399 ivory artefacts found in 64 sites) and the uneven research intensity may constitute the main limitation of our model. However, its structure makes it more robust to changes and easily adaptable and implemented. Moreover, the model predicted areas without known ivory artefacts, which may prove accurate in the future, particularly in North Africa, where the likelihood of an incomplete dataset is higher.

Due to the need to delimit the study area and computational constraints, critical regions in ivory exchange networks, such as the Levant (a source of Asian elephant ivory), Cyprus (an important hub in the Eastern Mediterranean), and Anatolia, were excluded and should be included in future studies. This may be the cause of the lower LCP density and predictivity of the eastern part of the study area (Greece and the Adriatic shores of the Balkans). Similarly, this applies to the northern part of the study area, which had active connections with northwestern and central Europe (Kristiansen 2022), forming a trans-Mediterranean and European network stretching from the British Isles and Scandinavia to the Levant, the Near East, and West Africa.

The authors do not declare any competing interest.

Acknowledgements: GB’s, MK’s, and JMML’s work was supported by the Germany Research Foundation and the Kiel University (CRC 1266) under grant number 290391021. MK received a post-doctoral fellowship grant from the Swiss National Science Foundation (SNSF) at the University of Basel (CH) under the project number TMPFP2_217358, project EXOCHAINS - Exploring Holocene Climate Change and Human Innovations across Eurasia. This paper is supported by the open access DEAL agreement.

Code availability: During review the code can be accessed here: https://gitlab.com/bilottigiacomo/marfil

Author contributions: Concept: JMML, GB, MK; data collection: GB, MK, JMML; analyses and code: GB; MK; visualization: GB, MK; writing original draft: GB, MK, JMML

Ardesia V, Cattani M, Marazzi M, Nicoletti F, Secondo M (2006) Gli scavi nell'abitato dell'età del Bronzo di Mursia, Pantelleria (TP): relazione preliminare delle campagne 2001-2005. Rivista di Scienze Preistoriche LVI. Firenze: Istituto italiano di preistoria e protostoria. Permalink: http://digital.casalini.it/10.1400/206423
Balbi M, Croci S, Petit EJ et al (2021) Least-cost path analysis for urban greenways planning: A test with moths and birds across two habitats and two cities. J Appl Ecol 58:632–643. https://doi.org/10.1111/1365-2664.13800
Baddeley A, Chang YM, Song Y, Turner R (2012) Nonparametric estimation of the dependence of a spatial point process on spatial covariates. Statist Interface 5:221-236
Baddeley A, Rubak E, Turner R (2016) Spatial point patterns: Methodology and applications with R. CRC Press Taylor & Francis Group, Boca Raton
Baddeley A, Turner R (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42:283–322
Barciela González V, García Atiénzar G, López Padilla JA, Hernández Pérez MS (2022) Producción y consumo de marfil en la Edad del Bronce: Cabezo Redondo (Villena, Alicante) como marco. Cuadernos de Prehistoria y Arqueología de la Universidad Autónoma de Madrid 48(2):79–106. https://doi.org/10.15366/cupauam2022.48.2.003
Bivand R (2023) rgrass: Interface Between 'GRASS' Geographical Information System and 'R'. R package version 0.3-8, https://CRAN.R-project.org/package=rgrass
Blake E (2014) Social Networks and Regional Identity in Bronze Age Italy. Cambridge: Cambridge University Press
Broodbank C (2013) The Making of the Middle Sea: A History of the Mediterranean from the Beginning to the Emergence of the Classical World. Oxford: Oxford University Press
Broodbank C (2016) The Transmitting Sea: A Mediterranean Perspective. In Kiriatzi E, Knappett C (Eds.), Human Mobility and Technological Transfer in the Prehistoric Mediterranean. British School at Athens Studies in Greek Antiquity, pp 18-30. Cambridge: Cambridge University Press. doi:10.1017/9781316536063.003
Camps G (1995) Les nécropoles mégalithiques de l’Afrique du Nord, in: L’Afrique du Nord antique et médiévale: Monuments funéraires, institutions autochtones, Trousset P (Dir.), Paris, Comité des Travaux historiques et scientifiques, p 17-31 (6ème Colloque international sur l’histoire et l’archéologie de l’Afrique du Nord, Pau 1993)
Chapman R (1990) Emerging complexity. The later prehistory of south-east Spain, Iberia and the West Mediterranean. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511735486
Chenorkian R (2000) « Hallebarde ». Encyclopédie berbère. DOI: https://doi.org/10.4000/encyclopedieberbere.1644
Ducke B, Suchowska P (2022) Exploratory Network Reconstruction with Sparse Archaeological Data and XTENT. J Archaeol Method Theory 29:508–539. https://doi.org/10.1007/s10816-021-09529-3
Fabrício Machado A, Miranda FR (2022) The potential distribution of Cyclopes didactylus, a silky anteater, reveals a likely unknown population and urgent need for forest conservation in Northeast Brazil. Journal of Tropical Ecology 38(6):454–461. Available at: https://doi.org/10.1017/S0266467422000372
Fjellström M, Seitsonen O, Wallén H (2022) Mobility in Early Reindeer Herding, in: Salmi AK (Ed.), Domestication in Action, Arctic Encounters. Springer International Publishing, Cham, pp 187–212. https://doi.org/10.1007/978-3-030-98643-8_7
García Sanjuán L, Luciañez Triviño M, Schuhmacher T, Wheatley D, Banerjee A (2013) Ivory craftsmanship, trade and social significance in the southern Iberian Copper Age: the evidence from the PP4-Montelirio Sector of Valencina de la Concepción (Seville, Spain). European Journal of Archaeology 16(4):610–635. http://doi.org/10.1179/1461957113Y.0000000037
González Prats A, Ruiz Segura E, Gil Fuensanta J, Seva Román R (1994) Cerámica roja monócroma anatólica en el poblado calcolítico de Les Moreres (Crevillente, Alicante, España). LVCENTVM XI-XIII (1992-94):7-38. http://dx.doi.org/10.14198/LVCENTVM1992-1994.11-13.01
GRASS Development Team (2017) Geographic Resources Analysis Support System (GRASS) Software, Version 7.2. Open Source Geospatial Foundation. Electronic document: https://grass.osgeo.org
GRASS Development Team (2022) Geographic Resources Analysis Support System (GRASS) Software, Version 8.2. Open Source Geospatial Foundation. https://grass.osgeo.org
Grandell J (1976) Doubly stochastic Poisson processes. Springer, Berlin
Guerrero Ayuso VM, Calvo Trias M, García Rosselló J, Gornés Hachero S (2007) Prehistoria de las Islas Baleares / Prehistory of the Balearic Islands: Registro Arqueológico y Evolución Social Antes de la Edad del Hierro / Archaeological Record and Social Evolution before the Iron Age. BAR international series 1690. https://doi.org/10.30861/9781407301280
Guilaine J (2004) Les Campaniformes et la Méditerranée. Bulletin de la Société préhistorique française 101(2):239–249. https://doi.org/10.3406/bspf.2004.12991
Guilaine J (2009) La Sicile et l’Europe campaniforme, in La Sicile et l’Europe campaniforme: la collection Veneroso à Sciacca, Guilaine J, Tusa S, Veneroso P (Dir.), Toulouse, Archives d’Ecologie préhistorique, p 135-195
Guilaine J (2015) Les hypogées protohistoriques de la Méditerranée. Arles: Errance. ISBN: 978-2-87772-544-6
Guilaine J (2018) Siret’s smile. Antiquity 92(365):1247-1259. http://doi.org/10.15184/aqy.2018.112
Hair PEH (1987) The “Periplus of Hanno” in the History and Historiography of Black Africa. History in Africa 14:43–66. https://doi.org/10.2307/3171832
Herzog I (2014) Least-cost Paths – Some Methodological Issues. Internet Archaeology 36. https://doi.org/10.11141/ia.36.5
Herzog I (2020) Spatial analysis based on cost functions. In: Gillings M, Hacıgüzeller P, Lock G (Eds.), Archaeological Spatial Analysis: A Methodological Guide. Routledge, Taylor & Francis Group, London and New York, pp 333-358
Herzog I (2022) Issues in Replication and Stability of Least-cost Path Calculations. Studies in Digital Heritage 5(2):131–155. Available at: https://doi.org/10.14434/sdh.v5i2.33796
Heyd V (2007) When the West meets the East: The Eastern Periphery of the Bell Beaker Phenomenon and its Relation with the Aegean Early Bronze Age. In: Between the Aegean and Baltic Seas. Ed. by I. Galanaki et al Aegaeum 27 (Liège), p 91-107
Hijmans RJ (2023) terra: Spatial Data Analysis. R package version 1.7-39, https://CRAN.R-project.org/package=terra
Hijmans RJ, Barbosa M, Ghosh A, Mandel A (2023) geodata: Download Geographic Data. R package version 0.5-8, https://CRAN.R-project.org/package=geodata
Iacono F (2019) The archaeology of Late Bronze Age interaction and mobility at the gates of Europe: people, things and networks around the southern Adriatic Sea. New York: Bloomsbury
Jarvis A, Reuter HI, Nelson A, Guevara E (2008) Hole-filled SRTM for the globe Version 4, available from the CGIAR-CSI SRTM 90m Database. http://srtm.csi.cgiar.org
Kempf M (2021) Take a seed! Revealing Neolithic landscape and agricultural development in the Carpathian Basin through multivariate statistics and environmental modelling. PLoS ONE 16(10):24. https://doi.org/10.1371/journal.pone.0258206
Kempf M, Guenther, G (2023) Point pattern and spatial analyses using archaeological and environmental data – A case study from the Neolithic Carpathian Basin. Journal of Archaeological Science: Reports, 47. https://doi.org/10.1016/j.jasrep.2022.103747
Kristiansen K (2022) Bronze Age Travellers, in Fernández-Götz M, Nimura C, Stockhammer PW, Cartwright R (eds), Rethinking Migrations in Late Prehistoric Eurasia. London: British Academy. https://doi.org/10.5871/bacad/9780197267356.003.0005
Kristiansen K, Larsson TB (2005) The Rise of Bronze Age Society: Travels, Transmission and Transformations. Cambridge: University Press. ISBN: 9780521585675
Krzyszkowska O (1990) Ivory and Related Materials (Classical Handbook 3, Bulletin Supplement 59). London: British Institute of Classical Studies
Lafrenz KA (2004) Tracing the source of the elephant and hippopotamus ivory from the 14th century B.C. Uluburun shipwreck: The archaeological, historical, and isotopic evidence. Dissertation, University of South Florida, Tampa. https://digitalcommons.usf.edu/etd/1122
Leidwanger J., Knappett C (2018) Maritime Networks, Connectivity, and Mobility in the Ancient Mediterranean. In J. Leidwanger, C. Knappett (eds), Maritime Networks in the Ancient Mediterranean World (pp. 1-21). Cambridge: Cambridge University Press. doi:10.1017/9781108555685.002
Lewis J (2021) Probabilistic Modelling for Incorporating Uncertainty in Least Cost Path Results: a Postdictive Roman Road Case Study. J Archaeol Method Theory 28, 911–924. https://doi.org/10.1007/s10816-021-09522-w
Liesau von Lettow-Vorbeck C, Schuhmacher T (2012) Un taller de marfil del Bronce argárico en el yacimiento de Fuente Álamo (Cuevas del Almanzora, Almería)
Luciañez Triviño M (2018) El marfil en la edad del cobre de la península ibérica. Una aproximación tecnológica, experimental y contextual a las colecciones ebúrneas del mega-sitio de Valencina de la Concepción - Castilleja de Guzmán (Sevilla). Dissertation, Universidad de Sevilla, Sevilla
Morillo León JM (2019) Exotic raw materials in Early Minoan Crete: cretaceous ivory and amber. A connection with the Iberian Peninsula? Proceedings of the 12th International Congress of Cretan Studies:1-12. Heraklion.
Morillo León JM (2021) Warriors, Sailors and Traders across the Sea: A glimpse on Mediterranean islands in the III millennium BC and the Bell Beakers phenomenon. In: Albertazzi G, Muti G, Saggio A, ISLANDS IN DIALOGUE (ISLANDIA). Proceedings of the First International Postgraduate Conference in the Prehistory and Protohistory of the Mediterranean Islands (Artemide)
Morillo León JM, Pau C, Guilaine J (2018) The proboscidean ivory adornments from the Hypogem of Padru Jossu (Sanluri, Sardinia, Italy) and the Mediterranean Bell Beaker. Zephyrus, LXXXII. https://doi.org/10.14201/zephyrus2018823563
Müller J, Van Willigen S (2001) New radiocarbon evidence for European Bell Beakers and the consequences for the diffusion of the Bell Beaker Phenomenon. In Nicolis F (ed.): Bell Beakers today. Pottery, people, culture, and symbols in Prehistoric Europe. Trento: Servicio Beni Culturali, pp. 59-80
Needham S (2009) Encompassing the Sea: ‘Maritories’ and Bronze Age maritime interactions. In P. Clark (Ed.), Bronze Age Connections: Cultural Contact in Prehistoric Europe (pp. 12–37). Oxbow Books. http://www.jstor.org/stable/j.ctt1cd0p3m.5
Nicoletti F (2001) Il commercio preistorico protostorico dell'ossidiana nel Mediterraneo ed il ruolo di Lipari e Pantelleria nel più antico sistema di scambio. In: Preistoria: Dalle coste della Sicilia alle Isole Flegree: Saggi. Palermo: Arnaldo Lombardi Editore, pp 259-270
Nocete F, Vargas JM, Schuhmacher TX, Banerjee A, Dindorf W (2013) The ivory workshop of Valencina de la Concepción (Seville, Spain) and the identification of ivory from Asian elephant on the Iberian Peninsula in the first half of the 3rd millennium BC. Journal of Archaeological Science 40:1579-1592. https://doi.org/10.1016/j.jas.2012.10.028
Novaline Jaga RM, Sundaram A, Natarajan T (1993) Wasteland development using geographic information system techniques. International Journal of Remote Sensing 14:3249e3257.
Orton CR (1982) Stochastic process and archaeological mechanism in spatial analysis. J Archaeol Sci 9 (1). https://doi.org/10.1016/0305-4403(82)90002-4
Pascual Benito JL (2012) El Taller de marfil del Bronce Pleno de la Mola d'Agres (Alicante). Iberia Archaeologica Band 16, Faszikel 1:173-198.
Pau C, Morillo León JM, Cámara Serrano JA, Molina González F (2018) Los objetos de adorno en marfil del yacimiento del Cerro de la Virgen (Orce, Granada). Complutum 29(2). https://doi.org/10.5209/CMPL.62581
Pebesma E (2018) Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal 10 (1):439-446. https://doi.org/10.32614/RJ-2018-009
Pebesma E, Bivand R (2023) Spatial Data Science: With Applications in R (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016
Pfälzner P (2016) The Elephants of the Orontes, Syria [Online], IV. http://journals.openedition.org/syria/5002 https://doi.org/10.4000/syria.5002
Renfrew C (1967) Colonialism and megalithismus. Antiquity 41:276–288.
Santana A, Arcos T (2009) La descripción de África en Plinio (Nat. V y VI). In: Candau JM, Gonzalez FJ, Chavez AL (eds), Libyae Lustrare Extrema. Realidad y literatura en la visión grecorromana de África: estudios en Honor del Profesor Jehan Desanges. Universidad de Sevilla, pp. 347-360.
Schuhmacher TX (2011) Elfenbein des Chalkolithikums und der Frühen Bronzezeit auf der iberischen Halbinsel: Ergebnisse einer interdisziplinären Forschungsprojekts. Elfenbein und Archäologie, INCETIVS Tagungsbeiträge 2004-2007, A. Banerjee, C. Eckmann, eds., RGZM 7, Mainz: 91-122.
Schuhmacher TX (2014) Dancing in the Dark? The westernmost "Cypriot" knot-headed pin from Aïn Smene (Morocco). Cuadernos de Prehistoria y Arqueología de la Universidad Autónoma de Madrid 40:41-50. http://doi.org/10.15366/cupauam2014.40.003
Siret L (1913). Questions de chronologie et d'ethnographie ibériques. Tome I: De la fin du quaternaire à la fin du bronze. Préface d'Émile Cartailhac. Journal des savants 11:458-459. www.persee.fr/doc/jds_0021-8103_1913_num_11_10_4134_t1_0458_0000_3
Souville G (1986) Témoignages sur l'âge du bronze au Maghreb occidental. In: Comptes rendus des séances de l'Académie des Inscriptions et Belles-Lettres, 130ᵉ année 1:97-114. https://doi.org/10.3406/crai.1986.14346
van Lanen RJ, Kosian MC, Groenewoudt BJ, Spek T, Jansma E (2015) Best travel options: Modelling Roman and early-medieval routes in the Netherlands using a multi-proxy approach, Journal of Archaeological Science: Reports 3:144-159. https://doi.org/10.1016/j.jasrep.2015.05.024
Verhagen P, Jeneson K (2012) A Roman Puzzle: Trying to Find the Via Belgica with GIS. In: Chrysanthi C, Murrieta Flores P, Papadopoulus C (eds) Thinking BeyondtheTool: Archaeological Computing and the Interpretive Process. Oxford: Archaeopress, pp. 123–130
Van de Noort R (2012) Exploring agency behind the Beaker phenomenon. The navigator's tale. In: Fokkens H, Nicolis F (eds), Background to the Beakers: Inquiries into Regional Cultural Backgrounds of the Bell Beakers Complex. Leiden: Sidestone Press, pp. 61–79
Verhagen P, Brughmans T, Nuninger L, Bertoncello F (2013) The Long and Winding Road: Combining Least Cost Paths and Network Analysis Techniques for Settlement Location Analysis and Predictive Modelling. In: Earl G, Sly T, Chrysanthi A, Murrieta-Flores P, Papadopoulos C, Romanowska I, Wheatley D (eds), CAA2012 Proceedings of the 40th Conference in Computer Applications and Quantitative Methods in Archaeology, Southampton, United Kingdom, 26-30 March 2012. Amsterdam: Pallas Publications, pp. 357-366
Wickham H, François R, Henry L, Müller K, Vaughan D (2023) dplyr: A Grammar of Data Manipulation. R package version 1.1.2, https://CRAN.R-project.org/package=dplyr

No competing interests reported.

tableivorieswithreferences.csv

Download PDF

Journal Publication

published 20 Jul, 2024

Read the published version in Archaeological and Anthropological Sciences →

Editorial decision: Revision requested
02 Feb, 2024
Reviews received at journal
12 Dec, 2023
Reviewers agreed at journal
18 Nov, 2023
Reviewers agreed at journal
18 Oct, 2023
Reviewers invited by journal
18 Oct, 2023
Editor assigned by journal
10 Oct, 2023
Submission checks completed at journal
05 Oct, 2023
First submitted to journal
04 Oct, 2023

You are reading this latest preprint version

Bronze Age communication and trade networks from Mediterranean Ivory records

Status:

Journal Publication

Version 1

Abstract

Figures

1 Introduction

2 Material and method

2.1 Archaeological background and site database

2.2 Environmental datasets

2.3 Model set-up

2.4 Network and prediction model of Ivory sites

3 Results and Discussion

3.1 Prehistoric ivory network model

3.2 Least Cost Path density

3.3 Predictive modelling results

3.4 Archaeological contextualization of the model

4 Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1