Bayesian mixing models as a tool to explore Bronze Age bitumen trade from Tell Lashkir (Erbil, Iraq)

Bitumen from north-eastern Mesopotamian seeps was exchanged across the Tigris and the Euphrates in antiquity. Nonetheless, recent research proposed a decline in its trade during the Bronze Age alongside the boost of central Mesopotamian (Hit) seeps trade. This paper presents new provenance studies from ancient bitumen recovered at Tell Lashkir, a 4th-2nd millennium BC settlement between the lower and upper Zab Tigris tributaries. The novel application of Bayesian mixing models to biomarker and isotopic data has narrowed down the identification of the bitumen putative sources and improved the apportionment of mixtures. This placed Tell Lashkir within the pathway of bitumen active trade routes in the Bronze Age. Thus, the population inhabiting Tell Lashkir participated in exchange networks which reached southern Mesopotamia, and traded with natural resources most probably originating from central and north-eastern Mesopotamia. From a methodological standpoint, the appraisal of the Bayesian mixing models in our study with modern and archaeological reference samples has provided a better understanding of the strengths and limitations of such an approach to apportion bitumen sources from several Bronze age territories.


Introduction
In ancient times, bitumen was a valuable commodity used for its stickiness, ignition, and waterproofing properties.A testimony of its importance for the first states in Mesopotamia was its use in monumental architecture, shipbuilding, water management, art pieces, magic and medicine (Forbes, 1964).The term for bitumen and its extensive derived vocabulary appear in cuneiform tablets attesting both commercial exchanges and political fights for its control (Stol, 2012).This natural resource played a significant role in the economic and political developments of the Bronze Age as it was not only exchanged in either raw or processed form, but its use in boat construction facilitated trade itself.The molecular and isotopic analyses of ancient bitumen samples have revealed that this product circulated across south-west Asia from early prehistory to medieval times (Connan, 1999;Connan and Van de Velde, 2010).This is thanks to the chemical characterization of putative sources (i.e., oil seeps) and the identification and comparison of their geochemical fingerprints.
As suggested by recent studies in the Ubaid period, the three main regions involved in the long distance trade of bitumen were southwestern Iran and central and northern Mesopotamia (Connan and Van de Velde, 2010).Nevertheless, during the Uruk period, bitumen from Hit (central Mesopotamia) also became prevalent across the Euphrates, thus possibly replacing previous trade networks at a time when the construction of monumental architecture would have demanded increasing quantities of this product (Connan and Van de Velde, 2010;Schwartz and Hollander, 2016).Texts from the Ur III period and molecular analyses suggest that, through the third millennium BC., Madga bitumen from north-eastern Mesopotamia still reached southern cities (Connan and Van de Velde, 2010).Nevertheless, the scarcity of molecular data across the Tigris in the Bronze Age (Connan et al., 2013(Connan et al., , 2018) ) prevents obtaining a clear picture of the extension of the bitumen trade networks at the moment of the formation of the first empires in history.Did settlements near the Tigris bank or across its tributaries use bitumen presenting the same geochemical fingerprints as those from the southern cities which they allegedly served?Were bitumen types more frequent in settlements across the Euphrates or the Persian Gulf present in the northern Tigris?The range of uses bitumen had in north-eastern Mesopotamia, where the Neo-Assyrian empire would develop during the Iron Age, might have been different, as building materials such as stone were more readily available than in the south (Forbes, 1964;Moorey, 1994).This, in turn, could have affected the preferred types of bitumen throughout the region.
To contribute further to the current understanding of Bronze Age bitumen trade across the Tigris this paper reports the study of the bitumen remains recovered at Tell Lashkir (hereafter TL).This is an archaeological site located roughly 10 km east from ancient Arbela and 100 km north from modern Kirkuk.Our study is the first that tackles the origins of ancient bitumen from the land around the lower and the upper Zab rivers, and offers the opportunity to understand better the structure of the exchange networks which supported the demand of this natural resource by the bustling Sumerian and Akkadian cities and states.
To identify the most probable sources of bitumen, we undertook qualitative and semi-quantitative molecular and isotopic analyses of black amorphous visible residues inside pottery vessels using gas chromatography coupled to Mass Spectrometry (GC-MS), and Elemental Analysis coupled to Isotope Ratio Mass Spectrometry (EA-IRMS).For the first time in the study of ancient bitumen, the identification of putative sources and its comparison with other published results was performed through Bayesian mixing models.To date, these have been successfully applied to quantify the proportions of different foodstuffs in archaeological pottery residues (Lucquin et al., 2018;Courel et al., 2020) and to study the composition of past human diets (Fernandes et al., 2015).This statistical tool was selected because of its capability to evaluate the possible existence of mixtures from different areas; its ability to suggest the possible broad location of unknown sources, and, to treat both reference data and results as probability distributions carrying a degree of uncertainty (Fernandes et al., 2018).To strengthen these claims, the validity of the model deployed in this study was tested using a significant corpus of published modern references from a single point of origin and previously identified archaeological samples.
In consequence, the four objectives of this study are: a) To test the potential of Bayesian models to better distinguish and apportion bitumen remains with different geographic origins.b) To assess the variability of potential bitumen fingerprint types within TL. c) To study the similarities between TL bitumen and archaeological bitumen from other contemporary sites.d) To estimate the possible original source of TL's archaeological bitumen.

Materials and methods
The site of Tell Lashkir (Fig. 1), is a chronologically long lived settlement with indices of Neolithic presence and clear Uruk, Ninevite 5 (sectors 3 and 6), Middle Bronze Age and Iron Age phases (Gómez-Bach et al., 2018, 2019;Molist et al., 2019a;2019b).Excavations at TL started in 2015 thanks to a collaboration between the Autonomous University of Barcelona, The Directorate of Antiquities of the Kurdistan regional Fig. 1.Location of Tell Lashkir (star), other archaeological sites mentioned in the text (rhomb) and the approximate regions of origin of the modern bitumen references used in the Bayesian models.Source of the data: (Connan and Nishiaki, 2003;Connan et al., 2005Connan et al., , 2013Connan et al., , 2018)).
government and the department of archaeology in the Salahaddin University (Molist, 2015).Since the start of the excavation, significant quantities of bitumen-bearing potsherds were detected in the site (Breu et al., 2019).Hitherto, four excavation campaigns have certified the presence of bitumen in almost all sectors, especially at the Ninevite 5 period, at the end of the Bronze Age and in the Iron Age (Molist et al., 2019b).
This paper presents the results of an initial study which included a total of 12 samples.Five samples were recovered from Sector 1; two samples were from sector 3; four samples from sector 6; and one sample from sector 11.They were all amorphous black visible residues attached to the interior of pottery sherds (see images of the samples in supplementary materials 1).Regarding their context of discovery, Sectors 1, 3 and 11 presented walls and architectural features associated with domestic areas.In Sector 1, placed at the start of the Middle Bronze Age (MBA), four samples were taken from a single large building presenting several storage areas constructed after 2059-1900 2σ cal BC and one sample originated from layers dated after this date.Sector 3 (Ninevite 5, 2887-2636 2σ cal BC, Molist et al., 2019a), containing up to six architectural phases, presented vestiges of earthen architecture related to rectangular buildings.In this sector, samples were collected from the most recent occupation phase.Finally, samples from sector 6 were obtained from layers containing exterior areas related to a nearby human installation placed in the first half of the third millennium cal.BC (Early, Middle Ninevite 5, Early Jezirah I-II, Amuq H and G (Molist et al., 2019b).

Extraction and analysis of archaeological bitumen
Roughly 0.5 g of the black amorphous crusts adhered to internal vessel walls were removed and ground with a pestle and mortar.The resulting samples were extracted with ca. 5 ml of dichloromethane.After 15 min of sonication and 10 min of centrifugation at 2000 rpm, the supernatant was removed.This procedure was repeated 3 times and the three extracts were combined.They were then dried and extracted with 5 ml of hexane.After 15 min of sonication and 10 min of centrifugation at 2000 rpm the supernatants were removed.This procedure was repeated 3 times.The precipitates (asphaltenes) were analysed by EA-IRMS.Lastly, the compound classes in the hexane soluble fraction were separated into the saturated (first), aromatic (second) and polar (third) lipid fractions using a silica gel column and a range of solvents (Bastow et al., 2007).
The isotopic analyses of the asphaltene fractions were undertaken using a Thermo FlashEA 1112 Series elemental analyser coupled to a Thermo Delta V Advantage isotopic ratio mass spectrometer with a Conflo III interface (EA-IRMS).The δ 13 C measurements are reported relative to the V-PDB standard.

Statistical analyses
A series of non-weighted concentration-independent Bayesian mixing models were built using the FRUITS v3.1beta software (Fernandes et al., 2018) to quantify the relative amounts of different bitumen types to a hypothetic mixture defined by specific geochemical ratios.An underlying assumption of this approach is that the source of any compounds extracted from the visible amorphous residues could only be the same components in the bitumen references.In all models, a Markov Chain Monte Carlo (MCMC) method using Gibbs sampling and the Metropolis-Hastings algorithm was run using the BUGS software.It performed 10000 burn-in steps and, afterwards, ran for 40000 iterations.Following Fernandes et al., (2018) the general minimum uncertainty of the model was set at 0.001 to account for possible unknown sources of uncertainty.
Several Bayesian mixing models were deployed to assess data from this study.The first model detected variability within bitumen samples from our site; the second evaluated the similarity of the geochemical fingerprints from TL to other contemporary archaeological sites, and, the third explored the possible geographic location of TL's bitumen sources.In all cases, the model outputs were interpreted as the percentage of the total bitumen by weight and no informative priors were used.
In the first model, the 8 measured geochemical ratios (Ts/Tm, Ts/ 30αβ, Tm/30αβ, 29αβ/30αβ, 30βα/30αβ, 31αβS/30αβ, the 30G/30αβ and 29ββ/29αα) were used as proxies and were described by their median and their standard deviation.Reference values for different bitumen types were selected from representative samples in TL after statistical evaluation.This included a Principal Component Analysis (PCA) using a correlation matrix and hierarchical clustering using a Euclidean similarity index and a paired group (UPGMA) algorithm (Past v4.03).When evaluating posterior distributions, unknown samples were attributed to the type most probably contributing the highest to the mixture.
In the second model, the frequently reported Ts/Tm, 29αβ/30αβ and 30G/30αβ ratios were used as proxies.The reference values and interpretation guidelines from the first model were maintained.
In the third model, potential sources were defined using the values of δ 13 C Asphaltenes , and the Ts/Tm and the G/C 30 αβH biomarker ratios as proxies.These were chosen because of their widespread availability in published studies of modern reference samples.To obtain the source values required by FRUITS (Table 1), we evaluated 94 modern published bitumen samples (Fig. 2) from Central, Southern, North-eastern and North-western Mesopotamia (Connan et al., 1998(Connan et al., , 2005(Connan et al., , 2006b(Connan et al., , 2013(Connan et al., , 2018;;Connan and Nishiaki, 2003) (see supplementary information 3).
The probability that a source contributed more than 50% to the mixture was calculated by integrating the probability curve between the 50% and 100% percentiles in the posterior distributions produced by the Bayesian mixing model.Next, a Linear Discriminant Analysis (LDA) was applied to these results.This step corrected biases undervaluing or overvaluing the probability that certain regions were the main source of bitumen in the mixture.The LDA found the region with the shortest Mahalanobis distance to each sample.This distance was calculated using a pooled within-group covariance matrix (Hammer et al., 2001)  PAST v4.03 statistical software.
To test the accuracy of the third model, 57 published modern reference samples and 86 previously published archaeological samples were analysed as unknowns.In the first case, 76.8% of the results were identical to those published in the literature.In the second case, the rate of success rose to 83.7%.The main sources of error resulted from southern and central Mesopotamian samples mistakenly attributed to north-eastern sources (12.8%)(see supplementary information 5 for more details).Given the high degree of agreement between the references and the outputs and bearing in mind the limitations of the model, well-preserved bitumen residues from TL were submitted to the same analysis.
Finally, the intersections of posterior probability distributions were calculated to quantify the total amount of a bitumen type in a group of samples.This, P(A ∩ B) = P(A)⋅P(B), has been successfully applied to both palaeodieary studies (Phillips et al., 2005(Phillips et al., , 2014;;Stock et al., 2018) and the analysis of lipid residues in archaeological artefacts (Breu et al., 2021).In this case, the overall quantity of TL's bitumen types in the early third millennium and early second millennium was calculated and used to create Figs. 5 and 7.

Molecular and isotopic data
The 12 analysed amorphous visible residues contained well-defined hopane fingerprints and 11 of them also included steranes (Fig. 3) (see additional chromatograms in supplementary information 2).Their identification confirms the petrogenic origin of the samples, which had been preliminarily labelled as bitumen based on their black colour and amorphous texture.
The outcome of the PCA and the groups resulting from the hierarchical clustering analysis at height = 0.26 suggested the existence of four distinct geochemical fingerprints within Tell Lashkir (Fig. 4).Furthermore, as differences between these four groups were maximised in samples 8, 9, 10 and 12, these were used to define four dominant fingerprint types: Type A (sample 10): it contained the lowest 29ββ/29αα and 31αβS/ 30αβ values combined with the highest 29αβ/30αβ and Tm/30αβ values.
In summary, the hopanes eluting before 30αβ in types D and C presented low ratios while those from compounds eluting afterwards were higher.The opposite was observed in the A and B geochemical fingerprints, and it represented one of the main sources of variation between archaeological bitumen in TL.The Ts/Tm ratio further differentiated the D and C from the A and B types, and the 29αβ/30αβ hopane and the 29ββ/29αα sterane ratios were the major sources of variation between D and C, and between A and B. The δ 13 C isotopic ratios of the asphaltene fractions did not aid in distinguishing between types A, C and D. Its only

Table 2
The sample's archaeological context (SU: stratigraphic unit) and the results of the GC-MS and EA-IRMS analyses with their associated standard deviation.significant variation corresponded with type B's asphaltene δ 13 C isotopic values being 0.5‰-0.8‰enriched compared with the other types.
To validate the groupings suggested by the hierarchical clustering analysis and to explore the presence of mixtures, the first Bayesian mixing model incorporated samples 8, 9, 10 and 12 as reference values for the A, B, C and D types, respectively.Results indicated that type A bitumen was the most abundant in the assemblage followed by D, B and C respectively (Fig. 5).Several samples could have resulted from a mixture of two types (Table 4).Sample 1 would have been majorly conformed by type B bitumen with a minor input of type A. Alternatively, sample 5 may have contained major amounts of type A and minor inputs of type B. Finally, sample 7 was interpreted as originating from a mixture of equal parts type D and A (Fig. 6).

Exploring trade through published analyses of archaeological bitumen
To place TL within the Mesopotamian Bronze Age trade networks, the presence of the A, B, C and D geochemical fingerprints was evaluated in 33 published bitumen remains from five archaeological sites located in the 2500-1100 cal.BC time bracket: Kavus ¸an Höyük (Connan et al., 2013), Salat Höyük (Connan et al., 2018), Umm an Namel, Tell F6 (Connan and Carter 2007) and Ra'S al-Jinz (Connan et al., 2005) (See supplementary information 6).Contemporary data from Bas ¸hur Höyük (Connan et al., 2018) and Tell Abraq (Van de Velde et al., 2017), was not modelled as their studies determined origins (Iran) outside the range of sources contemplated by the Bayesian mixing model.In this model, the three hopane ratios more consistently reported in the literature (Ts/Tm, 29αβ/30αβ and 30G/30αβ) were used.Nonetheless, accurate differentiation between types A and B was not clear given their similar Ts/Tm, 30G/30αβ and 29αβ/30αβ ratios.In the absence of other sources of discrimination such as the 29ββ/29αα sterane ratio or the Tm/30αβ and the 31αβS/30αβ hopane ratios, these two types were studied as indicative of the same source.The reference values for the A/B group were thus constructed by averaging samples 8 and 10.
The results of the second Bayesian mixing model indicate that bitumen similar to TL's type D was only present in sites located around the Persian Gulf.Upriver from TL, the model considered the probability for the presence of this geochemical fingerprint to be below 0.1% in Bashur Höyük and Salat Höyük and not higher than 2.5% in Kavus ¸an Höyük.Type C's fingerprint was not needed by the model to explain the composition of bitumen from sites other than TL.Only Umm an Namel presented a lower than 2.5% chance of it contributing more than 20% to a mixture with other sources.Finally, hopane ratios coherent with types A and B were reported in all sites both north and south from TL.Moreover, the model was not able to find a mixture of TL's types which could result in the values obtained either in Bashur Höyük or from sample 993a in Tell F6.This could indicate that these sites had access to trade routes transporting bitumen types which did not reach TL.Alternatively, they may have chemically altered bitumen in a way that made it dissimilar to that from TL. Finally, this types may also be unrepresented in our sample collection.

Exploring the geographic origin of archaeological bitumen
Given the absence of oleanane in our samples, it was improbable that TL bitumen originated from Iran or Azerbaijan (Connan and Deschesne, 1995;Connan et al., 1998;Bechtel et al., 2014).Furthermore, an origin in Oman, where modern references present δ 13 C Asphaltene values around − 34‰ (Connan et al., 2005) or Jordan, where δ 13 C Asphaltene values are placed at − 29‰ (Connan et al., 2006a), would be equally unlikely.Moreover, all reference values from sources in modern Syria and the Black sea report 30G/30αβ ratios above 0.5 (Connan and Nishiaki, 2003;Hauck et al., 2013;Fulcher et al., 2020), which none of TL's archaeological samples present.Other sources in Italy (Pennetta et al., 2020) have not been considered due to their absolute distance from

Table 3
The archaeological context and the results of the GC-MS and EA-IRMS analyses with their associated standard deviation.north-eastern Mesopotamia.Alternatively, modern sources across southern, central, northeastern, and north-western Mesopotamia (Fig. 1) presented values resembling those in samples from TL.At this point, the third Bayesian mixing model was used to better distinguish only between the aforementioned zones.The results (Table 5) indicate that bitumen from three different regions might have reached TL.
In types A (samples 5, 6, 7, 10) and D (samples 11 and 12), the probability that most of the bitumen in the residue originated from either north-eastern, north-western or southern Mesopotamia was always below 10%.Alternatively, the likelihood their origin was from

Table 4
Mean and standard deviation of the quantity of each bitumen type present in samples from TL according with the results of the first Bayesian mixing model.central Mesopotamia was always higher than 50% and the LDA analysis assigned them to this region.
In type B (samples 1 and 8), the highest probability was that most of the bitumen originated from north-eastern Mesopotamia (49% and 42% respectively).The LDA analysis assigned them to this region.Coherently, the model allocated significantly low probabilities to an origin from central, southern, or north-western Mesopotamia.
Type C (Sample 9) was placed in a separate category despite the high probability that this bitumen originated from central Mesopotamia (63%), an assumption supported by the LDA analysis.This is because the model did not fully reject the possibility that more than 50% of the residue originated from north-western Mesopotamia (24%).Given that no modern references presented similar biomarker values, it is likely that this sample is the outcome of a mixture of bitumen sourced from these two regions.

Strengths and caveats of Bayesian mixing models
Across our research into bitumen trade, Bayesian mixing models have provided us the capacity to treat archaeological data as probability distributions reflecting both the analytical results and their associated uncertainties.By assuming all samples are the product of a mixture, the model evaluated the full range of combinations which could explain the obtained values, including the possibility that they not resulted from a mixture in the first place.In our approach, modern oil seeps were grouped in regions which display coherent biomarker and isotopic ratios.Consequently, the Bayesian model would be able to approximate the location of unknown exploited seeps originating from the same crude oil reservoirs as other nearby seeps included in the modern reference set.
However, these models may conduce to erroneous interpretations when used as a black box.In our case, the complex composition of petrogenic organic matter implied that at least eight hopane and stearane ratios were necessary to distinguish the four bitumen types present

Table 5
The probability each region was the major source (>50%) of each bitumen sample from Tell Lashkir according to the results of the third Bayesian mixing model and its interpretation after the LDA analysis.in TL.When this was not possible in the second model, a clear separation between types A and B could not be achieved.Therefore, in cases where the model is most prone to error, results should be always carefully evaluated both across modern references and published studies from other archaeological sites.This is the case of samples 1 and 8, tentatively interpreted as originating from north-eastern Mesopotamia despite the model's bias towards producing false positives from this region.The archaeological reference samples from southern Mesopotamia which were misidentified by the model as from the North-East had Ts/Tm ratios exceeding 0.23 and δ 13 C Asphaltene values above − 27.5‰.This was never the case in samples 1 and 8, where the Ts/Tm ratios were lower and the δ 13 C Asphaltenes were placed below − 27.5‰.Similarly, all the misidentified archaeological reference samples from central Mesopotamia, which presented a δ 13 C Asphaltenes value like that of TL's bitumen, had different Ts/Tm ratios.While we can posit that it is safe to assume that samples 1 and 8 do not constitute false positives, the limitations of the model did require that this possibility was evaluated.
Overall, four major considerations should be taken into account when applying Bayesian mixing models to the apportionment of archaeological bitumen mixtures.First, the model's discrimination power rises significantly when higher numbers of hopane and sterane ratios are used as proxies.In our case, the first model was capable of distinguishing two different fingerprints (A and D) which were both assigned to the central Mesopotamian region by the third model.The second model could not distinguish between types A and B and, additionally, the uncertainties of the first model posteriors were significantly lower than those of the second and third models.This suggests that there is room for significant improvements in precision and accuracy as more hopane and sterane ratios are made available.
Second, despite bitumen sources mentioned by cuneiform tablets (Heimpel 2009) and the variability detected in archaeological samples are used to infer good coverage of potential past sources by the modern reference set, still, cases where archaeological hopane and sterane biomarkers do not match (ex: Connan 1999; Van de Velde 2017b; Connan and Oates 2018) point at potential gaps in the reference data.Partially mitigating this problem, Bayesian mixing models offer the opportunity to test whether unmatched samples could be the result of a mixture from known seeps, as it was the case with sample 9 in TL, but won't be able to compute cases where bitumen originated from a crude oil reservoir missing from the reference dataset.Alternatively, the use of intra-site and inter-site archaeological references will reflect the structure of bitumen trade in Bronze Age Mesopotamia regardless of the exact location of each seep, thus avoiding this issue completely.
Third, the Bayesian mixing models were most effective when deployed as ad hoc solutions supported by prior knowledge.A fair understanding of the most likely sources and a careful evaluation of the model's sensitivity were key to develop models capable to answer our archaeological questions.In the case of the third model, reducing the number of possible sources to four was key to ensure model convergence.Therefore, it should not be ruled out that the range of types detected in TL could be reformulated as additional samples are analysed.
Fourth, given that the model outputs are statistical distributions, they can be used to test specific archaeological questions.This makes possible to claim that the probability that bitumen extracted in the Persian Gulf was a major component in the analysed samples is lower than 6%.This can also be applied to test the similarity of different samples or groups of samples.In our case, the probability that a sample of bitumen from the Bronze Age layers contained major amounts (>50%) of the same type than a sample from the MBA layers is 11%.Finally, a posteriori addition of sample results from similar periods or archaeological contexts can be used to significantly improve the accuracy of the results (see error margins in Figs. 6 and 7).

Historic significance of TL's bitumen remains
Bearing in mind the relatively low sample size, our results so far suggest that bitumen from sectors 3 and 6 might have contained different geochemical fingerprints from sector 1 (Fig. 5).In the Ninevite 5 period, type D would have been the most abundant (38-47%), closely followed by B (23-35%) and A (22-34%).Later, in the MBA, the major type could have shifted to A (43-62%), closely followed by C (33-39%), and the presence of B (1-23%) and D (1-12%) would have been significantly reduced.Nonetheless, Type's A ubiquity in the three studied sectors suggests that the supply of this source of bitumen could have been somewhat stable across time.Major changes in the origin of the bitumen reaching TL would not be unexpected as sectors 3 and 6 A. Breu et al. predate the formation of the Assyrian imperial core (Ashur-Niniveh-Erbil) while sector 1 would be placed in a period when the first empires had started the integration of territories across Mesopotamia and beyond.
Neighbouring contemporary archaeological sites including Qalinj Agha (Abu al-Soof, 1966;1969), Tell Nader (Beuger 2016: 20), Erbil's citadel (Nováček et al., 2008), and later sites such as Kurd Qaburstan (Schwartz et al., 2017), the Darband-i Rania pass (MacGinnis et al., 2020), and previous landmarks in the region such as Shanidar Cave (Solecki et al., 2004), Gird Banahilk (Gómez-Bach et al. 2019) or Jarmo (Braidwood and Braidwood 1950) present considerable amounts of bitumen remains.This suggests that, in northern Mesopotamia, bitumen had been and was acquired, consumed, and traded by most human settlements in the region.Resembling the processes in motion in southern Mesopotamia, this was a product deeply intertwined in the Bronze age way of life.In TL, the detection of obsidian blades, moulds for metallurgical production and a diverse range of pottery types (Molist et al., 2019a) are coherent with the site's nature as a small settlement operating as logistic point through long trade routes crossing Mesopotamia.In this sense, TL would not be unique, as other archaeological sites including Bardastee (Eidem, 2015), Bash Tapa (Marti et al., 2015) or Kani Shaie (Tomé et al., 2015) exhibit similar characteristics (Molist et al., 2019a).Nonetheless, TL has offered us the opportunity to probe the range of bitumen types in circulation across the region.Results from our second mixing model indicate that TL's occupants were involved in Mesopotamian-wide bitumen trade networks, as exemplified by types A and B. In the early third millennium, they were able to acquire products which also reached the Persian Gulf (type D bitumen) but were not present upriver (Kavus ¸an Höyük, Salat Höyük and Bas ¸ur Höyük, Connan et al., 2013;2018).Even though the sample size is still low, the absence of oleanone, and hopane, sterane and asphaltene values from modern Azerbaijan or Iran should not be underestimated.This set of results foreshadows the tighter relationship this region would develop with southern Mesopotamia when, in the third millennium, the first written record of the city of Erbil appeared in an inscription from Nippur (MacGinnis 2014).
During the Early Dinastic Period, tablets from Girshu, Umma and Gudea indicate that expeditions consisting of up to one hundred men were sent to obtain copious quantities of bitumen upriver to either northern or central Mesopotamian seeps (Heimpel 2009).The analysis of bitumen remains from Susa suggests that, across time, the kingdom of Elam acquired bitumen from closer sources in the Susiana plain (Connan and Deschesne 1996).At TL, the outcomes of the third model indicate that central Mesopotamian seeps (types A and D) might have been contributing bitumen to the Tigris trade networks significantly more than their north-eastern Mesopotamian counterparts.Nonetheless, results from samples 1 and 8 suggest that settlements in the upper Tigris might still have played a significant role as producers of this raw material.In the Euphrates, analyses from Tell Brak (Connan and Oates 2018) and Mari (Connan and Deschesne 2007) suggest bitumen was also acquired from a combination of northern and central Mesopotamian seeps.
It is still difficult to assess the specific role bitumen played in TL.So far, their absence in architecture and their overwhelming presence as a coating in pottery suggests it might have been commonly used as an adhesive or as a waterproofing agent in pottery.The absence of petrographic studies prevents assessing whether the provenance of its ceramic containers could be tied with the trade of this product.Nonetheless, it is our interpretation that, given that the Bronze age buildings in sector 1 are coherent with storage spaces, it should not be ruled out that part of the detected bitumen was not used, but stored in large jars waiting for the right trading conditions.Rather than importing significant quantities of bitumen from a single source, the Bronze Age inhabitants of TL acquired it from several regions, a strategy which could have mitigated the risk of losing access to this raw material.Moreover, the possible detection of mixtures (samples 1, 5 and 7) suggests that bitumen could have been combined and recycled when necessary, a practice which was not uncommon (Schwartz and Hollander, 2000).

Conclusions
The novel implementation of Bayesian mixing models to the study of archaeological bitumen has been successful in facilitating and improving the apportionment of mixtures from various sources.The first analyses of bitumen from north-eastern Mesopotamia in the early third millennium (Ninevite 5 horizon) suggest that the population inhabiting Tell Lashkir owned products found in exchange networks which reached southern Mesopotamia.We propose that Tell Lashkir's inhabitants traded with bitumen originating from both nearby sources and distant regions including central Mesopotamia.While results support the widespread use of central Mesopotamian bitumen in this period, the putative origin of samples 1 and 8 is coherent with written accounts and geochemical analyses indicating that bitumen from the north-east also reached southern cities (Connan and Van de Velde, 2010).
Bitumen was instrumental to the origin of the first civilizations and states, as well as the development of the first global trade networks.As this paper has shown, the use of more refined statistical techniques to interpret chemical data on archaeological bitumen allows deciphering in exceptional detail the apportionment of ancient bitumen trade in southwestern Asia.Thus, the combined biomarker and Bayesian approach offers new insights into the origin of civilisation and the first states, one of the defining periods in the history of humankind.

Fig. 4 .
Fig. 4. a) Principal Components Analysis.Dashed ellipses group samples according to b. b) Dendrogram presenting possible groupings coherent with the calculation of the Euclidean similarity index.The dashed line cuts the dendrogram at height 0.26.

Fig. 5 .
Fig. 5. Boxplot reporting the normalised total quantity of each bitumen type from the studied assemblage in Tell Lashkir.

Fig. 6 .
Fig. 6.Boxplots present the percentage of bitumen by weight predicted by the first Bayesian mixing model in samples most probably containing a mixture of two types.

Fig. 7 .
Fig. 7. Boxplots present the percentage of bitumen by weight in a) the combination of six samples from the early 3r millennium cal BC b) the combination of three samples from the Middle Bronze Age.

Table 1
in the Median and standard deviation values representative of each region.
A.Breu et al.