New Strategies in Archaeometric Provenance Analyses of Volcanic Rock Grinding Stones

1. Introduction

Among the vast field of archaeometry, provenance studies stand out. These studies involve the analysis of artifacts, raw materials, and geological sources to determine the origin and movement of ancient objects. Provenance studies help unravel the mysteries of human history, trade, and cultural exchange. These studies have been developed on a number of materials, e.g., pottery (and their corresponding clays), stones, metals, mortars, glasses, glazes and pigments. The strategy to trace the source of the corresponding raw materials is specific for every type of material. Importance can be focused on specific mineral inclusions, relevant chemical compositions, isotopic ratios, etc.

In the particular field of stone provenance studies, there is also a large diversity of approaches that simply reflects the variety of types of stones. Marbles and volcanic rocks are possibly among the most provenanced materials both for archaeological and analytical reasons. Obsidians (volcanic glasses) aside, volcanic stones had specific uses as grinding tools like millstones and mortars [1,2,3] or as building materials such as vaults [4,5]. The number of plausible volcanic sources can be delimited by using geological, geographical and historical data. In some cases, the analyzed materials can bear a particular chemical and mineral content enabling to trace them back to a specific origin among the possible sources. These provenance studies are of capital importance to infer trading routes and to assess the status of a given archaeological site.

In the present paper we characterize petrographically, and geochemically volcanic stones retrieved from two Roman archaeological sites and we explore and compare conventional and new strategies to infer their corresponding volcanic sources. On the one hand, Iulia Libica, a roman municipium (2nd-3rd century CE) located in the Eastern Pyrenees (between Spain and France). On the other hand, Sidi Zahruni, a non-excavated late Roman large pottery center (5th-7th century CE) located in northern Tunisia (Figure 1a). Particular attention has been given to quantify petrographic parameters. The two sites have been selected for their role as trading centers [6,7] and their proximity to a known geological source of volcanic stone (Olot Volcanic Field and Pantelleria island respectively). Unsupervised and supervised statistical methods have been applied to the geochemical data to tackle the provenance of the characterized materials.

2. Materials and Methods

2.1. Archaeological Materials

2.1.1. Three Roman Millstones from Iulia Libica (Eastern Pyrenees)

The Roman city of Iulia Libica (currently Llívia) was placed in the most important route crossing the Eastern Pyrenees, at more than 1200 m high, its most splendorous period was along the 1st and 2nd centuries CE. Archaeological research in the area started in the nineteen-seventies and regular excavation campaigns have been carried out since 1997. From 2013 one of the excavated areas has been confirmed as the central square (forum) of the Roman city [8]. In the 2019 excavation campaign works focused on the external part of the eastern wall of the forum reaching its foundations. Within the infilling material of a ditch that runs parallel to the eastern wall, abundant construction materials (tegulae and imbrices) and pottery were retrieved, the last indicating a chronology around the change of era. Also, within this ditch, three rotary millstones were retrieved along with a terracotta artefact interpreted as a small sundial [9].

The three millstones (Figure 1b) correspond to the catillus, the biconical funnel that rotates on a fixed bell-shaped millstone (meta) in hourglass-type millstones. The morphology of the millstones corresponds to type Lattes B2g dated from the last decades of 1st century BCE to the late 1st CE [10]. Most of the similar millstones found at Lattes and other sites in southern France were produced from local basalts.

Small chips from every retrieved millstone were carefully collected from fracture surfaces using a small hammer and a chisel. The three samples (labelled L1, L2 and L3) were split into two parts, one to be used for thin-section preparation, the other one to be pulverized for chemical analysis.

2.1.2. Volcanic Stone from Sidi Zahruni (Nabeul, NE Tunisia)

Sidi Zahruni is a non-excavated site located near the modern town of Beni Khiar, 6 km north-east of Nabeul (ancient Neapolis), one of the best-documented amphora production areas in roman Africa Proconsularis. The site covers an area of 13 ha and hosts a huge pottery where amphorae, in particular African Keay 25 type, were massively produced [11] as attested by the tilled fields containing an excessive number of pottery shards. In 2012, a systematic surficial survey was performed on the site to determine the different levels of shard concentration and their corresponding identifiable pottery typologies [12]. The chronology of the site, largely based on the known periods of production of the identified pottery typologies can be placed within the 5th-7th centuries CE.

Apart from pottery shards, marble, volcanic scoria, mosaic and cocciopesto mortar fragments were also occasionally found. These were interpreted as evidence of at least two areas with luxurious houses. Regarding the vesiculated volcanic scoria blocks, it is clear that they are not local stones but part of the archaeological site. The larger volcanic blocks on site were fragmented and irregular, the larger ones with an approximate volume of 3 dm³. A curved surface on some of them indicated a probable use as millstone, though other uses cannot be excluded. Pottery fragments including small embedded volcanic chips were also found indicating their use as mortarium.

Seven small samples of highly vesiculated volcanic scoria (labelled Z1 to Z7) were obtained in the form of irregular fragments from the largest blocks found on the surface of the site using a small hammer and a chisel. Similarly to the samples from Iulia Libica, the seven samples from Sidi Zahruni were split into two parts, to be used for thin-section preparation and chemical analysis respectively (Figure 1c).

2.2. Characterization Methods

Standard thin sections of the samples were prepared at the Laboratori de Preparació de Làmines Primes, Universitat Autònoma de Barcelona (UAB). Petrographic descriptions were performed using a Nikon Eclipse E400POL microscope; pictures in plane polarized light (PPL) and cross polarized light (XPS) modes were taken using a DS-Fi2 High-Definition color camera head mounted on the microscope. Thin section observation was used to determine the fabric, to identify the mineralogy of phenocrysts and to describe their shapes and features. High-resolution (4800 ppp) full scans of the analyzed thin sections were obtained using an EPSON 4990 scanner, including both PPL and XPL modes. The scans were treated with image processing software ImageJ [13] to quantify porosity, matrix and phenocryst ratios including discrimination between the different minerals that appear as phenocrysts.

Petrographic information was complemented with elemental point-composition measurements on the phenocrysts imaged using a Merlin Field Emission Scanning Electron Microscope (FE-SEM) (Zeiss, Oberkochen, Germany) equipped with an energy-dispersive X-ray spectrometer (EDS) at the Servei de Microscòpia (UAB). Average elemental composition over different areas of the microcrystalline matrices were also measured. Measuring conditions were accelerating voltage of 15 kV, 1 nA beam current and an overall counting time of 60–100 s.

Multielemental bulk analysis was accomplished by means of energy dispersive X-ray fluorescence (EDXRF) for the samples with remaining powder after cutting the thin section. For sample preparation, 5 g of previously powdered sample were mixed with 0.4 g of a binding agent (Elvacite™) and then homogenized in an agate mortar [14]. The mixture was pressed at 10 tons to achieve a homogeneous pressed powder pellet of 4 cm in diameter. Different pellets were produced for the three samples from Iulia Libica and five of the Sidi Zahruni samples (Z1, Z2, Z4, Z5, Z7). The pellets were analyzed using a benchtop spectrometer S2 Ranger EDXRF system (Bruker AXS, GmbH, Germany) with a Pd X-ray tube (max. power 50W) and XFLASHTM Silicon Drift Detector (SDD) with a resolution of <129eV Mn-Kα1 at 100,000 cps. Samples were evaluated under vacuum at different excitation voltage conditions to properly excite elements of low, medium and high atomic number. Multielemental quantitative determination for Si, Al, Fe, Mg, Ca, Na, K, Ti, Mn, Nd, Cr, Ni, Zn, Rb, Sr, Zr, Ba, Y and Pb was made using a calibration designed for volcanic rocks and operated with the SPECTRA.EDX package (Bruker AXS, GmbH, Karlsruhe, Germany). Nine certified reference materials (CRMs) were used comprising several basalts (BR –CRPG, BE-N –IWG-GIT, BHVO2 –USGS, DC73303 –NCS), an andesite (AGV-1 –USGS), a diorite (DR-N –ANRT), a dolerite (WSE –IWG-GIT), a microgabbro (PM-S –IWG-GIT) and a syenite (STM-1–USGS).

Reference geochemical data corresponding to volcanic materials were extracted from the GEOROC database [15]. A starting dataset was obtained by setting a query by geography (latitude from 25°N to 54°N and longitude from 12°W to 45°E) to cover the Mediterranean Basin, North-Africa and a great part of Europe and then selecting only the records that declare the rock name as basalt, excluding modern eruptions. This produced a database of around 3600 records that was additionally restricted to chemical analyses that really correspond to the basalt area in a total alkali-silica (TAS) diagram [16] (with data recalculated on an anhydrous basis). Data regarding Fe was standardized recalculating all the analyses to total iron as Fe₂O₃. This left 1855 records that were further reduced to smaller homogeneous datasets containing a given set of analyzed chemical elements. This produced datasets G90, G60 and G30 containing roughly 90%, 60% and 30% of the original records. The data was labelled into different groups according to the corresponding tectonic setting and location (Figure 2). The original database (Table S1) and the three reduced datasets (Tables S2–S4) can be downloaded from the Supplementary Materials section.

Several statistical methods were applied combining both the geological reference geochemical values as well as the data obtained from the analyzed archaeological materials. Principal component analysis (PCA) was used to check the structure of the data, in particular to investigate if any clustering appears using the main components to distinguish different reference fields. The data was standardized to zero mean and unit variance before computing the PCA. A number of supervised methods were applied to the reference datasets, where all the compositional values were divided by the SiO₂ concentration as a way to standardize them. Different 8:2 splits were tested following a train-test approach. In a second step the trained models were used to obtain the predictions of provenance for the archaeological samples. The supervised approach was performed using the “Supervised Provenance Analysis” R Markdown files from Anglisano et al. [17] including different classification models: generalized linear models (GLM), random forest (RF), artificial neural network (ANN), weighted k-nearest neighbors (kkNN), linear discriminant analyses (LDA) and a stack of all these models as described in [18]. Four additional classification models were tested using Python code (see Supplementary Materials section) under the Jupyter Notebook interactive computing platform: gradient boosting (GB), Gaussian process classifier (GPC), Gaussian naive Bayes (GNB) and linear support vector machine (LSVM).

3. Results

3.1. Petrography

The three millstone samples from Iulia Libica exhibit a porphyritic and vesiculated texture with rather idiomorphic phenocrysts of two compositions (Figure 3a). The most abundant and larger phenocrysts correspond to clinopyroxene (cpx, augite as analyzed by SEM-EDS, see Figure 4 and Table 1), often appearing in clusters of polygonal idiomorphic crystals, some of them exhibiting occasional hourglass zoning (Figure 3d).

Basal sections reveal two cleavage directions at near 90° angles (Figure 3c). Smaller phenocrysts of olivine (ol) are also very common in the three samples. These are less idiomorphic (often rounded and occasionally corroded and dissolved into the groundmass) and characteristically with outer alteration rims containing iron oxide that surround the fresh olivine cores, indicating post-magmatic hydrothermal alteration processes (Figure 3b,c). The fresh olivine composition measured by SEM-EDS is Fo_75±2, while the composition in the altered rims is still essentially olivine but with an increased amount of Fe (19-29 wt % instead of 18 wt %). Small amounts (< 1 wt %) of Ca and Al are also occasionally detected in the altered rims. The groundmass is a fine-grained matrix formed by unoriented lath-shaped feldspar crystals, clinopyroxene and opaque equidimensional minerals (possibly Ti-bearing ilmenite or magnetite).

All the collected samples from Sidi Zahruni exhibit also a porphyritic and vesiculated texture with more or less idiomorphic phenocrysts (Figure 5) and a uniform mineral association. The identified minerals in form of phenocrysts are: i) zoned and twinned feldspars with rectangular and quadratic sections (Figure 5a). SEM-EDS analyses indicate that they correspond to plagioclases (pl), with composition ranging from 46-78 An% (labradorite and bytownite); ii) clinopyroxene (cpx) crystals with a slight brown colour (augite as analyzed by SEM-EDS, see Table 1 and Figure 4). These appear particularly idiomorphic, exhibiting mainly hexagonal sections or octagonal basal sections (Figure 5f), and occasionally showing the diagnostic two cleavage directions at near 90° angles. Hourglass zoning also appears in some crystals (Figure 5h); iii) olivine (ol), with rare polygonal sections and more common rounded shapes (Figure 5e), sometimes significantly corroded and partially dissolved into the groundmass. Its composition by SEM-EDS is Fo_72±4. Apart from these, in one of the samples (Z6) a large green-coloured (PPL) rounded crystal was spotted (Figure 5g). The SEM-EDS analysis revealed a clinopyroxene composition with increased amounts of Na and Fe indicating a member of the aegirine-augite solid solution (see Table 1). The different phenocrysts are embedded also in a fine-grained matrix mainly formed by lath-shaped feldspar crystals, clinopyroxene and opaque equidimensional minerals (possibly Ti-bearing ilmenite or magnetite).

The quantification of the petrographic components was performed on treated scan images of the petrographic thin sections (Figure 6) allowing the discrimination between porosity, matrix and the different types of phenocrysts. Porosity is around 24-39% of the analyzed area for Iulia Libica (L) samples and usually lower (6-30%) for Sidi Zahruni (Z) samples, although for sample Z5 exceeds 50%. Pore size histograms reveal that small pores (with areas below 1 mm²) are the most abundant in all the analyzed samples (Figure 7). The three samples from Iulia Libica display a wide range of pore sizes similarly to some of the samples from Sidi Zahruni. However, samples Z3 and Z4 have narrower pore size distributions. Circularity of the pores is within the range 0.9 to 0.5 for all the samples and the mean circularity values are around 0.7 or a bit lower for the samples with irregular or elliptic pore shapes (L2, Z2, Z6 ad Z7).

Regarding the mineral part of the samples, it is mainly formed by the fine-grained matrix in all the samples and the phenocrysts are only a minor component (3-16%), see Table 2. The proportion of matrix is higher for the samples from Sidi Zahruni (constituting 94% ± 2) compared to Iulia Libica samples (88% ± 4).

The only two types of phenocrysts that appear in Iulia Libica samples have a relative ratio of around 2:1 between clinopyroxenes (cpx) and olivine (ol). In contrast, in the basalts from Sidi Zahruni the abundance of these components is reversed and surpassed by plagioclase (pl) phenocrysts (Table 2). The pl:ol:cpx relative ratio in samples from Sidi Zahruni is around 4:1.5:1. The quantification of the crystals indicates that regardless of the sample and considered mineral, the frequency of sizes increases towards small values blending with those forming the microcrystalline matrix (Figure 8). The area of the phenocrysts only reaches rarely values above 1 mm².

The distribution of clinopyroxene and olivine crystals is less biased towards small crystals for the basalts retrieved from Iulia Libica, while the distribution corresponding to clinopyroxene exhibits a secondary peak centered around 0.7 mm² for these basalts. The distribution corresponding to olivine is relatively similar for all the samples and there are also (regardless of the archaeological site provenance) hints of a secondary peak for sizes around 0.4 mm². The crystal size of feldspar crystals (only present as phenocrysts in samples from Sidi Zahruni) is possibly the one that fits better to a monotonous exponential decay (to a second order exponential decay function to be precise).

Another quantified parameter (circularity) indicates higher mean values (<C>) for clinopyroxenes and olivine crystals from the basalts retrieved from Sidi Zahruni. However, feldspar phenocrysts from these basalts deviate markedly from circularity consistently with their common lath-shaped morphologies (see <C> in Figure 8).

3.2. Geochemistry

The average bulk composition of the analyzed basalts is shown in Table 3. Each set of samples have a homogeneous composition in agreement with the observed petrographic texture and mineralogical composition. All the analyzed samples have a similar amount of silica (~46-48 %). In contrast, there are slight but significant differences between L and Z samples regarding the concentration of the other elements. For instance, L samples are richer in K, Fe and P and poorer in Al, Ca, Na and Ti. However, the sum of alkali oxides (Na₂O + K₂O) is similar for both types of samples so that all the investigated samples can be classified as basalts (Figure 9) according to the TAS diagram [16]. The higher concentration of Al, Ca and Na in Z samples is possibly connected to the presence of plagioclase (Ca-Na aluminosilicates) as phenocrysts only in the samples retrieved from Sidi Zahruni. In contrast, the amount of Fe is slightly higher in the samples from Iulia Libica (~14% compared to ~12%) correlating with the observed Fe-rich alteration rimming in olivine phenocrysts. Finally, the trace element proportions are also dissimilar for both types of basalts. Those retrieved in Iulia Libia are richer in trace elements, except the amounts of V that are higher in the basalts from Sidi Zahruni.

The tectonic setting of volcanic lavas can be sometimes geochemically inferred using certain compositional ratios (e.g., Zr vs. TiO₂ or Zr vs. Zr/Y) [19,20]. The measured geochemical values indicate that the samples (both from Iulia Libica and Sidi Zahruni) correspond to intraplate basalts and this should help to determine the provenance of the samples.

3.3. Geochemical Comparison with Reference Materials

3.3.1. Elemental and PCA Biplots

Provenance determination using geochemical data is commonly undertaken using two-variable scatterplots that examine the relations between two chemical elements (or two simple chemical ratios). In the same plot several reference variation fields are also depicted. The choice of the two elements (or chemical ratios) is arbitrary and usually the only criterion is to find a combination of variables that allows discrimination between the different possible provenances.

In Figure 10 several examples of biplots are shown. When only a few possible provenances are considered (Figure 10a) and the analyzed samples lie within a given reference field, the provenance seems well determined. For instance, looking at Figure 10a the geological provenance of the samples retrieved from Sidi Zahruni appear to be attributable to the Pantelleria island (strait of Sicily) and those from Iulia Libica to the Hyblean Plateau (southeastern Sicily). In contrast, if the analyzed samples do not lie within any reference field it can be assumed that the samples do not belong to any of the considered provenances or alternatively that the reference fields perhaps are not well defined due to a number of potential reasons like a lack of a fully representative set of reference samples. For instance, one of the samples from Sidi Zahruni is not strictly within any reference field.

When the number of possible provenances is increased the degree of overlap between reference fields increases inevitably and it becomes very unlikely that the analyzed samples could be unambiguously assigned to a given provenance. Looking at Figure 10b, the basalts retrieved from Sidi Zahruni could have been quarried in Pantelleria but also in Sardinia and those from Iulia Libica, apart from the suggested Hyblean Plateau, could also belong to the Southern France basalts or to those from the Middle Atlas (Morocco). Also, using different biplots dissimilar provenances are sometimes suggested. Using the biplot from Figure 10c, the samples from Sidi Zahruni appear more scattered, within the reference fields of Pantelleria and partially within the Hyblean Plateau and Spanish basalts reference fields. On the other hand, the samples from Iulia Libica are now also within the reference fields of Middle Atlas, Southern France and Spain basalts but, strictly, not within the Hyblean Plateau reference field.

Using two-variable plots, a large number of compositional parameters are disregarded. Multivariate statistical analyses extend the examination of the data to all the analyzed elements. In particular, PCA condense the relevant information in a few principal components that can be represented also in biplots. However, one major difficulty is to have a uniform set of data. The reference data from GEOROC database contains heterogeneous data, since samples were prepared using different methods and measured using different equipment. The analyzed elements are not always the same set, in particular, the minor and trace elements can vary.

Two PCA were computed using two different sets of data. For the PCA shown in Figure 11, a total 1702 records (G90 dataset) from the preselected GEOROC database were combined with the 25 individual analyses for samples from Sidi Zahruni and 8 for samples from Iulia Libica. To get a homogeneous dataset the number of considered compositional variables was restricted only to nine (SiO₂, TiO₂, Al₂O₃, CaO, MgO, MnO, K₂O, Na₂O and P₂O₅), all of them being main components of the basalt samples (both geological reference and archaeological materials). This allowed to consider samples from 26 different sites/geological settings. The first two principal components have a joint explained variance ratio of ~50%. The obtained PCA reveals a positive correlation between Si and Al, and a negative correlation of them with Mn, P and Ti. These variables are the main contributors to the first principal component (PC1). On the other hand, there is also a very good positive correlation between K and Ca, and both are negatively correlated with Na. Both Sidi Zahruni and Iulia Libica appear in a similar area of the biplot with slightly negative PC1 values. Concerning the PC2 values, they are slightly positive for Sidi Zahruni samples and almost zero or slightly negative for Iulia Libica samples. However, the degree of overlap between different provenances observed in the two main components biplot is very high, therefore is not possible to obtain clear data about the origin of the analyzed samples. The samples from the most represented provenance (Anatolia, see yellow dots in Figure 11) spread covering almost all the point cloud formed by all the dataset. Reference samples that are geographically close to Iulia Libica (L), such as those labelled as Catalonia, Calatrava and Massif Central appear relatively close to the measured specimens of L samples but individual samples from many other reference sites plot in the same area. Instead of assigning provenances perhaps appears more clearly that some of the considered reference provenances can be excluded. Specifically, the Aeolian arc, the Campanian and the Lazio basalts, seem to bear compositions that have not affinity with those sampled at Sidi Zahruni and Iulia Libica. The Aegean arc, Algeria, Alps-Sardinia-Corsica and Etna provenances could also be discarded, except for some occasional individual samples.

An additional PCA was computed (Figure 12) with a much more complete set of compositional variables (20 including, besides the previously considered, Fe₂O₃, Cr, V, Ni, Cu, Zn, Rb, Sr, Y, Zr and Nb). In this case, the preselected GEOROC database was reduced to 550 records from 18 different geological sites constituting dataset G30. In this newly computed PCA the degree of overlap between different classes is not so strong. Besides, the samples from Sidi Zahruni (Z) and Iulia Libica (L) appear in a marginal area of the point cloud and therefore some provenances can be discarded with a certain confidence, especially for Iulia Libica basalts, which appear with higher PC1 values (see X symbols in Figure 12). Namely, the Aegean, Gibraltar and Tyrrhenian arcs form clusters without intersecting neither the Z or L clusters and from the Aeolian arc samples only one is in the area of Z and L samples; the Campanian, Etna, Liguria and North Anatolia samples are also far from the Z and L clusters but these provenances contain very few samples in the processed database. Considering that Z and L clusters appear relatively close, certain provenances match with the position of both clusters, that is the case of Catalonia, Egypt, Atlas and Sicily Channel rift clusters. Additionally, the samples from Suez match particularly well with those from Sidi Zahruni. Besides, other provenances cannot be completely excluded. Also, there are a number of issues that should be considered: the joint PC1 and PC2 explained variance ratio is only ~44%, some reference sites contain very few samples within the processed database and they could not be statistically representative of their location, some actually have disappeared from the PCA because GEOROC does not contain any analyzed sample including all the considered compositional variables.

3.3.2. Supervised Classification Methods

The supervised approach was applied to the reference data obtained for the GEOROC database. Three different databases were tested (G90, G60 and G30 as described in Section 2.2). From the different classes (i.e., provenances) defined for the data in these databases, only those bearing at least 9 samples were considered. Classes with a smaller number of features could not be properly trained and predicted by the machine learning algorithms. As one of the classes relevant to the present study contained very few samples, additional basalt analyses [24,25] were added to the database to keep the corresponding provenance into consideration. This class (X8) is labelled as Catalonia and corresponds to the Olot volcanic field, only at around 50 km south-east of Iulia Libica. An additional reason to retain this class is that the corresponding provenance was actually suggested from the results of some of the computed PCA. Finally, 18, 16 and 12 different classes were considered from G90, G60 and G30 datasets respectively.

Table 4, Table 5 and Table 6 present the results corresponding to individual runs of the supervised classification models using G90, G60 and G30 datasets respectively. The tables display the classification accuracy of three different random test sets for each trained classification model and the overall classification predictions for the ensemble of Z samples and L samples. The results appear very consistently homogeneous and pointing to a given provenance for Z samples, in contrast for L samples the results are less monotonous, and a variety of provenances is suggested.

The different models tested, regardless of the database used, indicate a very high probability that the samples from Sidi Zahruni (Z) are from the Sicily Channel rift (i.e., Pantelleria island) which is not far from the archaeological site (~115 kilometers). Using the G90 database and combining all the Z samples and all the trained models, the overall probability for class X21 (Sicily Channel rift) amounts 88.4% and only two more classes are suggested, Anatolia (X5) with a 11.1% and Egypt (X10) with only a 0.5%). Using the G60 database only two provenances are suggested, X21 captures an overall probability of 96% and X5 the remaining 4%. Finally, using G30 class X21 holds again almost all the probability (96.4%) and the rest is split between Anatolia (X5) with a 3.3% and Suez (0.3%).

Contrastingly, the different models do not show a similar agreement to assign a provenance to the samples from Iulia Libica (L) and there are even significant differences using different splits (see for instance L-sample predictions using the ANN model in Table 4). In any case, using G90 database and combining all the results, the probability spreads into eight different provenances including Anatolia (X5), 50.4%; Libya (X15), 17.5%; Adria Domain (X1), 7.9%; Hyblean Plateau (X13), 7.1%; Catalonia (X8), 5.8%; Atlas (X6), 5.0%, Massif Central (X17), 3.8% and Sicily channel rift (X21), 2.5%. One could think that at least one of the classes (X5) captures more that 50%, but checking out the results obtained using G60 (Table 5) the predominance of class X5 is lost. The overall results using G60 database indicate even more fractionated probabilities: Atlas (37.1%); Hyblean Plateau (25.4%); Anatolia (17.1%); Sicily Channel Rift (7.9%); Catalonia (6.3%); Egypt (5.4%); Adria Domain and Libya (both 0.4%). Finally, using the database containing the highest number of compositional variables (G30), the overall probabilities are spread again among different provenances, in this case: Atlas (30.8%); Anatolia (19.2%); Hyblean Plateau (17.1%); Catalonia (9.2%); Sicily Channel Rift (3.3%) and Egypt (0.8%).

4. Discussion

Regardless of the properties or techniques used, the attribution of provenance to archaeological volcanic samples requires an accurate quantification of a set of parameters, and also a dataset containing the same quantification for a representative collection of reference geological samples.

In this paper a quantification effort has been done both for petrographic and compositional properties. Some of the quantifiable petrographic properties could arguably not be representative of a given provenance, e.g., the porosity or the shape of the pores could vary enormously even within the rocks of a given volcanic eruption. However, the petrographic characterization has possibly the potential to assembly a set of discerning elements including mineralogy of the phenocrysts, relative ratios, size, morphology, presence of alterations, texture/mineralogy of the matrix, etc. The problem with petrography is that presently there is a lack of large reference datasets that should comprise a standard set of quantified parameters. Therefore, the petrographic approach is not ready to fulfill the requirements to be applied as a routine technique for provenance determination. However, it has the potential to be developed as such. In any case for well delimited cases such as a binary classification problem, the petrographic approach can be applied because the self-production of the required reference dataset is feasible. Besides this, the petrographic characterization is also very helpful to determine whether the set of sampled materials from a given archaeological site is homogeneous (likely from a single provenance) or not. It is worth to mention that a set of volcanic samples mineralogically and petrographically different could, by chance, exhibit chemical homogeneity.

In the present investigation, samples retrieved from Sidi Zahruni (Z) appear petrographically homogenous indicating a single geological supply of volcanic materials and the same can be said for the three millstones from Iulia Libica (L). The minerals that have been identified are those common in basalts such as clinopyroxenes, olivine and felspars, the later absent as phenocrysts in L samples. The absence of felspars (or the absence of leucite in both L and Z samples) could help to constrain the provenance. Other particular petrographic features that could be helpful are the alteration rims in olivine crystals (L samples) or the occasional presence of green aegirine-augite crystals (Z samples). Alteration rims in olivine is actually a common feature in many potential provenances, including Sardinia, Agde (the south end of Massif Central), Olot (Catalonia), Ustica (Italy) and Middle Atlas (Morocco) basalts [26]. The presence of aegirine-augite crystals has been described in Pantescan volcanic rocks and is typical in their more felsic rocks [27] sometimes intimately associated [28].

The chemical approach is currently much more ready than the petrographic approach to determine the provenance of archaeological volcanic materials. The quantification of the chemical composition is routinely applied by many researchers that publish their data and some initiatives have been undertaken to collect systematically these data to conform huge sets of geochemical reference databases like GEOROC. However, some difficulties arise from the use of different equipment and sample preparations methods by the contributors to the reference database. Also, the geochemical data can be expressed in different ways (as oxides or elements, wt%, ppm, including volatiles or not, assuming a given valence state for elements like Fe or distinguishing different states, etc.).

Despite carefully selecting and standardizing relevant sets of data, the limitations of the common geochemical classification tools have been illustrated. On the one hand biplots representing only two chemical elements (or two chemical elemental ratios) can be misleading because they disregard data that could be relevant. On the other hand, multivariate data analysis methods like PCA, despite considering all the features (i.e., the analyzed chemical elements), reveal that as the number of considered predefined classes increases it is more likely that they form overlapping clusters. For the studied archaeological samples retrieved from Iulia Libica (L) the suggested provenances using certain biplots are Hyblean plateau, Middle Atlas, Southern France and Spain. For those retrieved in Sidi Zahruni (Z) the suggested provenances are Pantelleria (recurrently), Sardinia and also the Hyblean Plateau. PCA is not really useful to decide which among these suggested provenances is the most probable. However, PCA can be used to discard some of the multiple considered provenances. The basalts from the Aeolian, Aegean, Gibraltar and Tyrrhenian arcs, the Campanian and those from Lazio and possibly Algeria, Alps-Sardinia-Corsica and Etna provenances could be discarded both for the samples from Iulia Libica and Sidi Zahruni.

Finally, the supervised machine learning models are also multivariate methods that take into consideration a high number of compositional variables, but they have the additional advantage of being designed to discriminate the different predefined reference classes (i.e., the provenances). The capacity to discriminate the different classes can be measured using different statistical indicators computed using the test sets. Along this line, as it was previously stated in a methodologically similar study [17], using supervised models it is very important to check for the convergence between the results using different classification models and training splits. Reliable provenance determination requires a significant agreement between the results using different strategies. In the present research, besides different models and splits, three different reference databases have been used (G90, G60 and G30) with average accuracies ranging from 0.55-0.80 (using G90) to 0.66-0.90 (using G30).

In the case of the samples retrieved from Sidi Zahruni (Z) the results point very monotonously at Sicily Channel Rift (X21 class) as the provenance of the samples. This means basalts extracted from Pantelleria or Linosa volcano islands, which are made of ocean island basalts known to be exploited in antiquity for the manufactures of lava grinding tools like millstones [3,29]. This origin (especially Pantelleria) had appeared recurrently among the provenances suggested by the common geochemical classification tools (biplots and PCA). Besides, the petrographic properties of Z samples also agree with those exhibited by Pantelleria basalts, e.g., weakly porphyritic texture, mineralogy and relative abundance of phenocrysts with pl > ol > cpx [29,30,31]. From the two basalt types described in Pantelleria [29] those with a higher affinity with the samples retrieved from Sidi Zahruni are the younger basalts with relatively low TiO₂ and P₂O₅ contents [30].

In contrast, for the samples from Iulia Libica there is no agreement between the models and usually, the predicted provenance appears statistically divided among different classes that hold fractions of provenance probability. Up to nine different classes (i.e., provenances) are suggested by the models and in general the probability average percentages combining all the predictions are very scattered (the highest value is 50.4% for class X5 (Anatolia) using database G90). According to [17] this heterogeneity of results is indicating that the provenance attribution is not robust. Indeed, in some cases even the results from a given classification algorithm vary significantly using different trained models (i.e., different train-test sets). Some classes (among others incidentally X5) are overrepresented within the reference databases and the class appears spread out covering most of the point cloud observed in the corresponding PCA. Therefore, it is not surprising that such classes could capture substantial fractions of the provenance probability. Using the geographical proximity as a criterion, classes X17 (Massif Central) and specially X8 (Catalonia) are the nearest locations to Iulia Libica. The models assign to these classes rather low probabilities and the two classes have in common a very low number of data within the reference databases. Even if training is enabled, a low number of data could not be fully statistically representative of the class. In fact, class X17 has only been considered by the models using database G90 because in the other datasets the number of data was too low to train the models to discriminate this class. The same would have occurred for class X8 if we had not supplemented the original database (produced using GEOROC) with additional data for this class. Ideally, all classes with a low number of reference samples should be supplemented. In particular, it would have been interesting to increase the reference samples for class X17 as the Massif Central (in particular cap d’Agde) is the presumed origin for the basaltic millstones produced in Lattes [10], typologically similar to those found in Iulia Libica. Such trading route would be consistent with presence of pottery from the southern Gaulish area in Iulia Libica [6]. With the available data and, comparing with the results predicted for Sidi Zahruni, the conclusion is that the provenance of Iulia Libica samples has not been successfully determined and this would imply that the provenance corresponds to a class not contained within the database or imperfectly represented within the considered reference databases.

To conclude, supervised machine learning methods have a number of advantages over other statistical methods. Conceptually, they are designed to learn the best way to discriminate the different classes. However, they are not foolproof, and they are subject to certain limitations linked to incomplete reference databases (unbalanced, underrepresented or absent classes). Despite the good indicators of overall performance (see accuracies in tables G90, G60 and G30) obtained after testing the trained models, a model will obviously never be able to predict a class not described in the training database and the performance of the model to predict underrepresented classes can be far from appropriate. For instance, from the confusion matrices obtained during the test step, it could be seen that classes X8 (Catalonia) and X17 (Massif Central) exhibit low average sensitivities per model (usually below 0.5 and sometimes even 0). Specifically, for class X8 the highest average sensitivity values for are 0.78 and 0.69 using LSVM and kkNN models respectively and incidentally these two models are the ones that statistically attribute higher provenance probabilities to class X8 compared to other models (see Table 4, Table 5 and Table 6). As for class X17 (only considered in runs using database G90) the highest average sensitivity is only 0.23, obtained using LSVM model and this is actually the only model that assigns a certain probability to this provenance.

5. Conclusions

Statistical approaches to establish provenance of volcanic rock tools require quantification of a set of properties and a large reference database. The quantified properties could be geochemical, mineralogical, or petrographic. However, currently only large reference geochemical datasets are available. To implement the statistical approach to petrographic properties would first require an agreement on the parameters to be quantified and a collective effort to produce the corresponding datasets.

The limitations of the common geochemical classification tools have been illustrated. On the one hand biplots using few variables can be misleading because relevant data could be omitted. As a result of this, different biplots can suggest different provenances. On the other hand, multivariate data analysis methods, like PCA, often exhibit a single point cloud and not separated clusters and then it is difficult to assign a given provenance to the unlabeled archaeological data.

The supervised machine learning approach applied to geochemical data has proven to be a powerful tool to discern the different reference clusters, and this enables provenance prediction for unlabeled samples. Among the different 18 considered provenances, the models assign systematically to all the samples from Sidi Zahruni a very high provenance probability to the class representing the basalts from the Sicily Channel rift. The provenance of Sidi Zahruni samples has been successfully established and not only because the strong agreement between the different supervised models but also due to the high petrographic consistency between the studied samples and the known characteristics of the Pantelleria basalts. Additionally, this provenance was among those suggested using common geochemical classification tools and it lies only ~115 kilometers from the archaeological site of Sidi Zahruni.

In contrast, the models assign different provenances to the samples from Iulia Libica and this would imply that the corresponding true provenance was not among the considered reference datasets or it was incompletely described. This highlights the limitations of the supervised approach and the need for a statistical approach including the use of different sister samples, different trained models and different classification algorithms to identify a homogeneous and therefore robust provenance prediction.