# A miner's perspective on Pb isotope provenances in the Western and Central Mediterranean

Francis Albarède, Janne Blichert-Toft, Liesel Gentelli, Jean Milot, Markos Vaxevanopoulos, Sabine Klein, Katrin Westner, Thomas Birch, Gil Davis, François de Callataÿ

► To cite this version:

Francis Albarède, Janne Blichert-Toft, Liesel Gentelli, Jean Milot, Markos Vaxevanopoulos, et al. A miner's perspective on Pb isotope provenances in the Western and Central Mediterranean. Journal of Archaeological Science, 2020, 121, pp.105194. 10.1016/j.jas.2020.105194. hal-02903203

HAL Id: hal-02903203

https://hal.science/hal-02903203

Submitted on 4 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A miner's perspective on Pb isotope provenances in the Western and Central Mediterranean

Francis Albarede<sup>a, *,</sup>, Janne Blichtert-Toft<sup>a</sup>, Liesel Gentelli<sup>a</sup>, Jean Milot<sup>a</sup>, Markos Vaxevanopoulos<sup>a</sup>, Sabine Klein<sup>b</sup>, Katrin Westner<sup>b</sup>, Thomas Birch<sup>c</sup>, Gil Davis<sup>d</sup>, François de Callatayj<sup>e, f</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Pb isotopes

Silver

Coins

Ores

Galena

## ABSTRACT

The conventional approach to ore provenance studies of ancient silver coins and artifacts has been to first analyze and then try to match them to published data about mining districts, a difficult task given our incomplete knowledge of these. While literary sources are useful to identify possible provenances, they potentially bias interpretations proper because of a variety of limitations of their time. Archeological evidence in the form of mining shafts, galleries, spoil heaps, and tools also provides a tangible and reliable record of mining, but dating such activity can be problematic and the record is inconsistent. Here we propose a new approach driven by Pb isotopic data rather than numismatic groups. Statistical analysis of Pb isotopic data is used to identify ore-defined isotopic clusters. This new method is based on an algorithm that predicts the number of isotopic clusters necessary to fulfill the simple condition that variance within isotopic clusters is minimized whereas inter-cluster variance is maximized. Since each cluster reflects a discrete geological episode within a particular environment broadly datable to a specific Pb model age, it can be identified as a potential source exploited by ancient miners. We explore the potential of this method in two examples using data from coins and ores respectively. In the first example, Roman Republican silver coins form three 'end-member' clusters sourced in mining districts with Cenozoic, Mesozoic, and Paleozoic Pb model ages. The example demonstrates how sources of silver used to mint coinage of the Roman Republic shifted within 50 years of the end of the Second Punic War in 201 BCE. In the second example using Aegean galena samples, Pb isotopes distinguish components with model ages datable to the Hercynian basement, the recent Aegean tectonic province, and Cyprus, noting that significant silver mining districts may remain unidentified in either Spain or the Aegean world. We further clarify a number of potential analytical issues and advocate that users of Pb isotopes for tracing archeological artifacts measure all four lead isotopes and inspect the 12 proposed isotope combinations in order to select those that provide the best representation of the data. We also emphasize that full advantage should be taken of the geologically informed parameters (model age and Th/U/Pb relationships) to identify the geological context of metal sources.

## 1. Introduction

Powerful techniques exist to assess the origin of metals used to produce artifacts such as votives, silverware, jewelry, and coinage which, together, represented a significant store of wealth in ancient economies. Among these techniques, lead isotope analysis (LIA) has

become prevalent for silver artifacts since this metal was most commonly extracted by cupellation of lead-bearing charges (galena, cerussite, anglesite, jarosite). The conventional approach to provenance appraisal is to determine the lead isotope composition of a sample taken from a particular object and compare it with the closest values from published Pb isotope data on ores and relevant metallurgical residues (i.

<footnote><sup>a</sup> Ecole Normale Supérieure de Lyon, CNRS, and Université de Lyon, France</footnote>

<footnote><sup>b</sup> Deutsches Bergbau-Museum Bochum, Germany</footnote>

<footnote><sup>c</sup> Centre for Urban Network Evolutions (UrbNet), Aarhus University, Hejbjerg, Denmark</footnote>

<footnote><sup>d</sup> Department of Ancient History, Macquarie University, NSW, Australia</footnote>

<footnote><sup>e</sup> Royal Library of Belgium, Brussels, Belgium</footnote>

<footnote><sup>f</sup> Ecole Pratique des Hautes Etudes, Paris, France</footnote>

<footnote>* Corresponding author.
E-mail address: albarede@ens-lyon.fr (F. Albarede).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** The four possible ternary plots involving ratios with the same reference (denominator) isotope showing how Pb isotope ratios evolve through time as a function of the tectonic, or geological, 'model' age $T_m$ and the $\mu$ ($^{238}\text{U}/^{204}\text{Pb} \propto \text{U/Pb}$) and $\kappa$ ($^{232}\text{Th}/^{238}\text{U} \propto \text{Th/U}$) parameters. Panels (a) to (d) correspond to normalization to $^{204}\text{Pb}$, $^{206}\text{Pb}$, $^{207}\text{Pb}$, and $^{208}\text{Pb}$, respectively. The blue and red stem-like segments correspond to evolution of the ratios with radiogenic ingrowth through geological time from 540 Ma ago to 0 (modern lead, represented as full circle symbols). The numbers labelling the lines indicate $\mu$ values ($^{238}\text{U}/^{204}\text{Pb} = 8$, 10, 12). The blue and red colors indicate different $\kappa$ values ($^{232}\text{Th}/^{238}\text{U} = 3.8$ and 4.2, respectively). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

e., slags or litharge, cf. Blichert-Toft et al., 2016; Stos-Gale and Gale, 2009). This method is adequate to determine whether connecting a particular artifact with a particular ore or archeological site is consistent with isotopic observations, but does not explore the full topology of the database. In addition, the effect of analytical uncertainties, especially on data produced prior to the 2000s, harms the reliability of isotopic connections. A recent review by Artioli et al. (2020) on Pb isotopes around the Mediterranean describes the databases and clarifies some of the issues related to different elements such as copper, lead, silver, and tin used to manufacture coinage and a variety of other objects.

Attempts to further the isotopic relationships between artifacts and potential metal sources have been made using Euclidian and Mahalanobis distances to identify all the ores representing a statistically acceptable provenance for a particular object (Birch et al., 2020a, 2020b; De Ceuster and Degryse, 2020; Delile et al., 2015; Westner et al., 2020). Yet, this approach also does not provide a holistic or synoptic overview of Pb isotope compositions of artifacts and ore populations. In other words, neither artifacts nor ores are considered parts of isotopic populations and, as they are mostly visible for coins, leaves numismatic identification as the undisputed tool for unraveling metal fluxes.

Here, we propose a complementary, data-driven approach to conventional provenance studies. Lead isotope data on ores may be considered as forming ore-defined clusters because geological processes acting during particular geological events and within particular environments lend themselves to a treatment as populations of events rather than isolated accidents (Albarede et al., 2012; Blichert-Toft et al., 2016). Silver-bearing mineralizations are, for example, typically associated with Cenozoic events in the Aegean-Anatolian region (Melfos and

Voudouris, 2017; Yigit, 2009) and in the Betic chains (Arribas and Tosdal, 1994; García-Hernández et al., 1980), whereas Pan-African (late Proterozoic and early Phanerozoic) domains of Africa are rich in gold but, with exceptions such as the Moroccan Anti-Atlas ore deposits (Baron et al., 2020; Cheilletz et al., 2002), devoid of major silver ores.

We describe a method to identify clusters in Pb isotope space taken from databases independent of arbitrary assignments, then test the method in two studies, one of coins (numismatics) and one of ores (geology). The nature and the goal of this approach are different from modern classification techniques derived from artificial intelligence, such as the multiple classifier approach of elemental abundances developed by Żabiński et al. (2020) and based on major and trace elements in slags. The four lead isotopes only provide a three-dimensional description of variability, but they are physically and geologically consistent and, because isotopes do not fractionate at high temperatures, are relatively robust with respect to metallurgical processes. In addition, the selection of the isotopic ratios best used to define the isotopic space will be rationalized with respect to environmental parameters (geological ages, U/Pb, and Th/U) with a description of their sensitivity to analytical issues.

Major and trace elements on the one hand and Pb isotope variations on the other follow very different rules because the underlying physics behind their variations are fundamentally different. In spite of multiple attempts to combine isotopic and major and trace element markers, a given major or trace element pattern remains intrinsically difficult to relate to a class of ores, a locality, or a particular metallurgical workshop, especially for potentially cupelled noble metals. Consequently, we have not attempted to consider the systematics of major and trace
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** The twelve binary plots corresponding to the orthogonal projections of the lines shown in Fig. 1 onto each face. For symbol and color codes refer to Fig. 1. The three plots showing red color lines only do not involve the radiogenic ingrowth of <sup>208</sup>Pb from Th. The ellipses centered on the points representing modern Pb show how uncertainties on mass bias correction affect results (here 0.002 per difference of a.m.u.). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

element data in this particular context.

## 2. Discriminant analysis

### 2.1. Background

#### 2.1.1. Lead has four isotopes

* <sup>206</sup>Pb (~25%) and <sup>207</sup>Pb (~21%), which are still being produced by the radioactive decay of <sup>238</sup>U and <sup>235</sup>U, respectively

* <sup>208</sup>Pb (~53%), which is still being produced by the radioactive decay of <sup>232</sup>Th

* the small (~1.3%) <sup>204</sup>Pb isotope, which is stable on the time scales of geological processes and hence not radiogenic.

The rate of radiogenic ingrowth decreases in the order <sup>207</sup>Pb, <sup>206</sup>Pb, and <sup>208</sup>Pb. Because <sup>235</sup>U has largely decayed away over the aeons, the modern ingrowth of <sup>207</sup>Pb is extremely slow and can for all practical intents and purposes be considered to be zero. The four isotopes can be represented in 3-dimensional (3D) space as a coordinate triplet (Fig. 1) which itself can be projected onto the respective 2-dimensional (2D) faces (Fig. 2). The linearity required for clean statistical tests requires that all ratios have the same isotope as the denominator, contrary to standard usage in numismatic and archeological literature. The choice of isotopes combined into ratios will be discussed further below.

#### 2.1.2. Principle

We use hierarchical cluster analysis, an algorithm which collects similar objects into clusters using objective distance criteria. The

algorithm starts by treating each observation as a separate cluster. It then identifies the two observations that are the closest and merges them into a cluster. This is repeated until a minimum number of clusters that satisfies a chosen statistical criterion is obtained. The algorithm identifies the 'tree' of successive clustering steps and can be stopped when either a prescribed number of clusters has been reached or a particular statistic, usually associated with a probability level, has reached a specified value. Different algorithms using different metrics ('distances') between clusters are available from the literature (Kettenring, 2006). One way of defining cluster analysis is that intra-cluster distances are minimized, whereas inter-cluster distances are maximized (Szekely and Rizzo, 2005; Székely and Rizzo, 2013).

#### 2.1.3. Data decorrelation (whitening)

Given the strong correlations among Pb isotope ratios involving <sup>204</sup>Pb (Albarède et al., 2004), we begin by removing the overall correlation structures by first centering the data around the mean and then rotating the data in 3-dimensional space until the correlation coefficients between variables are zero. Each variable is then stretched or contracted in such a way that standard deviation is equal to unity in each direction. Whitening requires eigen component analysis of the centered data.

#### 2.1.4. Selecting a representation

Any of the four Pb isotopes can be arbitrarily picked as the reference isotope (the denominator). Each 'raw Pb isotope' 3D space has a unique mapping into the 'geologically informed' 3D parameter space, the model age $T_m$ (usually the age of the last local major tectonic event), $\mu = \text{<sup>238</sup>U/\text{<sup>204</sup>Pb}$ (or short-hand, U/Pb) of the underlying crust, and $\kappa = \text{<sup>232</sup>Th/}$
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Cluster analysis applied to Pb isotope compositions of silver *denarii* and *victoriati* from the Roman Republic (Albarède et al., 2016; Desaulty et al., 2011; Westner et al., 2020) in three dimensions and assuming a maximum number (n) of 12 clusters.</figcaption></figure>

²³⁸U (or, short-hand, Th/U) (Albarede et al., 2012; see also Appendix A). Fig. 1 shows the four possible 3D plots corresponding to the four possible reference isotopes, while Fig. 2 shows a column-organized representation of the four possible projections on each face of Fig. 1, these being the twelve 2D plots achievable in LIA.

Each set of plots has its strengths and weaknesses in provenance studies:

* ²⁰⁴Pb normalization (Fig. 1a and 2, col. 1): although in theory, T_m, μ, and κ are cleanly separated by these plots, the abundance of this isotope is smaller by factors of 15–40 relative to the abundances of the other isotopes, which introduces correlated noise on the three ratios (Albarède et al., 2004). Steep linear arrays are commonly observed in both the y = ²⁰⁷Pb/²⁰⁴Pb vs x = ²⁰⁶Pb/²⁰⁴Pb plot (see, for example, some fields of literature data used as background in recent publications: Artioli et al. (2020); Eshel et al. (2019). Similar linear arrays also appear in z = ²⁰⁸Pb/²⁰⁴Pb vs x = ²⁰⁶Pb/²⁰⁴Pb plots. These arrays are not isochrons (see Appendix A) but rather reflect mass-dependent fractionation effects due to (i) weathering or

low-temperature diagenesis under unusually oxidizing conditions (Pb²⁺-Pb⁴⁺) (Fujii et al., 2011; Yang and Liu, 2015), (ii) incomplete yields during Pb purification chemistry, or (iii) incorrect correction of mass-spectrometer bias. Mass fractionation lines with slopes indistinguishable from the predicted theoretical slopes are visible in some pre-2000 literature data, particularly data on ores used as background in recent publications (Artioli et al., 2020; Birch et al., 2020a; Eshel et al., 2019) (Appendix A). These steep mass fractionation trends should not be confused with the shallow trends also visible in these papers for newly acquired data; in this case the trends often reproduce the Pb isotope growth curve, which, in the plots discussed above, is much less steep than the slope of mass fractionation curves. Such shallow trends reflect that radiogenic Pb and initial Pb occupy different mineralogical sites and are easily separated by weathering and laboratory leaching procedures (Frei and Kamber, 1995; Harlavan et al., 2009). Altogether, mass bias and mass fractionation lines are much less frequently observed in more recent data sets. Modern multiple-collector inductively-coupled plasma mass spectrometry (MC-ICP-MS) seems to have resolved the
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Same as Fig. 3 (n = 12) but as binary projections on selected faces of the ternary plots. (x) Volcanic rocks from the southern coast of Andalusia and the Alboran Ridge (Duggen et al., 2004); (+) galenas and slags from La Unión, Murcia (Baron et al., 2017). The bar shows the uncertainty associated with ±1‰ per a.m.u fractionation during either low-temperature natural processes or mass spectrometric analysis. See Supplementary Fig. S1 for a comprehensive plot of all possible pairs of ratios.</figcaption></figure>

analytical issues associated with instrumental mass bias owing to the addition of Tl as mass bias monitor. Double- and tripe-spike techniques (Galer, 1999) also solve both yield and instrumental mass bias issues, but double the analytical burden (Albarède et al., 2004) and have not so far received much attention in archeometry.

* $^{206}$Pb normalization (Fig. 1b and 2, col. 2): this family of plots has achieved preeminence in geochronology, notably for meteorites and lunar rocks (Tera and Wasserburg, 1972). It presents only weak correlation effects since the ratios of the large isotopes $^{207}$Pb/$^{206}$Pb and $^{208}$Pb/$^{206}$Pb vary within a small range with most of the variation being due to $^{204}$Pb/$^{206}$Pb, which is the most prone to recent geological effects and analytical mass bias. As a consequence of the narrow range of $^{207}$Pb/$^{204}$Pb for recent tectonic ages due to the

near-exhaustion of $^{235}$U, which for all practical intents and purposes can be considered extinct today, the $^{207}$Pb/$^{206}$Pb vs $^{204}$Pb/$^{206}$Pb plot is of little practical interest. It only shows that the actual number of degrees of freedom in Pb isotope space is two (the product µT_m and κ) rather than three. Because of the slow decay of $^{238}$U and $^{232}$Th, $^{208}$Pb/$^{206}$Pb is only barely affected by the tectonic age T_m and the $^{208}$Pb/$^{206}$Pb vs $^{204}$Pb/$^{206}$Pb, or even better, $^{207}$Pb/$^{206}$Pb, plots offer a direct depiction of Th/U variations in ore sources.

* $^{207}$Pb normalization (Fig. 1c and 2, col. 3) (Erel et al., 1994): the $^{206}$Pb/$^{207}$Pb and $^{208}$Pb/$^{207}$Pb vs $^{204}$Pb/$^{207}$Pb plots offer a powerful representation of µ variations. The abscissa does not vary much (Getty and Depaolo, 1995). The $^{208}$Pb/$^{207}$Pb vs $^{206}$Pb/$^{207}$Pb offers a $^{204}$Pb-independent assessment of κ values (Th/U ratios).
<caption>**Table 1**</caption>

Clusters obtained on 99 Pb isotope data on *denarii* and *victoriati* minted during the Roman Republic. Probability matrix showing whether two clusters may be samples from a same population. When *p* > 0.05 (italic), clusters cannot be considered significantly different.

<table><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.003</td><td>0.001</td></tr><tr><th>2</th><td>0.001</td><td>1</td><td>0.001</td><td>0.001</td><td>0.021</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.008</td><td>0.006</td></tr><tr><th>3</th><td>0.001</td><td>0.001</td><td>1</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.007</td></tr><tr><th>4</th><td>0.001</td><td>0.001</td><td>0.001</td><td>1</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.016</td><td>0.002</td><td>0.01</td><td>0.014</td></tr><tr><th>5</th><td>0.001</td><td>0.021</td><td>0.001</td><td>0.001</td><td>1</td><td>0.001</td><td>0.002</td><td>0.001</td><td>0.002</td><td>0.006</td><td>0.02</td><td>0.026</td></tr><tr><th>6</th><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>1</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.06</td><td>0.207</td><td>0.009</td></tr><tr><th>7</th><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.002</td><td>0.001</td><td>1</td><td>0.003</td><td>0.002</td><td>0.013</td><td>0.043</td><td>0.018</td></tr><tr><th>8</th><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.003</td><td>1</td><td>0.002</td><td>0.023</td><td>0.413</td><td>0.033</td></tr><tr><th>9</th><td>0.001</td><td>0.001</td><td>0.001</td><td>0.016</td><td>0.002</td><td>0.001</td><td>0.002</td><td>0.002</td><td>1</td><td>0.3</td><td>0.014</td><td>0.036</td></tr><tr><th>10</th><td>0.001</td><td>0.001</td><td>0.001</td><td>0.002</td><td>0.006</td><td>0.06</td><td>0.013</td><td>0.023</td><td>0.3</td><td>1</td><td>0.185</td><td>0.07</td></tr><tr><th>11</th><td>0.003</td><td>0.008</td><td>0.001</td><td>0.01</td><td>0.02</td><td>0.207</td><td>0.043</td><td>0.413</td><td>0.014</td><td>0.185</td><td>1</td><td>0.235</td></tr><tr><th>12</th><td>0.001</td><td>0.006</td><td>0.007</td><td>0.014</td><td>0.026</td><td>0.009</td><td>0.018</td><td>0.033</td><td>0.036</td><td>0.07</td><td>0.235</td><td>1</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Parameters *T*<sub>m</sub>, μ, and *k* for the 12-cluster solution. Clusters 4, 9, and 12 give Paleozoic model ages (Hercynian basement). A supergroup corresponding to young *T*<sub>m</sub> (Cenozoic ages) consists of clusters 1 and 6 and corresponds to the recent magmatic phenomenon associated with the Betic orogeny. Clusters 2, 3, 5, and 7 reveal a Mesozoic end-member unrelated to a well-identified orogenic event and likely correspond to a large-scale hydrothermal occurrence known throughout Western Europe.</figcaption></figure>

* $^{208}$Pb normalization (Fig. 1d and 2, col. 4): we show in an example below that the various representations based on this isotope are less noisy and less prone to mass biases. They offer smooth visual discrimination of μ and κ values in the $^{207}$Pb/$^{208}$Pb vs $^{206}$Pb/$^{208}$Pb plot. In contrast, *T*<sub>m</sub> is not easily isolated from any $^{208}$Pb-normalized plot.

We address cases below for which the $^{208}$Pb/$^{206}$Pb vs $^{207}$Pb/$^{206}$Pb plot and the $^{208}$Pb normalization provide strong constraints on μ variations almost independently of variations in *T*<sub>m</sub> and κ, but this may differ from one data set to another. What appears to be valid for all data sets is that plots involving $^{204}$Pb are the noisiest.

### 2.1.5. Implementation

The cluster algorithm used here is the commercial MATLAB™ package but equivalent packages exist in *R* and *Python*. We chose Euclidian distances described in the Ward algorithm and used by Birch

et al. (2020a, 2020b); Delile et al. (2014); Westner et al. (2020). The data were whitened before processing. The clusters were reordered by decreasing population. In general, increasing the number of clusters above a certain value only increases the number of outliers without modifying the main clusters. The quality of cluster identification is evaluated using the non-parametric test based on statistical energy (Székely and Rizzo, 2013; https://github.com/brian-lau/multdist). When too many data overlap, the cluster relationships can be made visible by contouring 2D histograms.

We also considered data sets in their geographic context. In order to smooth regional variability, and to some extent minimize analytical issues, the data were averaged and gridded over longitude-latitude windows and a maximum number of clusters was allowed. For better readability, the relative frequencies calculated for each cluster at the center of each longitude-latitude window were plotted as pie charts using a Unix executable code and GMT map software version 5.3 (Wessel et al., 2013). The cluster numbers are eventually reshuffled so that the
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Same as Fig. 4 but with *n* = 4. Such a small number clearly under-represents the variability and may be too small to produce truly compact groups. Note the alignments indicative of mixing. See Supplementary Fig. S2 for a comprehensive plot of all possible pairs of ratios.</figcaption></figure>

higher the cluster number, the smaller the number of the member samples. Low-rank cluster numbers therefore indicate well-populated clusters.

**3. Results and discussion**

We explored the potential of the new method described above in two examples using data sets of Roman Republican silver coins and Aegean ores.

**3.1. Roman Republican silver coinage**

We first tested the clustering method on a non-geolocalized data set using 99 Pb isotope data on *denarii* and *victoriati* from the recent literature (Albarède et al., 2016; Desaulty et al., 2011; Westner et al., 2020). These data sets are of high quality owing to the excellent control of mass

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

bias allowed for by MC-ICP-MS. Additionally, many analyses were replicated. Applying a precautionary strategy to test the results for both coarse and fine clustering, the number *n* of clusters was set to, respectively, 4 and 12.

A smooth spread of the data is obtained for the 208Pb/206Pb vs 207Pb/206Pb plot and for the 208Pb-normalized plots whether 2D or 3D (Fig. 3 and Fig. 4, Suppl. Fig.1). The results roughly define a ‘mixing triangle’, with three end-members forming the apexes of the triangle. For a maximum number of 12 clusters, these apexes are (i) cluster 1 (red circles), (ii) cluster 7 (dark-green squares), and (iii) cluster 12 (grey squares). All the other clusters fall within this triangle. As indicated by rather high *p*-values (Table 1), group 11 could probably be merged with other groups.

Model ages must not be considered a tool to accurately date the emplacement of ores, which, for most ores, remain non or poorly dated. Nevertheless, for the past six decades Pb model ages have repeatedly
<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Histograms of minting ages for each cluster number and the number of samples in parentheses. The options of a maximum number of clusters *n* = 12 (top, with only the 10 most populated clusters shown) and a minimum number of clusters *n* = 4 (bottom) are both represented. For *n* = 12, only 10 groups are represented because groups 11 and 12 (2 points each) clearly represent unusual or rogue data (see Figs. 5 and 6). Coins from clusters 1, 6, and 10 were minted at the time of the Second Punic War, whereas coins from clusters 2, 3, 5, and 7 were minted several decades after the war was over. This suggests that the original source of silver associated with Cenozoic volcanism mined at the time of the Second Punic War was replaced by silver ores of Mesozoic model age.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 8.** Ternary plot of 1781 Pb isotope literature data on galena ores from the Aegean realm and neighboring regions using a maximum number of 12 clusters (compare with Fig. 3).</figcaption></figure>

proven their chronological usefulness for conformable and many other
types of galenas (Albarede and Juteau, 1984; Cumming and Richards,
1975; Stacey and Kramers, 1975; Stanton and Russell, 1959). Experi-
enced geochemists can usually infer at a glance the age of the tectonic
province, e.g., Alpine, Hercynian, Pan-African, in which galenas are
emplaced from a simple reading of their ²⁰⁶Pb/²⁰⁴Pb ratios. The major
reservation about model ages is whether they date the ore emplacement

or the source material of the ore in question.

The clusters identified by Pb isotopes can be examined in plots of geochemical parameters (Fig. 5). Three end-members, i.e., the extreme points, can be defined. First, the 0–70 Ma model ages of end-member (i) (cluster 1), which we will refer to as the ‘Cenozoic’ component ($T_m$~0 Ma, μ = 9.75, and κ = 3.85). Second, end-member (ii) (cluster 7), which we will refer to as the ‘Mesozoic’ component ($T_m$ = 100 Ma, μ~ 9.8, and κ = 4.02). It is remarkable that the (i)-(ii) edge of the triangle is consistent with the Pb isotopic array defined by Miocene felsic and intermediate volcanic rocks from the southern coast of Andalusia and the Alboran Ridge (Duggen et al., 2004) (Fig. 3 and 4, Suppl. Fig.1). This agreement is better for Pb isotope ratios that do not involve ²⁰⁴Pb, which suggests some issues with analytical mass bias (compare with the ellipses of Fig. 2). Particularly radiogenic Sr in these volcanic rocks, e.g., ⁸⁷Sr/⁸⁶Sr > 0.706, indicates that the Rb/Sr ratio of their source is high and therefore that the rocks contain melt components from the underlying continental crust. End-member (ii) also coincides with the Pb isotope composition of galenas found by Baron et al. (2017) in a Roman exploitation at La Unión, Murcia (Fig. 3 and 4, Suppl. Fig.1). Third, end-member (iii) (cluster 12, 2 samples), which is less populated ($T_m$ > 400 Ma, μ = ~9.8, and κ = 4.02) and therefore involves Paleozoic basement, perhaps signaling the Roman quest for alternative silver sources such as Sierra Morena, the Central Sierra, Sardinia, and Southern Gaul, or simply contamination by Pb from country rocks. Model ages >400 Ma are found for the 118 BCE Licinus Crassus *denarius* from the Narbo mint in Gaul previously highlighted by Desaulty et al. (2011) and for a 112-111 BCE *denarius* discussed by Westner et al. (2020). This end-member will be referred to as the ‘Paleozoic’ component.

Mixing relationships appear as rather elongated clusters particularly
visible in ²⁰⁸Pb/²⁰⁶Pb vs ²⁰⁷Pb/²⁰⁶Pb or ²⁰⁸Pb/²⁰⁷Pb vs ²⁰⁶Pb/²⁰⁷Pb
<figure><img src="image_15.png" /><figcaption>**Fig. 9.** Same as Fig. 8 but as binary orthogonal projections onto selected faces of ternary plots. As the individual data points are too many to be easily readable, contours enclosing 80% of the points in each cluster were drawn. See Supplementary Fig. S3 for a comprehensive plot of all possible pairs of ratios.</figcaption></figure>

plots for *n* = 4 (Fig. 6, Suppl. Fig. 2). The nature of mixing remains ambiguous since such relationships may reflect either a geological process (mixtures of Pb from magmatic or sedimentary ores with Pb from country rocks, or multiple ore populations), an anthropogenic binary blend of ores, flans, or scraps (recycling) of different origins, or both.

For both options (cluster solution 1–12 and cluster solution 1–4), the histograms of minting times by cluster (Fig. 7) add a remarkable element of information. For cluster solution 1–12 and with few exceptions, clusters 1, 6, and 8 ($T_m$ = 0–100 Ma, $\mu$ = 9.70–0.75, and $\kappa$ = 3.85–3.95) comprise coins struck between 250 and 200 BCE, a time interval that overlaps with the Second Punic War, whereas clusters 2, 3, 5, and 7 ($T_m$ = 100–200 Ma, $\mu$ = 9.75–9.80, and $\kappa$ = ~4.0) are dominated by coins

struck at a later time (150-100 BCE). Whether outliers in histograms, such as in cluster 1, reflect a wrong cluster assignment or metal recycling is unclear. *p* values > 0.05 (Table 1) indicate that Pb isotopes from cluster 11 could be merged with clusters 6 or 8 and therefore in both cases populate the Punic War-like cluster even further. Cluster solution 1–4 confirms the dichotomy observed for cluster solution 1–12 (Fig. 7, bottom). As for cluster solution 1–4, clusters 1 and 4 are dominated by post-Second Punic War material, whereas cluster 2 comprises coins struck at the time of the Second Punic War or shortly after. Cluster 3 collects coins with mixed or intermediate isotopic signatures.

The mines with Mesozoic Pb model ages exploited by the Romans seem to have been supplemented or largely replaced by other sources of coinage silver a few decades after the Second Punic War, at which point
<figure><img src="image_16.png" /><figcaption>**Fig. 10.** Same as Fig. 9 with rows 1 to 4 corresponding to $^{204}$Pb to $^{208}$Pb normalization, but within the space of geologically informed parameters. Negative model ages hint at a mantle-derived (magmatic) component. Three supergroups of clusters are easily identified: the *Hercynian supercluster A* (clusters 6, 8, and 10), the *supercluster B* (clusters 4, 5, and 9), and the most populated *Aegean supercluster C* (clusters 1, 2, 3, and 7) stand out (see text). See *Supplemental Fig. S4* for a comprehensive plot of all possible pairs of ratios.</figcaption></figure>

the Roman authorities pushed for finding new ore districts. As suggested by Westner et al. (2020), input of Aegean silver is also an option which could be an echo of the defeat of Macedonians at Pydna in 168 BCE and the Roman ban on gold and silver mining (Livy 45.29.11) the destruction of Carthage, or the fall of Corinth in 146 BCE. This interpretation is difficult to apply to coins with Mesozoic Pb model ages as ores with such model ages are uncommon in the most productive districts of the Aegean world (see below).

## 3.2. A geolocalized Pb isotope database of literature galena samples from the Aegean realm and adjacent regions

This case differs from the previous one by allowing data to be plotted in terms of geographic coordinates. We used a database of 1781 Pb isotope measurements compiled by Delile et al. (2014) and Blichert-Toft

et al. (2016) which assembles data able to be geolocalized from the OXALID database and recent literature. We restricted the present example to the Aegean realm and Asia Minor (19–35°E and 32–45°N). We adopted the toponymy of the time whenever possible. It should be kept in mind that an unknown fraction of these data may be affected by relatively strong mass bias of analytical origin. Local averages are computed with averaging windows of 0.25° × 0.25° and a maximum number of 12 clusters. In addition to the 2D and 3D isotopic diagrams, a map with pie charts was produced.

Overall, the data forms a massive cloud with barely discernible structure (Fig. 8). Fig. 9 and Suppl. Fig. 3 contour fields that enclose 80% of the data belonging to each cluster. The small populations of the smallest clusters (11 and 12) suggest that 12 clusters, regardless of the reference isotope used, adequately describe the full isotopic variability. For the sake of brevity, the case *n* = 4 is not considered here, but the conclusions are in line with those of the case *n* = 12. In contrast to the study of silver coinage from the Roman Republic, the clusters here depend on the reference isotope, though only to a small extent. The reason for such a difference is likely to be ascribed to the abundance of data acquired before the 2000s. Nevertheless, some features are common to all the representations. From the perspective of provenance and geological setting, smooth representations are obtained in the $^{208}$Pb/$^{204}$Pb vs $^{207}$Pb/$^{204}$Pb plot (Fig. 9, 1c), the $^{208}$Pb/$^{207}$Pb vs $^{204}$Pb/$^{207}$Pb plot (Fig. 9, 3b), and the $^{207}$Pb/$^{208}$Pb vs $^{204}$Pb/$^{208}$Pb plot (Fig. 9, 4b), which, with the help of Fig. 2, suggests the presence of at least two different provinces with different $T_m$, µ (U/Pb), and µκ (Th/Pb), whereas κ (Th/U) is less discriminant (see also Fig. 10 and Suppl. Fig. 4).

Taking the $^{208}$Pb-normalized plots (Fig. 9) and the inferred $T_m$, µ, and κ parameters (Fig. 10) as an example, clusters 6 (magenta), 8 (pink), and 10 (grey) form a conspicuous Paleozoic ‘supercluster’ A ($T_m$ > 200 Ma, µ = 9.8, and κ = 4.0), which corresponds to the Paleozoic basement (Hercynian). Supercluster B ($T_m$ = 0–100 Ma, µ = 9.4–9.6, and κ = 3.8–4.0) is composed of clusters 4, 5, and 9 and is best represented in Cyprus. The Aegean supercluster C gathers the most populated clusters 1, 2, and 3, plus cluster 7 ($T_m$ = 0–100 Ma, µ = 9.6–10.0, and κ = 3.85–4.0): $T_m$ model ages are younger (0–100 Ma) than those of supercluster A, while µ values are higher than those of supercluster B. Similar superclusters based on $T_m$ and µ are obtained when the normalizing isotope is changed, but the cluster numbering within the superclusters may differ because of a slightly different clustering tree. Again, *p* values < 0.01 indicate that the clusters are statistically significant. Negative model ages observed in some clusters are associated with a $^{206}$Pb/$^{204}$Pb more radiogenic than modern Pb. These samples are best accounted for by the presence of a mantle-derived Pb component (Hofmann, 1997).

The map of Fig. 11 highlights the overall geodynamic features of the Aegean, in particular the volcanic arc of the Cyclades, and the boundary between Rhodope and the Balkan Range. The Aegean supercluster C includes the major silver-producing districts (Thrace, Macedonia, and Attica). The districts are characterized by ages $T_m$ < 100 Ma, µ values from 9.6 to 9.9, and κ values from 3.80 to 3.95. Clusters 1 and 2 spread over Macedonia, Thrace, and Northwestern Asia Minor (Mysia). Attica and Lavrion account for a substantial part of the well-populated cluster 1. Some plots (e.g., Fig. 8, 4a) suggest that the largest cluster 1 may split up between Thrace and Macedonia on the one hand, and the geologically younger Lavrion on the other hand, which warrants further, more detailed and precise Pb isotope studies of these mining districts. In the Cyclades Islands and Attica, cluster 3, with many negative model ages suggestive of a mantle-derived magmatic component, is well separated from cluster 1 in Attica by the West Cyclades detachment (Grasemann et al., 2012; Melfos and Voudouris, 2017) and closely follows the active and modern volcanic arc all the way to Samos (Jolivet et al., 2013; Nomikou et al., 2013; Pe-Piper and Piper, 2007). Lead isotopes clearly separate ores from Sifnos from ores from Lavrion. In the Cyclades again, clusters 2 and 7 occupy an inner position with respect to cluster 3. The
<figure><img src="image_17.png" /><figcaption>**Fig. 11.** Map with pie charts of Pb isotope compositions from the Aegean realm and adjacent areas averaged over a 0.25° × 0.25° grid and with a maximum number of 12 clusters. Isotopic ratios are normalized to $^{208}$Pb. The most prominent features are the active volcanic arc (cluster 3), the inactive arc with prominent pre-200 Ma model ages (Peloponnese and Crete, clusters 6, 8, and 10), Cyprus (clusters 5 and 9), and the most populated, Ag ore-rich Aegean districts of Thrace, Macedonia, and Attica (clusters 1 and 2). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

Peloponnese and Crete are part of the inactive arc with a dominant contribution of old model age clusters 8 and 10. Clusters 4, 6, and 8 with Mesozoic and Late Paleozoic ages (100–300 Ma) dominate north of Thrace and Macedonia (Balkan proper). Cyprus is dominated by clusters 5 and 9 (supercluster B), and is accounted for by volcanogenic massive sulfide (VMS) Cu ores deposited by hydrothermal fluids on recent seafloor (Hannington, 2014) distinguishing them from the majority of mineralizations around the Aegean.

As discussed above, if the μ values of the Mesozoic end-member present in the coinage of the Roman Republic are compatible, the κ values are likely to be incompatible (κ > 4.0 for the Roman Republic coins), but this broad statement warrants the production of isotopic data conforming to modern standards of precision and accuracy. If Aegean silver represents this Mesozoic component, it must be from districts for which the production is either small (e.g., Sifnos), overlooked, or underestimated by historians.

## 4. Conclusions

We have reviewed all the possible representations of Pb isotope data in 2D and 3D space and have predicted which information on the geological environment can be derived from each plot. We encourage analysts to (1) measure all four Pb isotopes on their samples, (2) use consistent representations of Pb isotope ratios and select that or those combinations of plots that provide the most informative perspective, and (3) represent growth curves and μ and κ lines on at least part of their figures. The geological framework of Pb isotope distributions is an

invaluable source of information and sheds new light on the origins of ores used to manufacture silver artifacts and coins.

We further devised an algorithm that predicts the number of isotopic clusters necessary to fulfill the simple condition that variance within each cluster is minimized while inter-cluster variance is maximized. We found that, if the geological context is reasonably well understood, defining clusters of Pb isotope compositions in the proper 3D space provides a data-driven understanding of metal provenance independent of assumptions based on literary sources.

A first application of the new method presented here to silver *denarii* and *victoriati* from the Roman Republic using literature data (Albarède et al., 2016; Desaulty et al., 2011; Westner et al., 2020) defines clusters that may represent three end-members corresponding to Betic, Hercynian, and a Mesozoic mining district. Wherever mixing is present, however, it is not clear whether it is of natural or anthropogenic origin.

A second application is the case of Aegean ores which separate Pb components from the Hercynian basement, the recent Aegean tectonic province, and Cyprus. A strong connection can be made between the background geology and Pb isotope compositions of galenas, which confirms that Pb isotopes betray the underlying geodynamics of ore formation more than a specific provenance.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
<figure><img src="image_18.png" /></figure>

(caption on next column)

<caption>**Fig. 12.** Lead isotope ratios measured for Northern Italian awls, pins, and borers of the 4th millennium BCE (Artioli et al., 2020). The figures have been redrawn to emphasize the different references (1–10), which can be found in the original publication. Blue line: natural radiogenic ingrowth. Red line: mass-dependent fractionation resulting from low-temperature geological processes, incomplete yield during Pb purification, or improperly corrected instrumental mass bias. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</caption>

## Acknowledgments

This is a contribution of Advanced Grant 741454-SILVER-ERC-2016-ADG ‘Silver Isotopes and the Rise of Money’ given to FA by the European Research Council. Gilberto Artioli is gratefully acknowledged for sharing the Pb isotope database used in his publication. We gratefully acknowledge constructive suggestions by Yigal Erel and an anonymous reviewer.

## Appendix A

White’s (2015) textbook is a technically proficient reference for Pb isotope geochemistry, while Pollard’s (2018) text provides a more accessible perspective oriented towards archeometry. For measured isotope ratios subscripted 0, the model age $T_m$ and the parameters $\mu$ and $\kappa$ obey the following equations

$$x = \left( \frac{^{206}\text{Pb}}{^{204}\text{Pb}} \right)_0 = \left( \frac{^{206}\text{Pb}}{^{204}\text{Pb}} \right)_{SK,T_0} + \mu (e^{\lambda_{238U}/T_0} - e^{\lambda_{238U}/T_m})$$

$$y = \left( \frac{^{207}\text{Pb}}{^{204}\text{Pb}} \right)_0 = \left( \frac{^{207}\text{Pb}}{^{204}\text{Pb}} \right)_{SK,T_0} + \frac{\mu}{137.79} (e^{\lambda_{235U}/T_0} - e^{\lambda_{235U}/T_m})$$

$$z = \left( \frac{^{208}\text{Pb}}{^{204}\text{Pb}} \right)_0 = \left( \frac{^{208}\text{Pb}}{^{204}\text{Pb}} \right)_{SK,T_0} + \mu \kappa (e^{\lambda_{232Th}/T_0} - e^{\lambda_{232Th}/T_m})$$

where $T_0 = 3.8$ Ga is the beginning of the second stage, subscript SK, $T_0$ represents the value of each ratio at that age (Albarede and Juteau, 1984; Stacey and Kramers, 1975), and $\lambda_i$ is the decay constant of nuclide i. Each triplet $T_m$, $\mu$, and $\kappa$ corresponds to a triplet of isotopic ratios. Other ratios with a different normalization isotope can be used, but the triplet $T_m$, $\mu$, and $\kappa$ remains unchanged.

For young model ages, typically $\le 600$ Ma, radiogenic ingrowth is approximately linear with time, which leads to simple constraints on the arrays formed by isotopic data. In the $^{207}\text{Pb}/^{204}\text{Pb}$ vs $^{206}\text{Pb}/^{204}\text{Pb}$ diagram, the slope of isochrons is approximately equal to 0.4. In the $^{208}\text{Pb}/^{204}\text{Pb}$ vs $^{206}\text{Pb}/^{204}\text{Pb}$ diagram, this slope should not deviate much from about 1.0. Slopes defined by an array of points deviating significantly from these values in both plots signal either the inmixing of old Pb or some degree of mass-dependent mass fractionation in the ore-forming process or during analytical processing.

The effect of radiogenic ingrowth is readily distinguished from the bias created by analytical issues and low-temperature processes. Incomplete yields during Pb purification on ion-exchange resins (Baker et al., 2005; Blichert-Toft et al., 2003) and inadequately corrected instrumental mass bias (Albarède et al., 2004) spread ratios along arrays with predictable slopes in both the $^{207}\text{Pb}/^{204}\text{Pb}$ vs $^{206}\text{Pb}/^{204}\text{Pb}$ and $^{208}\text{Pb}/^{204}\text{Pb}$ vs $^{206}\text{Pb}/^{204}\text{Pb}$ plots. For sake of simplicity, let us assume that mass-dependent fractionation is linear, i.e., each measured ratio $r$ (meas) departs from its original value $R$ (true) according to $r = R(1 + \Delta m f)$, where $\Delta m$ is the mass difference between the masses of isotopes at the numerator and denominator, and $f$ is the fractionation factor per a.m.u. (atomic mass unit). The term ‘mass fractionation’ itself refers to the relative change of the isotope ratio during the measurement process and is equal to $(r-R)/R = \Delta m f$. Technically, there is no actual physics in this
equation (Albarède et al., 2004), which is no more than a first-order expansion of *r* variations with *Δm*: to order one, mass fractionation for a mass difference *Δm* of 2 is twice the mass fractionation for a mass difference of 1. For example, the slope of mass fractionation lines in the $y = 207\text{Pb}/204\text{Pb}$ vs $x = 206\text{Pb}/204\text{Pb}$ diagram is

$$ \frac{y - Y}{x - X} = \frac{207 - 204}{206 - 204} \times \left( \frac{207\text{Pb}/204\text{Pb}}{206\text{Pb}/204\text{Pb}} \right)_{\text{true}} = 1.5 \left( \frac{207\text{Pb}}{206\text{Pb}} \right)_{\text{true}} $$

while in the $z = 208\text{Pb}/204\text{Pb}$ vs $x = 206\text{Pb}/204\text{Pb}$ diagram the slope is

$$ \frac{z - Z}{x - X} = \frac{208 - 204}{206 - 204} \times \left( \frac{208\text{Pb}/204\text{Pb}}{206\text{Pb}/204\text{Pb}} \right)_{\text{true}} = 2 \left( \frac{208\text{Pb}}{206\text{Pb}} \right)_{\text{true}} $$

where lower case and upper case variables represent the measured and 'true' values of the ratios, respectively. The values of these slopes are $1.5 \times 15.5/18.7 \approx 1.2$, i.e., three times the slope of isochrons, and $2 \times 39/18.7 \approx 4$, i.e., four times the isochron slope (Table A1), respectively. These numbers do not change with the age of ore sources considered in the Mediterranean region (<600 Ma), but may be different in very old (i.e., Proterozoic and Archean) terranes. Choosing an exponential rather than a linear fractionation law results in only very minor differences in slopes.

Finally, insufficient precision due to a small ion beam signal on the 'small' isotope $204\text{Pb}$ and mass interference due to the presence of $204\text{Hg}$, both leading to poor sensitivity of $204\text{Pb}$, may occasionally affect isotopic ratios whenever $204\text{Pb}$ is used as the reference isotope. This problem is not uncommon for Pb isotope data acquired on ancient thermal ionization mass spectrometers and with more recent quadrupole and single-collector ICP-MS. This issue can be identified by the large spread of ratios of the major isotopes $206\text{Pb}$, $207\text{Pb}$, and $208\text{Pb}$ relative to $204\text{Pb}$ compared to ratios involving only the major isotopes. This is the case for the standard plots of $y = 207\text{Pb}/204\text{Pb}$ vs $x = 206\text{Pb}/204\text{Pb}$ and $z = 208\text{Pb}/204\text{Pb}$ vs $x = 206\text{Pb}/204\text{Pb}$, in which the spread of $y/x = 207\text{Pb}/206\text{Pb}$ and $z/x = 208\text{Pb}/206\text{Pb}$, respectively, remains quite narrow in spite of 'noisy' $x$, $y$, and $z$.

These estimates are summarized in Table A1. Lead isotope arrays with non-geological slope pairing according to this table therefore should be considered with utmost caution.

Examples of slopes of arrays resulting from either mass-dependent fractionation or radiogenic ingrowth can be found in Fig. 1 of Artioli et al.'s (2020) review, which has been redrawn here as Fig. 12. References associated with each symbol color can be found in the original publication.

## Appendix B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jas.2020.105194.

<caption>Table A1</caption>

Slopes predicted for isochrons (radiogenic ingrowth) vs other processes. Mass-dependent fractionation may result from natural processes at low temperature or from analytical artifacts, such as incomplete yield of Pb purification and inadequate correction of mass-spectrometer mass bias. Small signals (poor counting statistics) or isobaric interferences at mass 204 also result in spurious isotopic arrays. Similar predictions can be made for other pairs of isotopic ratios.

<table><thead><tr><th>x</th><th>206Pb/204Pb</th><th>206Pb/204Pb</th><th>207Pb/206Pb</th></tr><tr><th>y</th><th>207Pb/204Pb</th><th>208Pb/204Pb</th><th>208Pb/206Pb</th></tr></thead><tbody><tr><td><strong>Process</strong></td><td></td><td></td><td></td></tr><tr><td>Radiogenic ingrowth (isochron)</td><td>~0.05</td><td>~1.3</td><td>1.0</td></tr><tr><td>Low-temperature process, low separation yield, mass bias issues</td><td>~1.2</td><td>~4</td><td>5.0</td></tr><tr><td>Counting statistics, isobaric interference at mass 204</td><td>~0.8</td><td>~2.1</td><td>2.5</td></tr></tbody></table>

## References

Albarède, F., Blichert-Toft, J., Rivoal, M., Telouk, P., 2016. A glimpse into the Roman finances of the Second Punic War through silver isotopes. Geochemical Perspective Letters 2, 127–137.

Albarede, F., Desaulty, A.M., Blichert-Toft, J., 2012. A geological perspective on the use of Pb isotopes in Archaeometry. Archaeometry 54, 853–867.

Albarede, F., Juteau, M., 1984. Unscrambling the lead model ages. Geochem. Cosmochim. Acta 48, 207–212.

Albarède, F., Telouk, P., Blichert-Toft, J., Boyet, M., Agranier, A., Nelson, B., 2004. Precise and accurate isotopic measurements using multiple-collector ICPMS. Geochem. Cosmochim. Acta 68, 2725–2744.

Arribas, A., Tosdal, R.M., 1994. Isotopic composition of Pb in ore-deposits of the Betic Cordillera, Spain - origin and relationship to other European deposits. Econ. Geol. 89, 1074–1093.

Artioli, G., Canovaro, C., Nimis, P., Angelini, I., 2020. LIA of Prehistoric Metals in the Central Mediterranean Area: a Review. Archaeometry (in press).

Baker, J., Stos, S., Waight, T., 2005. Lead isotope analysis of archaeological metals by multiple-collector inductively coupled plasma mass spectrometry. Archaeometry 48, 45–56.

Baron, S., Rico, C., Antolinos Martin, J.A., 2017. The Cabeko del Pino workshops complex (Sierra Minera de Cartagena-La Union, Murcia) and the mining activities organisation at Carthago Noua at the end of the Roman Republic. Cross-

contributions of archaeology and geochemistry. Arch. Español Arqueol. 90, 147–169.

Baron, S., Souhassou, M., Fauvelle, F.X., 2020. Medieval silver production around sijilmasa, Morocco. Archaeometry 62, 593–611.

Birch, T., Kemmers, F., Klein, S., Seitz, H., Höfer, H., 2020a. Silver for the Greek Colonies: Issues, Analysis and Preliminary Results from a Large-scale Coin Sampling Project. In: Metallurgy in Numismatics, vol. 6.

Birch, T., Westner, K., Kemmers, F., Klein, S., Höfer, H.E., Seitz, H.M., 2020b. Retracing Magna Graecia's silver: coupling lead isotopes with a multi-standard trace element procedure. Archaeometry 62, 81–108.

Blichert-Toft, J., Weis, D., Maerschalk, C., Agranier, A., Albarède, F., 2003. Hawaiian hot spot dynamics as inferred from the Hf and Pb isotope evolution of Mauna Kea volcano. G-cubed 4. https://doi.org/10.1029/2002GC000340.

Blichert-Toft, J., Delile, H., Lee, C.T., Stos-Gale, Z., Billström, K., Andersen, T., Hannu, H., Albarède, F., 2016. Large-scale tectonic cycles in Europe revealed by distinct Pb isotope provinces. G-cubed 17, 3854–3864.

Cheilletz, A., Levresse, G., Gasquet, D., Azizi-Samir, M., Zyadi, R., Archibald, D.A., Farrar, E., 2002. The giant Imiter silver deposit: neoproterozoic epithermal mineralization in the Anti-Atlas, Morocco. Miner. Deposita 37, 772–781.

Cumming, G.L., Richards, J.R., 1975. Ore lead isotope ratios in a continuously changing earth. Earth Planet Sci. Lett. 28, 155–171.

De Ceuster, S., Degryse, P., 2020. A 'match-no match' numerical and graphical kernel density approach to interpreting lead isotope signatures of ancient artefacts. Archaeometry. https://doi.org/10.1111/arcm.12552.
Delile, H., Blichert-Toft, J., Goiran, J.-P., Keay, S., Albarède, F., 2014. Lead in ancient Rome's city waters. Proc. Natl. Acad. Sci. Unit. States Am. 111, 6594-6599.

Delile, H., Blichert-Toft, J., Goiran, J.-P., Stock, F., Arnaud-Godet, F., Bravard, J.-P., Brückner, H., Albarede, F., 2015. Demise of a harbor: a geochemical chronicle from Ephesus. J. Archaeol. Sci. 53, 202-213.

Desaulty, A.M., Telouk, P., Albalat, E., Albarede, F., 2011. Isotopic Ag-Cu-Pb record of silver circulation through 16th-18th century Spain. Proc. Natl. Acad. Sci. U.S.A. 108, 9002-9007.

Duggen, S., Hoernle, K., van den Bogaard, P., Harris, C., 2004. Magmatic evolution of the Alboran region: the role of subduction in forming the western Mediterranean and causing the Messinian Salinity Crisis. Earth Planet Sci. Lett. 218, 91-108.

Erel, Y., Harlavan, Y., Blum, J.D., 1994. Lead isotope systematics of granitoid weathering. Geochem. Cosmochim. Acta 58, 5299-5306.

Eshel, T., Erel, Y., Yahalom-Mack, N., Tirosh, O., Gilboa, A., 2019. Lead isotopes in silver reveal earliest Phoenician quest for metals in the west Mediterranean. Proc. Natl. Acad. Sci. Unit. States Am. 116, 6007-6012.

Frei, R., Kamber, B.S., 1995. Single mineral PbPb dating. Earth Planet Sci. Lett. 129, 261-268.

Fujii, T., Moynier, F., Agranier, A., Ponzevera, E., Abe, M., 2011. Nuclear field shift effect of lead in ligand exchange reaction using a crown ether. Proceedings in Radiochemistry A Supplement to Radiochimica Acta 1, 387-392.

Galer, S.J., 1999. Optimal double and triple spiking for high precision lead isotopic measurement. Chem. Geol. 157, 255-274.

García-Hernández, M., López-Garrido, A., Rivas, P., Sanz de Galdeano, C., Vera, J.A., 1980. Mesozoic Palaeogeographic Evolution of the External Zones of the Betic Cordillera.

Getty, S.R., Depaolo, D.J., 1995. Quaternary geochronology using the U- Th- Pb method. Geochem. Cosmochim. Acta 59, 3267-3272.

Grasemann, B., Schneider, D.A., Stöckli, D.F., Iglseder, C., 2012. Miocene bivergent crustal extension in the Aegean: evidence from the western Cyclades (Greece). Lithosphere 4, 23-39.

Hannington, M.D., 2014. 13.18 - volcanogenic massive sulfide deposits. In: Holland, H. D., Turekian, K.K. (Eds.), Treatise on Geochemistry, second ed. Elsevier, Oxford, pp. 463-488.

Harlavan, Y., Erel, Y., Blum, J.D., 2009. The coupled release of REE and Pb to the soil labile pool with time by weathering of accessory phases, Wind River Mountains, WY. Geochem. Cosmochim. Acta 73, 320-336.

Hofmann, A.W., 1997. Mantle geochemistry: the message from oceanic volcanism. Nature 385, 219-229.

Jolivet, L., Faccenna, C., Huet, B., Labrousse, L., Le Pourhiet, L., Lacombe, O., Lecomte, E., Burov, E., Denèle, Y., Brun, J.-P., 2013. Aegean tectonics: strain localisation, slab tearing and trench retreat. Tectonophysics 597, 1-33.

Kettenring, J.R., 2006. The practice of cluster Analysis. J. Classif. 23, 3-30.

Melfos, V., Voudouris, P., 2017. Cenozoic metallogeny of Greece and potential for precious, critical and rare metals exploration. Ore Geol. Rev. 89, 1030-1057.

Nomikou, P., Papanikolaou, D., Alexandri, M., Sakellariou, D., Rousakis, G., 2013. Submarine volcanoes along the Aegean volcanic arc. Tectonophysics 597, 123-146.

Pe-Piper, G., Piper, D., 2007. Late Miocene igneous rocks of Samos: the role of tectonism in petrogenesis in the southeastern Aegean, 291. Geological Society, London, Special Publications, pp. 75-97.

Pollard, M., 2018. Beyond Provenance. Leuven University Press.

Stacey, J.S., Kramers, J.D., 1975. Approximation of terrestrial lead isotope evolution by a two-stage model. Earth Planet Sci. Lett. 26, 207-221.

Stanton, R., Russell, R.D., 1959. Anomalous leads and the emplacement of lead sulfide ores. Econ. Geol. 54, 588-607.

Stos-Gale, Z.A., Gale, N.H., 2009. Metal provenancing using isotopes and the Oxford archaeological lead isotope database (OXALID). Archaeological and Anthropological Sciences 1, 195-213.

Székely, G.J., Rizzo, M.L., 2005. Hierarchical clustering via joint between-within distances: extending ward's minimum variance method. J. Classif. 22, 151-183.

Székely, G.J., Rizzo, M.L., 2013. Energy statistics: a class of statistics based on distances. J. Stat. Plann. Inference 143, 1249-1272.

Tera, F., Wasserburg, G.J., 1972. U-Th-Pb systematics in three Apollo 14 basalts and the problem of initial Pb in lunar rocks. Earth Planet Sci. Lett. 14, 281-304.

Wessel, P., Smith, W.H.F., Scharroo, R., Luis, J., Wobbe, F., 2013. Generic mapping tools: improved version released. EOS Trans. AGU 94, 409-410.

Westner, K.J., Birch, T., Kemmers, F., Klein, S., Höfer, H.E., Seitz, H.-M., 2020. Rome's rise to power. Geochemical analysis of silver coinage from the Western Mediterranean (4th to 2nd centuries BCE). Archæometry 62, 577-592.

White, W.M., 2015. Isotope Geochemistry. John Wiley & Sons.

Yang, S., Liu, Y., 2015. Nuclear volume effects in equilibrium stable isotope fractionations of mercury, thallium and lead. Sci. Rep. 5, 12626.

Yigit, O., 2009. Mineral deposits of Turkey in relation to Tethyan metallogeny: implications for future mineral exploration. Econ. Geol. 104, 19-51.

Żabiński, G., Gramacki, J., Gramacki, A., Miśta-Jakubowska, E., Birch, T., Disser, A., 2020. Multi-classifier majority voting analyses in provenance studies on iron artefacts. J. Archaeol. Sci. 113, 105055.