New Strategies in Archaeometric Provenance Analyses of Volcanic Rock Grinding Stones: Examples from Iulia Libica (Spain) and Sidi Zahruni (Tunisia)

Abstract

Archaeometry can help archaeologists in many ways, and one of the most common archaeometric objectives is provenance analysis. Volcanic rocks are often found in archaeological sites as materials used to make grinding tools such as millstones and mortars or as building materials. Petrographic characterization is commonly applied to identify their main mineralogical components. However, the provenance study of volcanic stones is usually undertaken by comparing geochemical data from reference outcrops using common descriptive statistical tools such as biplots of chemical elements, and occasionally, unsupervised multivariate data analysis like principal component analysis (PCA) is also used. Recently, the use of supervised classification methods has shown a superior performance in assigning provenance to archaeological samples. However, these methods require the use of reference databases for all the possible provenance classes in order to train the classification models. The existence of comprehensive collections of published geochemical analyses of igneous rocks enables the use of the supervised approach for the provenance determination of volcanic stones. In this paper, the provenance of volcanic grinding tools from two archaeological sites (Iulia Libica, Spain, and Sidi Zahruni, Tunisia) is attempted using data from the GEOROC database through unsupervised and supervised approaches. The materials from Sidi Zahruni have been identified as basalts from Pantelleria (Italy), and the agreement between the different supervised classification models tested is particularly conclusive. In contrast, the provenance of the materials from Iulia Libica remained undetermined. The results illustrate the advantages and limitations of all the examined methods.

1. Introduction

Among the vast field of archaeometry, provenance studies stand out [1]. These studies involve the analysis of artifacts, raw materials, and geological sources to determine the origin and movement of ancient objects [2]. Provenance studies help unravel the mysteries of human history, trade, and cultural exchange. This type of research has been developed on a number of materials (see [2] and references therein), including pottery (and their corresponding clays) [3,4], stones [5,6,7], metals (e.g., [8,9]), mortars/plasters (e.g., [10,11]), glasses [12], glazes (e.g., [13,14]) and pigments (e.g., [15,16]). The strategy to trace the source of the corresponding raw materials is specific to every type of material. Importance can be focused on specific mineral inclusions, relevant chemical compositions, isotopic ratios, etc. The most convenient and useful characterization techniques vary depending on the material under study and the tracing strategy.

In the particular field of stone provenance studies, there is also a large diversity of approaches that simply reflect the variety of rock types [17,18,19,20,21,22]. Marbles and volcanic rocks are probably among the materials whose provenance has been most frequently investigated, and this is for both archaeological and analytical reasons. Obsidians (volcanic glasses) aside, volcanic rocks had specific uses. They were often used as grinding tools like millstones and mortars [23,24,25] or as building materials such as vault stones [26,27]. 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 the ability to trace them back to a specific origin among the possible sources [28,29,30,31]. 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 petrographically and geochemically characterize volcanic stones retrieved from two Roman archaeological sites. We explore and compare conventional and new strategies to infer their corresponding volcanic sources. The novelty lies (i) in the exhaustive quantification of petrographic data as opposed to the common non- or poorly quantified petrographic descriptions and, (ii) more importantly, in the innovative use of supervised machine learning algorithms instead of the common unsupervised exploratory geochemical data analysis to tackle the provenance of the characterized materials. The studied materials were retrieved from Iulia Libica, a Roman municipium (2nd–3rd century CE) located in the Eastern Pyrenees (between Spain and France) and Sidi Zahruni, a non-excavated late Roman large pottery center (5th–7th century CE) located in northern Tunisia (Figure 1a). The two sites have been selected for their role as trading centers [32,33] and their proximity to a known geological source of volcanic stone (Olot Volcanic Field and Pantelleria island, respectively).

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 during the 1st and 2nd centuries CE. Archaeological research in the area started in the 1970s, and regular excavation campaigns have been carried out since 1997. Since 2013, one of the excavated areas has been confirmed as the central square (forum) of the Roman city [34]. In 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, with the latter indicating a chronology around the change of era. Also, within this ditch, three rotary millstones were retrieved along with a terracotta artifact interpreted as a small sundial [35].

The three millstones (Figure 1b) correspond to three thin catilli, the biconical funnel that rotates on a fixed bell-shaped millstone (meta) in rotary hand-mills. The morphology of the millstones corresponds to type Lattes B2g, dating from the last decades of 1st century BCE to the late 1st CE [36]. 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 (labeled L1, L2 and L3) were split into two parts, one to be used for thin-section preparation and 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 northeast 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 center where amphorae, in particular African Keay 25 type, were massively produced [37], 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 [38]. 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, they are clearly non-local stones but part of the archaeological site. The volcanic blocks on site were fragmented and irregular, with the larger ones having 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 (labeled 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 [39] to quantify porosity, matrix and phenocryst ratios, including discrimination between the different minerals that appear as phenocrysts. The modal relative proportions of the components were formally measured as area percentages. However, the relative surface measurements are representative of the volumetric relative proportions, considering that the analyzed samples are relatively isotropic and fine-grained. The different phenocryst compositions were analyzed separately using the “analyze particles” command within ImageJ to obtain size-frequency distributions. The crystals were fit to ellipses, and several shape parameters were computed; among others, circularity (C) defined as 4π·area/(perimeter)².

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 was 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 largest remaining sample parts after thin-section preparation. For EDXRF sample preparation, 5 g of previously powdered sample was mixed with 0.4 g of a binding agent (Elvacite™) and then homogenized in an agate mortar [40]. 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, Karlsruhe, Germany) with a Pd X-ray tube (max. power 50 W) and XFLASHTM Silicon Drift Detector (SDD) with a resolution of <129 eV 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, P, Mn, V, Cr, Ni, Cu, Zn, Rb, Sr, Y, Zr and Nb was made using a calibration designed for volcanic rocks and operated with the SPECTRA.EDX package (Bruker AXS, GmbH, Karlsruhe, Germany). These elements were always present above detection limits in all the analyzed samples. 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). Calibration lines indicate very good correlation (r² > 0.997) between measured and expected concentrations, with absolute errors ranging from 2 to 9% for all the elements. Achieved precision is higher than the intrinsic compositional variation of the analyzed rocks.

Reference geochemical data corresponding to volcanic materials were extracted from the GEOROC database [41]. 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 [42] (with data recalculated on an anhydrous basis). Data regarding Fe were standardized, and all the analyses were recalculated 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 were labeled into different classes 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. Firstly, principal component analysis (PCA) was used to check the structure of the data. This method, widely used in provenance studies, reduces the number of relevant variables by defining new variables (or components) as a linear combination of the original ones (i.e., the compositional variables). The two components that accumulate the largest percentages of variance (PC1 and PC2) were represented in a biplot to investigate if any clustering appeared and to check if the clusters correlate with any of the reference fields. As usual, data were standardized to zero mean and unit variance before computing the analysis. PCA is categorized as an unsupervised method, meaning that the procedure does not consider the class labels. In contrast, supervised machine learning methods use labeled reference datasets to train models to predict the class of unlabeled data (i.e., the archaeological samples). Training is usually performed using 80% of the reference data, and then the performance of the trained model is tested using the remaining 20%. The split between train and test subsets is made by using a random seed. This procedure is well-known in machine learning, but to our knowledge, it has never been applied to provenance analyses of grinding stones. All the compositional values were divided by the SiO₂ concentration as a way to standardize them. Different 8:2 splits were tested (using different random seeds) following the train–test approach. The supervised approach was performed using the “Supervised Provenance Analysis” R Markdown files from Anglisano et al. [43], 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 [44]. 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 to glomerophyric and vesiculated texture with rather idiomorphic phenocrysts of two compositions embedded within a microlitic groundmass (Figure 3a). The most abundant and largest phenocrysts correspond to clinopyroxene (cpx), with augite composition as analyzed by SEM-EDS, see Figure 4 and Table 1. Cpx often appears as clusters of polygonal idiomorphic crystals, some of them exhibiting occasional hourglass sector zoning (Figure 3b).

Smaller phenocrysts of olivine (ol) are also very common in the three samples. These are anhedral (often rounded and occasionally corroded, showing that they partly dissolved into the liquid now solidified as the groundmass). The ol phenocrysts systematically exhibit an iddingsite alteration rim, made of iron oxides and clay minerals, resulting from post-magmatic hydrothermal alteration (typical for olivine in many volcanic rocks) (Figure 3c,d). The fresh olivine composition measured by SEM-EDS is Fo_75±2, while the obtained composition in the altered rims is still essentially olivine but with larger Fe contents (19–29 wt% instead of 18 wt%). Possibly, the measured spectra, even if focused on the rims, have contributions from iddingsite and nearby olivine. 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 mostly unoriented lath-shaped plagioclase microlites, clinopyroxene and opaque equidimensional minerals (possibly ilmenite and Ti-bearing magnetite).

All the collected samples from Sidi Zahruni exhibit a porphyritic and vesiculated texture with anhedral to idiomorphic phenocrysts and microphenocrysts within a microlitic groundmass (Figure 5). The mineral association is uniform and consists of the following: (i) Zoned and twinned plagioclase (pl) crystals with rectangular and quadratic sections (Figure 5a). SEM-EDS analyses indicate that they have compositions ranging from 46 to 78 An% (labradorite and bytownite). (ii) Clinopyroxene (cpx) with a slight brown color (augite as analyzed by SEM-EDS; see Table 1 and Figure 4). These augite crystals appear particularly idiomorphic, exhibiting mainly hexagonal sections or octagonal basal sections (Figure 5e), often showing the diagnostic two cleavage directions at nearly 90° angles. Hourglass sector zoning also appears in some crystals (Figure 5f). (iii) Olivine (ol) with rare polygonal sections and more common rounded shapes (Figure 5g), sometimes significantly corroded. Its composition by SEM-EDS is Fo_72±4. Apart from these, in one of the samples (Z6), a large green-colored (PPL) rounded crystal was spotted (Figure 5h). The SEM-EDS analysis revealed Na- and Fe-rich clinopyroxene, a member of the aegirine–augite solid solution (see Table 1). The different phenocrysts are also embedded in a fine-grained matrix mainly formed by lath-shaped plagioclase crystals, clinopyroxene and opaque equidimensional minerals (possibly ilmenite or Ti-bearing magnetite).

The treated images from the scanned petrographic thin sections (e.g., Figure 6) allowed the discrimination between porosity, matrix and the different types of phenocrysts and their quantification. Porosity is around 24%–39% for Iulia Libica (L) samples and usually lower (6%–30%) for Sidi Zahruni (Z) samples, although, for sample Z5, it 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 similar to some of the samples from Sidi Zahruni. However, samples Z3 and Z4 have narrower pore size distributions. The circularity of the pores is within the range of 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 and Z7).

The mineral part of the samples is mainly formed by the fine-grained microlitic matrix in all the samples, and the phenocrysts are only a minor component (3%–12%), see Table 2. The matrix/phenocrysts content ratio is higher for the Sidi Zahruni samples (constituting 20 ± 8) compared to the Iulia Libica samples (8 ± 3).

The two types of phenocrysts that appear in Iulia Libica samples have a relative content ratio of about 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 a mineral, the frequency of sizes increases towards small values forming a continuous range of sizes that connect the microphenocrysts (area < 0.25 mm²) to the microlites of the matrix (Figure 8). The area of the phenocrysts only rarely reaches values above 1 mm².

The distribution of clinopyroxene and olivine crystals is less biased towards small crystals for the basalts 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 plagioclase 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 for Sidi Zahruni samples. However, plagioclase phenocrysts from these basalts deviate markedly from circularity consistently with their common lath-shaped morphologies (see <C> in Figure 8).

The petrographic characterization indicates that all the investigated samples (L and Z) are volcanic rocks exhibiting basaltic compositions and are classifiable as such according to the international classification of igneous rocks [46].

3.2. Geochemistry

The average bulk composition of the analyzed basalts is shown in Table 3. Each set of samples has 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 all of the investigated samples can be classified as basalts (Figure 9) according to the TAS diagram [42]. The higher concentration of Al, Ca and Na in Z samples is the consequence of a magma richer in these elements and where plagioclase crystallized earlier, forming the pl phenocrysts only observed in the basalts from Sidi Zahruni. In contrast, the amount of Fe is slightly higher in the samples from Iulia Libica (~14% compared to ~12%). Finally, the trace element proportions are also dissimilar for both types of basalts. Those retrieved in Iulia Libica are richer in trace elements, except for the amount of V, which is higher in the basalts from Sidi Zahruni.

The tectonic setting of volcanic lavas can sometimes be geochemically inferred using certain compositional ratios (e.g., Zr vs. TiO₂ or Zr vs. Zr/Y) [47,48]. According to these criteria, the measured geochemical values indicate that the samples (both from Iulia Libica and Sidi Zahruni) correspond to intraplate basalts, which could help determine the provenance of the samples. However, the classification of volcanic districts as intraplate can be challenging for complex tectonic settings like those of some Mediterranean volcanic districts.

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 from Sidi Zahruni appears to be attributable to the Pantelleria island (strait of Sicily) and those from Iulia Libica to the Hyblean Plateau (southeastern Sicily). In contrast, for samples not lying within any reference field, it would be assumed that they do not belong to any of the considered provenances or, alternatively, that the reference fields perhaps are not well defined due to a lack of a fully representative set of reference samples. For instance, one of the samples from Sidi Zahruni strictly is not 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 will be unambiguously assigned to a given provenance. Looking at Figure 10b, the basalts 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 those from the Middle Atlas (Morocco). Also, using different biplots, divergent 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 plotted 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 condenses the relevant information in a few principal components that can also be represented in biplots. There are a few precedents of the use of this statistical tool in stone provenance studies (see, e.g., [30,51,52]). However, one major difficulty is to have a uniform set of data. The reference data from the GEOROC database is heterogeneous since samples were prepared using different methods and measured using different equipment [53]. The analyzed elements are not always the same set; in particular, the minor and trace elements can vary.

Two PCAs were computed using two different sets of data. For the PCA shown in Figure 11, a total of 1702 records (G90 dataset) from the preselected GEOROC database were combined with the 25 individual analyses from Sidi Zahruni samples and the 8 from Iulia Libica. To obtain this homogeneous dataset, the number of considered compositional variables was restricted to the most represented within the preselected database (SiO₂, TiO₂, Al₂O₃, CaO, MgO, Mn, 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 us 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 SiO₂ and Al₂O₃ and a negative correlation between them with Mn, P₂O₅ and TiO₂. 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₂O and CaO, and both are negatively correlated with Na₂O. Both Sidi Zahruni and Iulia Libica samples 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, it is not possible to conclude on 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 the G90 dataset. Reference samples that are geographically close to Iulia Libica (L), such as those labeled Catalonia, Calatrava and Massif Central, are plotted near the L samples, but individual samples from many other reference sites are plotted in the same area. Instead of assigning provenances, it 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 no 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. Most of the discarded provenances agree with the fact that both Z and L basalts are intraplate. From the discarded volcanic districts, only Algeria is considered intraplate.

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 as strong as in the previous PCA. Moreover, 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 better 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 either the Z or L clusters (as expected since these volcanic arcs are not in an intraplate tectonic setting), 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 plot relatively close to each other, some provenances match the position of both clusters. This is the case for Catalonia, Egypt, Atlas and Sicily Channel Rift clusters (all of them considered intraplate volcanic districts). Additionally, the samples from the Suez Rift province match particularly well with those from Sidi Zahruni. Moreover, 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 volcanic districts are underrepresented (e.g., <6 samples) within the G30 database (see Table S4), and they could be non-statistically representative of their location, and some potential provenance locations have actually 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 [43,44] was applied to the reference data obtained from the GEOROC database. Three different databases were tested (G90, G60 and G30, as described in Section 2.2). Datasets G30 and G60 only differ in two features (Cu and Zn values are not considered in G60; compare Tables S3 and S4). From the different classes (i.e., provenances) defined for the data in these three databases, only those bearing at least nine samples were considered. Classes with a smaller number of samples cannot 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 [54,55] were added to the database to keep the corresponding provenance into consideration. Indeed, this class (X8), labeled Catalonia, corresponds to the Olot volcanic field, only about 50 km southeast of Iulia Libica. An additional reason to retain this class is that the corresponding provenance was, in fact, suggested from the results of some of the computed PCAs. 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 point to a given provenance for the Z samples. In contrast, for the 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) come from the Sicily Channel Rift (i.e., Pantelleria island), which is not far from the archaeological site (~115 km). Using the G90 database and combining all the Z samples and all the trained models, the overall probability for class X21 (Sicily Channel Rift) reaches 88.4%, and only two more classes are suggested, Anatolia (X5) with 11.1% and Egypt (X10) with only 0.5%. Using the G60 database, only two provenances are suggested: X21 captures an overall probability of 96%, and X5 captures the remaining 4%. Finally, using G30, class X21 holds almost all the probability (96.4%) again, and the rest splits between Anatolia (X5) with 3.3% and Suez (X22) withonly 0.3%.

Contrastingly, the different models do not show a similar agreement to assign provenance to the samples from Iulia Libica (L), and there are even significant differences using different seeds to define the train/test splits (see, for instance, L-sample predictions using the ANN model in Table 4). In any case, using the 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 a bit more than 50%. However, after checking the results obtained using G60 (Table 5), we see that the predominance of class X5 is lost. The overall results using the 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 spread again among different provenances: 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 accurate quantification of a set of parameters and a dataset containing the same quantification for a representative collection of reference geological samples.

In this paper, a quantification effort has been undertaken for both 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 combine 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 mentioning that a set of volcanic samples that are mineralogically and petrographically different could exhibit chemical homogeneity.

In the present investigation, the samples 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 plagioclase (absent as phenocryst in L samples), clinopyroxenes and olivine. The absence of plagioclase phenocrysts (as 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 of olivine (L samples) or the occasional occurrence of green aegirine–augite (Z samples). Alteration rims in olivine are actually a common feature in many olivine basalts worldwide and have many potential provenances, including Sardinia, Agde (to the south of Massif Central), Olot (Catalonia), Ustica (Italy) and Middle Atlas (Morocco) basalts [56]. The presence of aegirine–augite crystals has been described in volcanic rocks from Pantelleria and is typical of more felsic rocks from this island [57], with which the basalts are intimately associated [58].

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 who publish their data, and some initiatives have been undertaken to systematically collect these data to create huge geochemical reference databases like GEOROC. However, some difficulties arise from the use of different equipment and sample preparation methods by the contributors to the reference database. Nevertheless, part of the solution is to build larger and larger datasets, which contributes to reducing random errors (improving precision), and the average of different systematic errors from different laboratories could also increase accuracy. Also, an extra complication of large datasets is that 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.). All these points could imply additional transformation steps before using the data, but they should not represent a real issue for well-managed databases.

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 the 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 can also be used to discard some of the multiple considered provenances. Provenances that can be discarded both for the samples from Iulia Libica and Sidi Zahruni include the basalts from the Aeolian, Aegean, Gibraltar and Tyrrhenian arcs and those from Campania, Lazio and possibly Algeria, Alps–Sardinia–Corsica and Etna. The fact that the suggested provenances correspond to those of intraplate basalts (intraplate nature was also geochemically deduced for L and Z samples) and that the discarded reference sites are essentially those linked to volcanic arcs (i.e., linked to a subduction zone) attest the robustness of both GEOROC database and the PCA method. However, PCA is not really useful in deciphering which, among the suggested intraplate provenances, is the most probable.

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 between different classes can be measured using different statistical metrics computed using the test sets. Along this line, as was previously stated in a methodologically similar study [43], using supervised models is very important to check for the convergence between the results using different classification models and training several times using different splits. Reliable provenance determination requires a significant agreement between the results using these different strategies. In the present research, besides different models and splits, three different reference databases have been used (G90, G60 and G30) with accuracies ranging from 0.55–0.80 (using G90) to 0.66–0.90 (using G30).

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

In contrast, for the samples from Iulia Libica, there is no agreement between the prediction models, and the predicted provenance usually appears statistically divided into 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 [43], this heterogeneity indicates that the provenance attribution is not robust. Indeed, in some cases, even the results from a given classification algorithm vary significantly when using different trained models (i.e., different train–test sets). Some classes (among others, incidentally X5, Anatolia) are overrepresented within the reference databases and may spread out, covering most of the point cloud in the corresponding PCA (e.g., yellow dots in Figure 11 and Figure 12). Therefore, it is not surprising that such classes could capture substantial fractions of the provenance probability. Using geographical proximity as a criterion, classes X17 (Massif Central) and especially X8 (Catalonia) are the nearest locations to Iulia Libica. The models assign to these classes rather low probabilities, and both classes have a very low amount of data within the reference databases. Even if training is enabled, a low number of data may not be fully statistically representative of the class. In fact, class X17 has only been considered by the models using database G90 because the number of X17-class data in the other datasets 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 [36], typologically similar to those found in Iulia Libica. Such a trading route would be consistent with the presence of pottery from the southern Gaulish area in Iulia Libica [32]. With the available data and comparing with the results obtained for Sidi Zahruni, the conclusion is that the provenance of Iulia Libica samples has not been successfully determined, and this implies that it corresponds to a class imperfectly represented or absent from the reference databases used.

To summarize, 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 limitations linked to incomplete reference databases (unbalanced, underrepresented or missing classes). Despite the good indicators of overall performance (see accuracies in Table 4, Table 5 and Table 6) 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. One of the statistical metrics used to check how well a model performs for a particular class is sensitivity, which is defined as the ratio between well-predicted samples and the total samples within the class. 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 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 the LSVM model, and this is actually the only model that assigns a certain probability to this provenance.

5. Conclusions

Statistical approaches to establish the 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. Implementing 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 and statistical classification tools have been illustrated. On the one hand, biplots using a 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 for discerning different reference clusters, and this enables provenance prediction for unlabeled samples. Among the 18 different 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, not only because of the strong agreement between the different supervised models but also because of 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 km from the archaeological site of Sidi Zahruni.

In contrast, the models assign different provenances to the samples from Iulia Libica, and this suggests that the corresponding provenance was missing in the considered reference datasets or that it was incompletely described within them. This highlights the limitations of the supervised approach and the need for a statistical approach making use of different sister samples, different trained models and different classification algorithms to be able to identify a homogeneous and, therefore, robust provenance prediction.

Author Contributions

Conceptualization, L.C.; fieldwork and sample retrieval, L.C., B.F. and C.C.; experimental work—petrographic analyses, L.C. and R.D.F.; experimental work—electron microcopy, R.D.F.; experimental work—geochemical analyses, Á.P.M., A.A., R.D.F. and I.Q.; code writing, A.A. and L.C.; formal analyses of data, L.C., I.Q. and Á.P.M.; writing—first draft preparation, L.C.; writing—review and editing, all authors. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Most of the data presented in this study have been made publicly available. Any other data are available on request from the corresponding author.

Acknowledgments

We are grateful to the Tunisian Institut National du Patrimoine for allowing access and sampling at the Sidi Zahruni site. We thank Néjia Laridhi-Ouazaa, Siwar Baklouti, Aureli Álvarez, Marta Prevosti and Ramon Járrega for their support during fieldwork in Tunisia. Toni Lozano and Albert Ruiz are acknowledged for giving training on Python and machine learning. A.P.M. is a Serra Húnter Lecturer. L.C. and A.P.M. are members of the Grup de Recerca Conservació-Restauració de Patrimoni (AGAUR, 2021 SGR 00089). Finally, we would like to thank the editor as well as three anonymous reviewers for their insightful remarks and comments that definitely have helped to improve the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. (a) Location of the archaeological sites of Iulia Libica (Spain) and Sidi Zahruni (Tunisia). (b) The three sampled millstones from Iulia Libica. (c) The seven small volcanic samples retrieved from Sidi Zahruni.

Figure 2. Map of the Mediterranean area, including the location of all the selected data from the geochemical database GEOROC, colored or labeled according to the corresponding data group. The archaeological sites where the studied basalts were retrieved have also been highlighted.

Figure 3. Petrographic images of the samples from Iulia Libica illustrating their volcanic origin and basaltic composition. (a) PPL (left) and XPL (right) micrograph of the same area, showing abundant clinopyroxene (cpx) and olivine (ol) phenocrysts within a microlitic, vesiculated groundmass. (b) XPL micrograph of a cpx aggregate including hourglass sector zoning crystals. (c) XPL micrograph of a subhedral, corroded ol crystal, including an iddingsite alteration rim. (d) XPL micrograph including a basal section of cpx showing typical two cleavage directions at near 90° (left) and a corroded and altered ol crystal (right).

Figure 4. Pyroxene ternary diagram [45] showing the compositions of the pyroxene phenocrysts in samples from Iulia Libica (blue) and Sidi Zahruni (red).

Figure 5. Petrographic images of the samples from Sidi Zahruni illustrating their volcanic origin and basaltic composition. (a–d) General aspect of representative samples showing plagioclase (pl) phenocrysts (a) and microphenocrysts (a–d), vesicles (mostly on c,d), olivine (ol) and clinopyroxene (cpx) microphenocrysts (mostly on a,d) within a microlithic groundmass viewed in XPL. (e) Basal section of a cpx microphenocryst, PPL. (f) cpx hourglass sector zoning microphenocryst, XPL. (g) Detail of an ol microphenocryst, XPL. (h) Section of a sodic, partly corroded, cpx macrocryst showing a brighter reaction rim with the groundmass, PPL.

Figure 6. (a) L3, Iulia Libica, and (b) Z3, Sidi Zahruni, thin section images. To the left are the original XPL scans, and to the right are the corresponding colored images, including porosity (red), clinopyroxene (green), olivine (blue), plagioclase (black) and the microcrystalline matrix (white).

Figure 7. Porosity size histograms from the samples from Iulia Libia (top) and Sidi Zahruni (bottom). Next to each histogram, an inset shows the treated image with pores in black.

Figure 8. Crystal size (area) histograms for phenocrysts and microphenocrysts in basalts from (a) Iulia Libica (L) and (b) Sidi Zahruni (Z): the quantified minerals are clinopyroxene and olivine for L and Z, and plagioclase for Z samples. <C> indicates the average circularity values of the phenocrysts.

Figure 9. TAS diagram according to the IUGS classification scheme [42] showing the data plots for Iulia Libica (blue) and Sidi Zahruni (red) sample analyses.

Figure 10. Average compositions of samples from Iulia Libica (black dots) and Sidi Zahruni (open dots) plotted in different discrimination diagrams: (a) Sr vs. TiO2; (b) Nb vs. TiO2; (c) Zr/Nb vs. TiO2. Reference compositional fields from [28,49,50].

Figure 11. PCA biplot of factor scores for the first two principal components for the selected data from GEOROC (G90 dataset), Iulia Libica and Sidi Zahruni sample analyses have also been plotted. Inset: PCA biplot of the original variables.

Figure 12. PCA biplot of factor scores for the first two principal components for the selected data from GEOROC (G30 dataset), Iulia Libica and Sidi Zahruni sample analyses have also been plotted. Inset: PCA biplot of the original variables.

Table 1

Average structural formulae of the analyzed (SEM-EDS) clinopyroxenes calculated on the basis of six oxygens. Eight augite crystals were measured for L samples and Z samples, as well as a single additional augite–aegirine crystal, which was found in sample Z6.

	Iulia Libica (L) Samples	Sidi Zahruni (Z) Samples
	Augite	Augite	Augite–Aegirine
Si	1.85 ± 0.05	1.87 ± 0.05	2.03
Ti	0.04 ± 0.01	0.04 ± 0.01	0.01
Al	0.17 ± 0.03	0.16 ± 0.06	0.02
Fe	0.33 ± 0.03	0.38 ± 0.05	0.55
Mn	n.d.	n.d.	0.03
Mg	0.78 ± 0.04	0.82 ± 0.05	0.57
Ca	0.76 ± 0.01	0.74 ± 0.07	0.75
Na	0.03 ± 0.01	0.01 ± 0.01	0.12

Table 2

Quantification of the different components identified in samples from Iulia Libica (L) and Sidi Zahruni (Z). Average values are also shown. Percentages were measured as area %, but they can also be read as vol%.

Sample	Phenocrysts			Matrix	Porosity
	cpx (%)	ol (%)	pl (%)	(%)	(%)
L1	7.5	4.6	0.0	64.3	23.6
L2	3.4	1.9	0.0	56.1	38.6
L3	5.8	2.4	0.0	66.3	25.5
L	6 ± 2	3 ± 2	-	62 ± 6	29 ± 8
Z1	1.3	2.1	3.5	86.9	6.2
Z2	0.5	1.4	4.1	64.1	29.9
Z3	1.1	1.1	3.6	81.0	13.2
Z4	0.5	0.7	1.4	78.9	18.5
Z5	0.5	0.5	2.2	46.0	50.8
Z6	0.7	0.8	3.0	76.5	19.0
Z7	0.4	0.7	1.5	74.0	23.4
Z	0.7 ± 0.4	1.0 ± 0.5	3 ± 1	73 ± 14	23 ± 14

Table 3

Average composition (m) and corresponding standard deviations (σ) from N specimens from L and Z samples of the investigated archaeological basalts, including major (wt%) and trace elements (ppm).

	L1		L2		L3		Z1		Z2		Z4		Z5		Z7
N	3		2		3		5		5		5		5		5
Statistics	m	σ	m	σ	m	σ	m	σ	m	σ	m	σ	m	σ	m	σ
SiO2 (%)	47.63	0.25	48.01	0.11	47.04	0.21	47.59	0.08	48.06	0.10	47.50	0.17	46.43	0.22	46.76	0.24
Al2O3 (%)	14.73	0.10	13.51	0.24	13.54	0.09	16.61	0.11	16.13	0.08	16.02	0.04	15.58	0.12	15.21	0.17
Fe2O3 (%)	14.73	0.10	13.51	0.24	13.54	0.09	11.82	0.03	11.89	0.03	12.89	0.04	12.30	0.12	12.32	0.06
MgO (%)	5.61	0.17	5.60	0.40	6.23	0.10	5.59	0.02	6.11	0.05	5.10	0.05	4.99	0.03	4.79	0.08
CaO (%)	9.48	0.08	9.74	0.21	10.24	0.21	11.89	0.02	11.57	0.06	10.90	0.03	11.79	0.06	10.81	0.03
Na2O (%)	2.80	0.14	2.74	0.08	2.30	0.04	3.47	0.04	3.28	0.03	3.80	0.05	3.39	0.03	3.42	0.03
K2O (%)	1.77	0.03	1.77	0.06	1.64	0.03	1.02	0.01	1.01	0.01	1.06	0.01	0.92	0.01	1.01	0.01
TiO2 (%)	2.37	0.02	2.27	0.05	2.31	0.04	2.63	0.01	2.55	0.01	3.01	0.02	2.69	0.05	2.83	0.01
P2O5 (%)	0.76	0.01	0.85	0.01	0.72	0.01	0.51	0.01	0.51	0.00	0.60	0.00	0.55	0.01	0.61	0.01
Mn (ppm)	1376	14	1415	47	1686	89	1324	9	1352	9	1440	5	1332	20	1385	8
V (ppm)	268	11	266	18	277	18	333	17	316	19	379	22	364	23	358	22
Cr (ppm)	201	1	221	13	337	27	136	4	128	8	82	6	95	3	84	5
Ni (ppm)	181	2	204	9	280	11	83	0	88	0	57	0	61	0	56	0
Cu (ppm)	60	2	68	3	59	5	33	2	84	2	38	3	54	1	65	2
Zn (ppm)	123	2	126	7	133	9	85	1	96	1	97	1	91	2	96	1
Rb (ppm)	43	2	42	3	39	1	14	1	13	1	15	1	14	1	13	1
Sr (ppm)	898	7	938	60	849	41	496	1	596	1	471	1	500	2	476	2
Y (ppm)	29	0	31	2	31	1	22	0	22	1	27	0	23	1	27	1
Zr (ppm)	239	1	252	15	227	12	135	1	136	1	148	1	151	1	150	1
Nb (ppm)	79	3	80	8	71	3	30	3	31	1	32	4	25	3	30	3

Table 4

Accuracies computed using the predictions of the trained models on three different test subsets (from dataset G90) and corresponding average probabilities of predicted class for Sidi Zahruni and Iulia Libica samples.

Model	Seed	Accuracy	Class 1 Predictions
			Sidi Zahruni Samples	Iulia Libica Samples
GLM	1	0.6049	X21 (80%) X5 (20%)	X5 (62.5%) X1 (25%) X15 (12.5%)
	2	0.6292	X21 (68%) X5 (32%)	X5 (50%) X1 (37.5%) X15 (12.5%)
	3	0.5957	X21 (60%) X5 (40%)	X5 (50%) X1 (25%) X6 (12.5%) X15 (12.5%)
RF	1	0.7629	X21 (100%)	X5 (75%) X6 (25%)
	2	0.7508	X21 (100%)	X5 (62.5%) X6 (37.5%)
	3	0.7447	X21 (100%)	X6 (50%) X5 (37.5) X15 (12.5%)
ANN	1	0.6292	X21 (100%)	X15 (50%) X5 (37.5%) X1 (12.5%)
	2	0.5805	X21 (96%) X5 (4%)	X13 (75%) X21 (25%)
	3	0.5745	X21 (56%) X5 (44%)	X1 (62.5%) X13 (25%) X5 (12.5%)
kkNN	1	0.7660	X21 (100%)	X5 (75%) X8 (12.5%) X21 (12.5%)
	2	0.7264	X21 (100%)	X5 (37.5%) X8 (37.5%) X15 (12.5%) X21 (12.5%)
	3	0.7508	X21 (100%)	X5 (62.5%) X8 (37.5%)
LDA	1	0.5653	X21 (60%) X5 (40%)	X5 (87.5%) X15 (12.5%)
	2	0.5593	X21 (52%) X5 (48%)	X5 (62.5%) X1 (25%) X15 (12.5%)
	3	0.5471	X21 (52%) X5 (48%)	X5 (75%) X15 (12.5%) X21 (12.5%)
Stack	1	0.7538	X21 (100%)	X5 (100%)
	2	0.7295	X21 (100%)	X5 (87.5%) X15 (12.5%)
	3	0.7204	X21 (100%)	X5 (87.5%) X6 (12.5%)
GB	1	0.6770	X21 (92%) X5 (8%)	X5 (62.5%) X15 (37.5%)
	2	0.6794	X21 (88%) X10 (12%)	X13 (50%) X15 (37.5%) X5 (12.5%)
	3	0.6914	X21 (80%) X5 (16%) X10 (4%)	X13 (62.5%) X8 (25%) X5 (12.5%)
GPC	1	0.6733	X21 (100%)	X5 (50%) X15 (50%)
	2	0.6962	X21 (92%) X5 (8%)	X5 (50%) X15 (37.5%) X21 (12.5%)
	3	0.6962	X21 (100%)	X5 (50%) X15 (50%)
GNB	1	0.8254	X21 (100%)	X5 (50%) X15 (50%)
	2	0.8014	X21 (100%)	X5 (50%) X15 (50%)
	3	0.8014	X21 (100%)	X5 (50%) X15 (50%)
LSVM	1	0.5670	X21 (80%) X5 (20%)	X1 (37.5%) X5 (37.5%) X8 (25%)
	2	0.5861	X21 (96%) X5 (4%)	X17 (87.5%) X1 (12.5%)
	3	0.5981	X21 (100%)	X8 (37.5%) X5 (25%) X17 (25%) X6 (12.5%)

¹ Class key: X1: Adria Domain; X5: Anatolia; X6: Atlas; X8: Catalonia; X10: Egypt; X13: Hyblean Plateau; X15: Lybia; X17: Massif Central; X21: Sicily Channel Rift.

Table 5

Accuracies computed using the predictions of the trained models on three different test subsets (from dataset G60) and corresponding average probabilities of predicted class for Sidi Zahruni and Iulia Libica samples.

Model	Seed	Accuracy	Class 1 Predictions
			Sidi Zahruni Samples	Iulia Libica Samples
GLM	1	0.7282	X21 (88%) X5 (12%)	X6 (62.5%) X13 (37.5%)
	2	0.7538	X21 (80%) X5 (20%)	X6 (87.5%) X13 (12.5%)
	3	0.7795	X21 (88%) X5 (12%)	X6 (50%) X13 (50%)
RF	1	0.8359	X21 (100%)	X6 (87.5%) X13 (12.5%)
	2	0.8256	X21 (100%)	X6 (100%)
	3	0.8718	X21 (100%)	X6 (100%)
ANN	1	0.6974	X21 (100%)	X6 (62.5%) X13 (37.5%)
	2	0.7538	X21 (100%)	X6 (100%)
	3	0.7590	X21 (100%)	X6 (100%)
kkNN	1	0.8205	X21 (96%) X5 (4%)	X5 (50%) X8 (37.5%) X13 (12.5%)
	2	0.8154	X21 (100%)	X5 (50%) X8 (37.5%) X13 (12.5%)
	3	0.8718	X21 (96%) X5 (4%)	X5 (50%) X8 (37.5%) X21 (12.5%)
LDA	1	0.6564	X21 (100%)	X13 (100%)
	2	0.7231	X21 (100%)	X13 (100%)
	3	0.7231	X21 (100%)	X13 (100%)
Stack	1	0.8359	X21 (96%) X5 (4%)	X5 (50%) X6 (37.5%) X13 (12.5%)
	2	0.8205	X21 (100%)	X5 (50%) X6 (37.5%) X13 (12.5%)
	3	0.8718	X21 (96%) X5 (4%)	X5 (50%) X6 (37.5%) X21 (12.5%)
GB	1	0.8135	X21 (80%) X5 (20%)	X13 (62.5%) X6 (37.5%)
	2	0.8016	X21 (80%) X5 (20%)	X13 (100%)
	3	0.7619	X21 (80%) X5 (20%)	X13 (75%) X6 (25%)
GPC	1	0.7500	X21 (100%)	X21 (50%) X5 (37.5%) X15 (12.5%)
	2	0.7381	X21 (100%)	X21 (62.5%) X5 (25%) X10 (12.5%)
	3	0.7302	X21 (100%)	X5 (37.5%) X10 (37.5%) X21 (25%)
GNB	1	0.8413	X21 (100%)	X5 (37.5%) X10 (37.5%) X21 (25%)
	2	0.8690	X21 (100%)	X5 (37.5%) X10 (37.5%) X21 (25%)
	3	0.8651	X21 (100%)	X5 (37.5%) X10 (37.5%) X21 (25%)
LSVM	1	0.6825	X21 (100%)	X6 (87.5%) X1 (12.5%)
	2	0.6786	X21 (100%)	X6 (50%) X8 (37.5%) X13 (12.5%)
	3	0.6666	X21 (100%)	X6 (50%) X8 (37.5%) X13 (12.5%)

¹ Class key: X1: Adria Domain; X5: Anatolia; X6: Atlas; X8: Catalonia; X10: Egypt; X13: Hyblean Plateau; X15: Lybia; X21: Sicily Channel Rift.

Table 6

Accuracies computed using the predictions of the trained models on three different test subsets (from dataset G30) and corresponding average probabilities of predicted class for Sidi Zahruni and Iulia Libica samples.

Model	Seed	Accuracy	Class 1 Predictions
			Sidi Zahruni Samples	Iulia Libica Samples
GLM	1	0.8416	X21 (100%)	X13 (75%) X8 (25%)
	2	0.8614	X21 (100%)	X15 (50%) X8 (37.5%) X13 (12.5%)
	3	0.8614	X21 (96%) X5 (4%)	X15 (87.5%) X8 (12.5%)
RF	1	0.8713	X21 (100%)	X6 (62.5%) X5 (12.5%) X13 (12.5%) X21 (12.5%)
	2	0.8416	X21 (100%)	X6 (87.5%) X21 (12.5%)
	3	0.8713	X21 (100%)	X6 (100%)
ANN	1	0.8119	X21 (100%)	X6 (75%) X5 (12.5%) X15 (12.5%)
	2	0.7525	X21 (100%)	X6 (100%)
	3	0.8317	X21 (100%)	X8 (50%) X6 (37.5%) X5 (12.5%)
kkNN	1	0.8911	X21 (92%) X5 (8%)	X5 (37.5%) X6 (25%) X8 (12.5%) X13 (12.5%) X21 (12.5%)
	2	0.8812	X21 (92%) X5 (8%)	X6 (50%) X5 (12.5%) X8 (12.5%) X13 (12.5%) X21 (12.5%)
	3	0.9010	X21 (92%) X5 (8%)	X5 (62.5%) X6 (25%) X13 (12.5%)
LDA	1	0.7921	X21 (100%)	X13 (100%)
	2	0.7723	X21 (96%) X22 (4%)	X13 (100%)
	3	0.8317	X21 (100%)	X13 (100%)
Stack	1	0.8713	X21 (100%)	X5 (37.5%) X6 (37.5%) X13 (12.5%) X21 (12.5%)
	2	0.8416	X21 (92%) X5 (8%)	X6 (50%) X5 (12.5%) X8 (12.5%) X13 (12.5%) X21 (12.5%)
	3	0.8812	X21 (92%) X5 (8%)	X5 (62.5%) X6 (25%) X13 (12.5%)
GB	1	0.8271	X21 (80%) X5 (20%)	X6 (50%) X5 (25%) X10 (12.5%) X15 (12.5%)
	2	0.8272	X21 (80%) X5 (20%)	X6 (50%) X5 (25%) X10 (12.5%) X15 (12.5%)
	3	0.7894	X21 (80%) X5 (12%) X22 (8%)	X15 (62.5%) X5 (37.5%)
GPC	1	0.7594	X21 (100%)	X5 (37.5%) X15 (37.5%) X21 (25%)
	2	0.7368	X21 (100%)	X15 (62.5%) X5 (37.5%)
	3	0.7820	X21 (100%)	X15 (62.5%) X5 (37.5%)
GNB	1	0.8571	X21 (100%)	X15 (62.5%) X5 (37.5%)
	2	0.8947	X21 (100%)	X15 (62.5%) X5 (37.5%)
	3	0.8571	X21 (100%)	X15 (62.5%) X5 (37.5%)
LSVM	1	0.7669	X21 (100%)	X6 (50%) X8 (37.5%) X13 (12.5%)
	2	0.7444	X21 (100%)	X6 (62.5%) X8 (37.5%)
	3	0.6617	X21 (100%)	X6 (37.5%) X8 (37.5%) X13 (25%)

¹ Class key: X1: Adria Domain; X5: Anatolia; X6: Atlas; X8: Catalonia; X10: Egypt; X15: Lybia; X21: Sicily Channel Rift; X22: Suez.

Supplementary Materials

The following supporting information can be downloaded at www.mdpi.com/article/10.3390/min14070639/s1. Table S1: Database contains 1855 records obtained by selecting analyzed basalt rocks within the area contained within latitude [25° N–54° N] and longitude [12° W to 45° E] from GEOROC database [41]. Table S2: Dataset G90 contains all records from Table S1, which includes the following set of analyzed elements: Si, Ti, Al, Ca, Mg, Mn, K, Na and P. Table S3: Dataset G60 contains all records from Table S2, which also includes the following set of analyzed elements: V, Ni, Rb, Sr, Y, Zr and Nb. Table S4: Dataset G30 contains all records from Table S3, which includes additional analyses on the Cu and Zn contents. The Python code for this study is attached: [https://www.mdpi.com/article/10.3390/min14070639/s1].

References

1. Liritzis, I.; Laskaris, N.; Vafiadou, A.; Karapanagiotis, I.; Volonakis, P.; Papageorgopoulou, C.; Bratitsi, M. Archaeometry: An Overview. Sci. Cult.; 2020; 6, pp. 49-98.

2. Kolb, C.C. Provenance Studies in Archaeology. Encyclopedia of Global Archaeology; Smith, C. Springer: New York, NY, USA, 2014; pp. 6172-6181. ISBN 978-1-4419-0465-2

3. Waksman, Y. Provenance Studies: Productions and Compositional Groups. The Oxford Handbook of Archaeological Ceramic Analysis; Alice, H. Oxford University Press: Oxford, UK, 2017; pp. 148-161.

4. Pollard, A.M.; Heron, C.; Gillard, R.D. The Geochemistry of Clays and the Provenance of Ceramics. Archaeological Chemistry; Pollard, A.M.; Heron, C.; Gillard, R.D. The Royal Society of Chemistry: Cambridge, UK, 2008; pp. 98-143. ISBN 978-0-85404-262-3

5. Lancaster, L.; Sottili, G.; Marra, F.; Ventura, G. Provenancing of Lightweight Volcanic Stones Used in Ancient Roman Concrete Vaulting: Evidence from Rome. Archaeometry; 2011; 53, pp. 707-727. [DOI: https://dx.doi.org/10.1111/j.1475-4754.2010.00565.x]

6. Maniatis, Y.; Herz, N.; Basiakos, Y. The Study of Marble and Other Stones Used in Antiquity: Asmosia III Athens: Transactions of the 3rd International Symposium of the Association for the Study of Marble and Other Stones Used in Antiquity; Archetype Publications Ltd.: London, UK, 1995.

7. Flügel, E.; Flügel, C. Applied Microfacies Analysis: Provenance Studies of Roman Mosaic Stones. Facies; 1997; 37, pp. 1-48. [DOI: https://dx.doi.org/10.1007/BF02537369]

8. García de Madinabeitia, S.; Gil Ibarguchi, J.I.; Santos Zalduegui, J.F. IBERLID: A Lead Isotope Database and Tool for Metal Provenance and Ore Deposits Research. Ore Geol. Rev.; 2021; 137, 104279. [DOI: https://dx.doi.org/10.1016/j.oregeorev.2021.104279]

9. Levine, M.A. Determining the Provenance of Native Copper Artifacts from Northeastern North America: Evidence from Instrumental Neutron Activation Analysis. J. Archaeol. Sci.; 2007; 34, pp. 572-587. [DOI: https://dx.doi.org/10.1016/j.jas.2006.06.015]

10. Falkenberg, J.; Kaplan, U.; Mutterlose, J. Understanding the Provenance and Production Process of Historic Mortars—A Novel Approach Employing Calcareous Nannofossils. Archaeol. Anthr. Sci.; 2023; 15, 138. [DOI: https://dx.doi.org/10.1007/s12520-023-01840-2]

11. Miriello, D.; Barca, D.; Pecci, A.; De Luca, R.; Crisci, G.M.; López Luján, L.; Barba, L. Plasters from Different Buildings of the Sacred Precinct of Tenochtitlan (Mexico City): Characterization and Provenance. Archaeometry; 2015; 57, pp. 100-127. [DOI: https://dx.doi.org/10.1111/arcm.12074]

12. Gratuze, B. Provenance Analysis of Glass Artefacts. Modern Methods for Analysing Archaeological and Historical Glass; John Wiley & Sons: Hoboken, NJ, USA, 2013; pp. 311-343. ISBN 9781118314234

13. Di Febo, R.; Casas, L.; Capelli, C.; Cabella, R.; Vallcorba, O.; Di Febo, R.; Casas, L.; Capelli, C.; Cabella, R.; Vallcorba, O. Catalan Imitations of the Ligurian Taches Noires Ware in Barcelona (18th–19th Century): An Example of Technical Knowledge Transfer. Minerals; 2018; 8, 183. [DOI: https://dx.doi.org/10.3390/min8050183]

14. Fornacelli, C.; Briano, A.; Chiarantini, L.; Bianchi, G.; Benvenuti, M.; Giamello, M.; Kang, J.S.; Villa, I.M.; Talarico, F.M.; Hodges, R. Archaeometric Provenance Constraints for Early Medieval Sparse Glazed Pottery from Donoratico (Livorno, Italy). Archaeometry; 2021; 63, pp. 549-576. [DOI: https://dx.doi.org/10.1111/arcm.12633]

15. Velliky, E.C.; MacDonald, B.L.; Porr, M.; Conard, N.J. First Large-Scale Provenance Study of Pigments Reveals New Complex Behavioural Patterns during the Upper Palaeolithic of South-Western Germany. Archaeometry; 2021; 63, pp. 173-193. [DOI: https://dx.doi.org/10.1111/arcm.12611]

16. Doménech-Carbó, A.; Doménech-Carbó, M.T.; Valle-Algarra, F.M.; Gimeno-Adelantado, J.V.; Osete-Cortina, L.; Bosch-Reig, F. On-Line Database of Voltammetric Data of Immobilized Particles for Identifying Pigments and Minerals in Archaeometry, Conservation and Restoration (ELCHER Database). Anal. Chim. Acta; 2016; 927, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.aca.2016.04.052]

17. Ouazaa, N.L.; Casas, L.; Álvarez, A.; Fouzai, B.; Moreno-Vide, M.; Vidal, L.; Sihem, R.; Sonzogni, C.; Borschneck, D. Provenance of Marble Sculptures from the National Museum of Carthage (Tunisia). J. Archaeol. Sci.; 2013; 40, pp. 1602-1610. [DOI: https://dx.doi.org/10.1016/j.jas.2012.10.035]

18. Casas, L.; Di Febo, R.; Parcerisa, D. Petrographic Markers for Archaeometric Identification of Montjuïc Sandstone, the Flagship Stone of Barcelona (NE Spain). Minerals; 2020; 10, 154. [DOI: https://dx.doi.org/10.3390/min10020154]

19. Bigazzi, G.; Ercan, T.; Oddone, M.; Özdoḡan, M.; Yeḡingil, Z. Application of Fission Track Dating to Archaeometry: Provenance Studies of Prehistoric Obsidian Artifacts. Nucl. Tracks Radiat. Meas.; 1993; 22, pp. 757-762. [DOI: https://dx.doi.org/10.1016/0969-8078(93)90172-Z]

20. Re, A.; Giudice, A.L.; Angelici, D.; Calusi, S.; Giuntini, L.; Massi, M.; Pratesi, G. Lapis Lazuli Provenance Study by Means of Micro-PIXE. Nucl. Instrum. Methods Phys. Res. B; 2011; 269, pp. 2373-2377. [DOI: https://dx.doi.org/10.1016/j.nimb.2011.02.070]

21. Malfilatre, C.; Hallot, E.; Boulvais, P.; Poujol, M.; Chauvin, A.; Gapais, D.; Dabard, M.-P.; Bourquin, S.; Pallix, D. Fingerprinting the Provenance of Building Stones: A Case Study on the Louvigné and Lanhélin Granitic Rocks (Armorican Massif, France). Bull. Société Géologique Fr.; 2014; 185, pp. 13-31. [DOI: https://dx.doi.org/10.2113/gssgfbull.185.1.13]

22. Malfilatre, C.; Boulvais, P.; Dabard, M.-P.; Bourquin, S.; Hallot, E.; Pallix, D.; Gapais, D. Petrographical and Geochemical Characterization of Comblanchien Limestone (Bourgogne, France): A Fingerprint of the Building Stone Provenance. Comptes Rendus Geosci.; 2012; 344, pp. 14-24. [DOI: https://dx.doi.org/10.1016/j.crte.2011.12.002]

23. Williams-Thorpe, O.; Thorpe, R.S. Geochemistry and Trade of Eastern Mediterranean Millstones from the Neolithic to Roman Periods. J. Archaeol. Sci.; 1993; 20, pp. 263-320. [DOI: https://dx.doi.org/10.1006/jasc.1993.1018]

24. Peacock, D.P.S. The Roman Millstone Trade: A Petrological Sketch. World Archaeol.; 1980; 12, pp. 43-53. [DOI: https://dx.doi.org/10.1080/00438243.1980.9979779]

25. Santi, P.; Foresta Martin, F.; Spatafora, F.; de Vita, S.; Renzulli, A. Volcanic Grinding Tools in Ustica Island (Tyrrhenian Sea, Italy): Local Production vs. Import of Morgantina-Type Millstones in the Hellenistic-Roman Period. Minerals; 2020; 10, 389. [DOI: https://dx.doi.org/10.3390/min10050389]

26. Lancaster, L.C.; Sottili, G.; Marra, F.; Ventura, G. Provenancing of Lightweight Volcanic Stones Used in Ancient Roman Concrete Vaulting: Evidence from Turkey and Tunisia. Archaeometry; 2010; 52, pp. 949-961. [DOI: https://dx.doi.org/10.1111/j.1475-4754.2009.00509.x]

27. Jackson, M.; Marra, F. Roman Stone Masonry: Volcanic Foundations of the Ancient City. Am. J. Archaeol.; 2006; 110, pp. 403-436. [DOI: https://dx.doi.org/10.3764/aja.110.3.403]

28. Antonelli, F.; Lazzarini, L. Mediterranean Trade of the Most Widespread Roman Volcanic Millstones from Italy and Petrochemical Markers of Their Raw Materials. J. Archaeol. Sci.; 2010; 37, pp. 2081-2092. [DOI: https://dx.doi.org/10.1016/j.jas.2010.02.008]

29. Di Bella, M.; Italiano, F.; Martinelli, M.C.; Mazzoleni, P.; Quartieri, S.; Tigano, G.; Tripodo, A.; Sabatino, G. Archeometric Characterization of Prehistoric Grindstones from Milazzo Bronze Age Settlement (Sicily, Italy). Archaeol. Anthr. Sci.; 2018; 10, pp. 1571-1583. [DOI: https://dx.doi.org/10.1007/s12520-017-0483-8]

30. Falcone, F.; Dionisio, A.; Castorina, F.; Tufo, A.; Francis, R.E.; Stoppa, F. Archaeometry of a Roman Millstone from Santa Maria Arabona, Manoppello (Abruzzo, Central Italy). Minerals; 2021; 11, 948. [DOI: https://dx.doi.org/10.3390/min11090948]

31. Santi, P.; Renzulli, A.; Gullo, R. Archaeometric Study of the Hopper-Rubber and Rotary Morgantina-Type Volcanic Millstones of the Greek and Roman Periods Found in the Aeolian Archipelago (Southern Italy). Eur. J. Mineral.; 2013; 25, pp. 39-52. [DOI: https://dx.doi.org/10.1127/0935-1221/2013/0025-2250]

32. Guàrdia, J.; Carreras, C.; de Soto, P. Commercial Circuits of Iulia Lybica (Llívia): Notions Gleaned from Finds of Pottery and Marble. Rev. D’arqueologia Ponent; 2017; pp. 147-168. [DOI: https://dx.doi.org/10.21001/rap.2017.27.6]

33. Járrega Domínguez, R.; Prevosti Monclús, M. La Producción de Ánforas En El Territorium de Neapolis Durante La Antigüedad Tardía El Taller de Sidi Zahruni (Nabeul, Túnez); BAR Publishing: Oxford, UK, 2022.

34. Guàrdia Felip, J. The City and the Forum of Iulia Libica (Llívia, Cerdanya). Cuad. Arqueol. Univ. Navar.; 2021; 29, pp. 111-128.

35. Carreras, C.; Farré, E. Un Model de Rellotge Solar Trobat a Iulia Libica (Llívia). Rev. D’arqueologia Ponent; 2021; 31, pp. 115-128.

36. Py, M. Meules d’époque Protohistorique et Romaine Provenant de Lattes. Lattara, Recherches sur l’Économie Vivrière des Lattarenses; Py, M. Association Pour La Recherche Archéologique En Languedoc Oriental: Lattes, France, 1992; Volume 5, pp. 183-232.

37. Baklouti, S.; Maritan, L.; Casas, L.; Laridhi Ouazaa, N.; Jàrrega, R.; Prevosti, M.; Mazzoli, C.; Fouzaï, B.; Larabi Kassaa, S.; Fantar, M. Establishing a New Reference Group of Keay 25.2 Amphorae from Sidi Zahruni (Nabeul, Tunisia). Appl. Clay Sci.; 2016; 132–133, pp. 140-154. [DOI: https://dx.doi.org/10.1016/j.clay.2016.05.027]

38. Prevosti, M.; Járrega, R.; Casas, L.; Fouzai, B.; Baklouti, S.; Laridhi-Ouazaa, N. Sidi Zahruni, Una Alfarería En El Ager de Neapolis (Túnez): Aglomeración o Vicus circa Villam?. Dialogues D’histoire Ancienne; 2022; 481, pp. 265-313. [DOI: https://dx.doi.org/10.3917/dha.481.0265]

39. Schneider, C.A.; Rasband, W.S.; Eliceiri, K.W. NIH Image to ImageJ: 25 Years of Image Analysis. Nat. Methods; 2012; 9, pp. 671-675. [DOI: https://dx.doi.org/10.1038/nmeth.2089]

40. Marguí, E.; Queralt, I. Sample Preparation for X-ray Fluorescence Analysis. Encyclopedia of Analytical Chemistry; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2024; pp. 1-29. ISBN 9780470027318

41. Lehnert, K.; Su, Y.; Langmuir, C.H.; Sarbas, B.; Nohl, U. A Global Geochemical Database Structure for Rocks. Geochem. Geophys. Geosystems; 2000; 1, 1012. [DOI: https://dx.doi.org/10.1029/1999GC000026]

42. Le Bas, M.J.; Le Maitre, R.W.; Streickeisen, A.; Zanettin, B. A Chemical Classification of Volcanic Rocks Based on the Total Alkali-Silica Diagram. J. Petrol.; 1986; 27, pp. 745-750. [DOI: https://dx.doi.org/10.1093/petrology/27.3.745]

43. Anglisano, A.; Casas, L.; Queralt, I.; Di Febo, R. Supervised Machine Learning Algorithms to Predict Provenance of Archaeological Pottery Fragments. Sustainability; 2022; 14, 11214. [DOI: https://dx.doi.org/10.3390/su141811214]

44. Anglisano, A.; Casas, L.; Anglisano, M.; Queralt, I. Application of Supervised Machine-Learning Methods for Attesting Provenance in Catalan Traditional Pottery Industry. Minerals; 2020; 10, 8. [DOI: https://dx.doi.org/10.3390/min10010008]

45. Morimoto, N.; Fabries, J.; Ferguson, A.K.; Ginzburg, I.V.; Ross, M.; Seifert, F.A.; Zussman, J.; Aoki, K.; Gottardi, G. Nomenclature of Pyroxenes. Miner. Mag.; 1988; 52, pp. 535-550. [DOI: https://dx.doi.org/10.1180/minmag.1988.052.367.15]

46. Le Maitre, R.W.; Streckeisen, A.; Zanettin, B.; Le Bas, M.J.; Bonin, B.; Bateman, P. Igneous Rocks: A Classification and Glossary of Terms: Recommendations of the International Union of Geological Sciences Subcommission on the Systematics of Igneous Rocks; 2nd ed. Cambridge University Press: Cambridge, UK, 2002; ISBN 9780521619486

47. Pearce, J.A. Trace Element Characteristics of Lavas from Destructive Plate Boundaries; Thorpe, R.S. John Wiley and Sons: Chichester, UK, 1982.

48. Pearce, J.A.; Norry, M.J. Petrogenetic Implications of Ti, Zr, Y, and Nb Variations in Volcanic Rocks. Contrib. Mineral. Petrol.; 1979; 69, pp. 33-47. [DOI: https://dx.doi.org/10.1007/BF00375192]

49. Weller, D.; Orlandini, O.; Benton, J.; Schirmer, C.; LoBue, L.; Culotta, S. Provenancing the Stone Tools of Volubilis, Morocco: A Socio-Economic Interpretation of Stonework Lithologies. J. Archaeol. Sci. Rep.; 2021; 38, 103105. [DOI: https://dx.doi.org/10.1016/j.jasrep.2021.103105]

50. Antonelli, F.; Lazzarini, L.; Luni, M. Preliminary Study on the Import of Lavic Millstones in Tripolitania and Cyrenaica (Libya). J. Cult. Herit.; 2005; 6, pp. 137-145. [DOI: https://dx.doi.org/10.1016/j.culher.2004.10.005]

51. Orange, M.; Le Bourdonnec, F.-X.; Bellot-Gurlet, L.; Lugliè, C.; Dubernet, S.; Bressy-Leandri, C.; Scheffers, A.; Joannes-Boyau, R. On Sourcing Obsidian Assemblages from the Mediterranean Area: Analytical Strategies for Their Exhaustive Geochemical Characterisation. J. Archaeol. Sci. Rep.; 2017; 12, pp. 834-844. [DOI: https://dx.doi.org/10.1016/j.jasrep.2016.06.002]

52. Germinario, L.; Hanchar, J.M.; Sassi, R.; Maritan, L.; Cossio, R.; Borghi, A.; Mazzoli, C. New Petrographic and Geochemical Tracers for Recognizing the Provenance Quarry of Trachyte of the Euganean Hills, Northeastern Italy. Geoarchaeology; 2018; 33, pp. 430-452. [DOI: https://dx.doi.org/10.1002/gea.21666]

53. Gluhak, T.M.; Schwall, C. Provenance Analyses of the Volcanic Rock Grinding Stones from the Greek Colony of Selinunte, Sicily (Italy)—Constraints and Possibilities. Archaeometry; 2015; 57, pp. 246-268. [DOI: https://dx.doi.org/10.1111/arcm.12086]

54. Gisbert Pinto, G. Estudi Petrològic i Geoquímic Del Vulcanisme Recent de La Garrotxa. Master’s Thesis; Universitat de Barcelona: Barcelona, Spain, 2009.

55. Enrique, P.; Toribio, V. Subalkali Basalt Flows from the Ser Basin (Olot Volcanic Field). Geogaceta; 2009; 47, pp. 129-132.

56. Williams-Thorpe, O. Provenancing and Archaeology of Roman Millstones from the Mediterranean Area. J. Archaeol. Sci.; 1988; 15, pp. 253-305. [DOI: https://dx.doi.org/10.1016/0305-4403(88)90066-0]

57. White, J.C.; Parker, D.F.; Ren, M. The Origin of Trachyte and Pantellerite from Pantelleria, Italy: Insights from Major Element, Trace Element, and Thermodynamic Modelling. J. Volcanol. Geotherm. Res.; 2009; 179, pp. 33-55. [DOI: https://dx.doi.org/10.1016/j.jvolgeores.2008.10.007]

58. Ferla, P.; Meli, C. Evidence of Magma Mixing in the ‘Daly Gap’ of Alkaline Suites: A Case Study from the Enclaves of Pantelleria (Italy). J. Petrol.; 2006; 47, pp. 1467-1507. [DOI: https://dx.doi.org/10.1093/petrology/egl015]

59. Renzulli, A.; Santi, P.; Gambin, T.; Serrano, P.B. Pantelleria Island as a Centre of Production for the Archaic Phoenician Trade in Basaltic Millstones: New Evidence Recovered and Sampled from a Shipwreck off Gozo (Malta) and a Terrestrial Site at Cádiz (Spain). J. Archaeol. Sci. Rep.; 2019; 24, pp. 338-349. [DOI: https://dx.doi.org/10.1016/j.jasrep.2019.01.019]

60. Civetta, L.; D‘Antonio, M.; Orsi, G.; Tilton, G.R. The Geochemistry of Volcanic Rocks from Pantelleria Island, Sicily Channel: Petrogenesis and Characteristics of the Mantle Source Region. J. Petrol.; 1998; 39, pp. 1453-1491. [DOI: https://dx.doi.org/10.1093/petroj/39.8.1453]

61. Gail, A.M.; Don, R.B. Experimental Constraints on Depths of Fractionation of Mildly Alkalic Basalts and Associated Felsic Rocks; Pantelleria, Strait of Sicily. Contrib. Mineral. Petrol.; 1986; 93, pp. 251-264. [DOI: https://dx.doi.org/10.1007/bf00371327]