Strontium is a heavy alkaline earth metal with four naturally occurring isotopes: ⁸⁴Sr (~0.56%), ⁸⁶Sr (~9.87%), ⁸⁷Sr (~7.04%) and ⁸⁸Sr (~82.53%). While three of these isotopes are stable, strontium-87 (⁸⁷Sr) is radiogenic and is produced by the radioactive decay of rubidium-87 (⁸⁷Rb) with a half-life of ~4.88 × 10¹⁰ years (Bentley, 2006). Due to this process, the abundance of ⁸⁷Sr in the geosphere is gradually increasing through time, thus altering the ratio of ⁸⁷Sr to the stable isotopes of strontium. The isotopic evolution of strontium in the geosphere began about 4.5 ± 0.1 × 10 (four and a half billion) years ago with a primordial ⁸⁷Sr/⁸⁶Sr ratio of about 0.699 (Faure, 1986). The isotopic composition of strontium in a rock or mineral that contains rubidium is dependent on the age and the initial Rb/Sr ratio of that rock or mineral, and the utility of the isotopic variation of Sr for age determination of bedrock and provenance of minerals and sediments has long been recognized in geology (Nebel et al., 2011). Because of the utility of Rb-Sr decay systems in the estimation of geological age, the variation of Sr concentrations and Rb/Sr ratios in rock types are well studied in many parts of the world (Capo et al., 1998). As a corollary, Sr isotope ratios of geological formations have also been extensively studied, expressed using the ratio of ⁸⁷Sr relative to ⁸⁶Sr corrected for any mass-dependent fractionation by normalization to a fixed ⁸⁶Sr/⁸⁸Sr (0.1194) (Bataille et al., 2020). Accordingly, as a general principle, older bedrock typically has a higher ⁸⁷Sr/⁸⁶Sr. Very old rocks (>100 Ma) with original high Rb/Sr have higher ⁸⁷Sr/⁸⁶Sr (>0.710), while recently formed rocks (< 1–10 Ma) with low Rb/Sr ratios have lower ⁸⁷Sr/⁸⁶Sr (< 0.706) (Bentley, 2006). Rocks enriched in alkali metals, alumina, and silica (e.g., clay and shale), which are the oldest rocks on the Earth's crust, have high initial Rb/Sr and tend to have high ⁸⁷Sr/⁸⁶Sr, as well. Likewise, certain granitic rocks with high initial Rb/Sr, which were formed very early during our planet's formation have high ⁸⁷Sr/⁸⁶Sr, although the isotopic ratios of igneous granites may be highly variable. Oceans have a homogenous ⁸⁷Sr/⁸⁶Sr (0.7092) at present, which has varied between ~0.707–0.709 over geological time (McArthur et al., 2001). Carbonates, which have precipitated from seawater, reflect the seawater ⁸⁷Sr/⁸⁶Sr at the time of their formation and approximate the ocean water ⁸⁷Sr/⁸⁶Sr (Farkaš et al., 2025).

A geological map of bedrock types in a region can be used as a preliminary guide to produce a coarse map of expected ⁸⁷Sr/⁸⁶Sr of the basement rock, but the ⁸⁷Sr/⁸⁶Sr in bedrock is rarely the same as the ⁸⁷Sr/⁸⁶Sr in topsoil (Maurer et al., 2012; Bataille et al., 2020). Besides geological age and the initial Rb/Sr, a variety of geological processes related to igneous activity, metamorphism, weathering, and sedimentation play a significant role on the transfer of Sr from the Earth's crust to the surface (White, 2023). These factors influence the isotopic ratios of geological formations and cause their ⁸⁷Sr/⁸⁶Sr to be significantly different than the expected ratios based on Rb-decay. For example, as a general principle, the ⁸⁷Sr/⁸⁶Sr of basalts, which are young volcanic rocks with very low Rb/Sr, can go as low as 0.702; however, because the basaltic magmas thrusted out onto the surface may interact extensively with old sialic rocks, volcanic rocks may exhibit very high ⁸⁷Sr/⁸⁶Sr (Köksal et al., 2004). Research on various volcanic regions around the world has shown that the chemical and isotopic compositions of igneous rocks differ greatly in oceanic vs. continental regions, reflecting a unique set of circumstances and processes (e.g., multiple eruptions in stratovolcanoes, contamination by old sialic rocks, and blending of magmas) in each volcanic province (Prytulak and König, 2025).

Although bedrock (“parent rock”) is the primary source of Sr in soils, many natural processes like weathering, erosion, alluviation, and the mixing effect of rivers contribute to and alter the Sr isotopic ratio of the Earth's surface in each location (Probst et al., 2000). While differential weathering rates of various rocks and minerals result in an uneven contribution from the constituents of parent rock to the Sr budget of surface soils, atmospheric factors like sea-spray, rainfall, and aeolian dust also alter the ⁸⁷Sr/⁸⁶Sr in topsoil (Whipkey et al., 2000; Evans et al., 2010; Hartman and Richards, 2014). Rivers and streams have a high concentration of Sr (~110 ppm), which is derived from both the atmosphere and the weathering of rocks into solution in fluvial basins (Bentley, 2006). There are major factors that cause stream water ⁸⁷Sr/⁸⁶Sr to be different than the underlying bedrock including differential precipitation rates due to climatic and seasonal conditions and differential weathering rates of rocks (Bataille and Bowen, 2012; Bataille et al., 2020) (e.g., rapidly dissolving carbonates and evaporites contribute disproportionately high amounts of Sr). Most importantly, erosion and sedimentation by rivers are crucial for an accurate understanding of ⁸⁷Sr/⁸⁶Sr distribution in surface soils (Bentley, 2006): at higher elevations, there is a stronger correlation between stream water ⁸⁷Sr/⁸⁶Sr and bedrock, because elevated regions erode faster, generate more surface soil, and contribute more to river sediments than the soils at lower elevations. At lower elevations, on the other hand, the correlation between local bedrock and river ⁸⁷Sr/⁸⁶Sr is blurred due to the mixing of sediments along the course of rivers as they flow from higher elevations to the floodplains.

Considering all the processes that contribute to the Sr budget in an ecological zone, an important distinction has emerged between geological ⁸⁷Sr/⁸⁶Sr and bioavailable ⁸⁷Sr/⁸⁶Sr in the utilization of Sr isotope ratio analysis in archaeology in recent years. Accordingly, bioavailable ⁸⁷Sr/⁸⁶Sr is the term that is used for the isotopic ratio of the strontium in the ecosystem that can be absorbed by biological organisms (plants, animals, humans) in their natural habitat (Hartman and Richards, 2014). For an accurate understanding of bioavailable ⁸⁷Sr/⁸⁶Sr, it is necessary to understand the pathways that Sr follows from the Earth to the living organisms (Maurer et al., 2012). Due to its properties as a divalent (+2 valence) element and its similar atomic size, strontium replaces calcium in living organisms (Montgomery, 2010). Sr enters the food chain through plants first and it is then transferred into animals and humans who consume plants (and animals). Unlike the fractionation rate of light stable isotopes across the trophic levels of the food chain, which constitutes the main principle for paleodietary reconstruction in bioarchaeology, the radiogenic ⁸⁷Sr does not fractionate across trophic levels (Flockhart et al., 2015; for mass-dependent fractionation in stable isotope ⁸⁸Sr, see Knudson et al., 2010; Lewis et al., 2017). Plants absorb the atmospheric Sr through their leaves and the Sr in soil minerals (sourced from bedrock and surface soils) by root uptake, which implies that plants with different root depths absorb Sr from different substrata (Poszwa et al., 2004). The main source of Sr in terrestrial animals is plants, and ⁸⁷Sr/⁸⁶Sr in animal tissues is expected to reflect that of the local plants in their food catchment area but stream and lake water ⁸⁷Sr/⁸⁶Sr may contribute and skew (especially in carbonate-rich environments) the ⁸⁷Sr/⁸⁶Sr in animal tissues (Lewis et al., 2017; Bataille et al., 2020). ⁸⁷Sr/⁸⁶Sr in human tissues is also sourced from dietary intake (including plants, animals, and water), and in theory the ⁸⁷Sr/⁸⁶Sr in various parts of the human skeleton that develop over the lifetime of an individual reflects the averaged ⁸⁷Sr/⁸⁶Sr of the constituents of human diet at the time of the development of the sampled skeletal/dental tissue, which is the foundational premise for the utility of ⁸⁷Sr/⁸⁶Sr as a mobility tracer in archaeology (Montgomery, 2010).

The use of bioavailable ⁸⁷Sr/⁸⁶Sr as a mobility tracer in archeological animal and human populations is founded upon two theoretical principles. As reviewed above, the first principle is that the primary source of ⁸⁷Sr/⁸⁶Sr for a biological organism is the bioavailable ⁸⁷Sr/⁸⁶Sr in its food catchment area (see reviews in Price et al., 2002; Bentley, 2006). Accordingly, ⁸⁷Sr/⁸⁶Sr in the tissues of sedentary fauna that have a local feeding range (e.g., rodents and other small mammals) incorporate the bioavailable ⁸⁷Sr/⁸⁶Sr of the locality where they feed. Larger mammals and migratory taxa with a wider feeding range, on the other hand, incorporate ⁸⁷Sr/⁸⁶Sr from multiple localities (Toncala et al., 2020). For domesticated species, socio-economic factors and cultural choices further complicate the picture as animals may be grazed intentionally at different localities, varying altitudes, and may even be traded across long distances (Hammer and Arbuckle, 2017). In the case of human populations from hunter-gatherers to complex societies, multiple cultural and economic factors contribute to dietary intake including the proportion of plants and domesticated and hunted animals in diet, imported foods, and locations of available freshwater resources (for an Sr mixing model in human diet, see Montgomery et al., 2007). Therefore, in human mobility studies, it is among best practices to sample a variety of archaeological plants and animals, as well as stream water and groundwater, to assess the numerical range of ⁸⁷Sr/⁸⁶Sr for the local food catchment area (“local range”) as part of baseline determination (see discussion in Lengfelder et al., 2019; Snoeck et al., 2020; Frank et al., 2021a).

The second principle in mobility studies utilizing ⁸⁷Sr/⁸⁶Sr hinges on the chronological development of skeletal tissues and dental enamel across the lifetime of animal and human individuals (Müller et al., 2025). In humans, enamel formation of permanent molars begins in gestation stage and is completed in adolescence (Hillson, 1996; AlQahtani et al., 2010). Because the crown of each molar completes formation in about 3-4 years in succession and enamel does not undergo turnover in later life (unlike bone), the enamel of each tooth preserves a record of the ⁸⁷Sr/⁸⁶Sr at the geographical location where the human individual was present during the crown formation of that specific tooth (Montgomery, 2010). With this understanding, a series of enamel samples from the first, second, and third molars of a human individual can reveal their lifetime mobility from the time of birth to the beginning of adulthood. When the ⁸⁷Sr/⁸⁶Sr of enamel samples are compared with the “local range” ⁸⁷Sr/⁸⁶Sr at the burial location of the sampled individual and/or when compared with bone samples from the same individual, which reveal an averaged ⁸⁷Sr/⁸⁶Sr of their food catchment area in the last decade of life before death, migrant and mobile individuals can be identified with confidence (e.g., Haverkort et al., 2008; Giblin et al., 2013). In animal management studies, incremental enamel samples from the tall crowns of herbivores yield the ⁸⁷Sr/⁸⁶Sr at the location where the animal was feeding during the period the sampled teeth were developing, thus revealing seasonal grazing patterns (e.g., Balasse et al., 2002; Chazin et al., 2019) or imported animals (Arnold et al., 2016). In any case, the key to a correct interpretation of ⁸⁷Sr/⁸⁶Sr data and archaeological/historical inferences drawn from isotopic analysis is a robust baseline map of ⁸⁷Sr/⁸⁶Sr signatures for the geographies pertinent to the research questions being pursued.

The first approach, kriging, is a geostatistical algorithm that uses a limited set of empirical data from sampled regions (“cells” or polygons) to interpolate the value of a variable over a continuous spatial field within a geographically constrained area (e.g., Willmes et al., 2018; Lugli et al., 2022; Tarrant et al., 2024). Where datasets are geographically scattered, a continuous surface of predicted values can be created with this method, i.e., geographical areas where no empirical data exists for the bioavailable ⁸⁷Sr/⁸⁶Sr can be assigned estimated ⁸⁷Sr/⁸⁶Sr values. The presence of high-resolution geological ⁸⁷Sr/⁸⁶Sr values and the density of empirical bioavailable ⁸⁷Sr/⁸⁶Sr data points increase the accuracy of the kriging method in a study region. Kriging utilizes either deterministic interpolation, in which the similarity of nearby cells is determined by a mathematical formula, or geostatistical interpolation, in which this formula is enhanced by an element of probability built into the model (Wheatley and Gillings, 2002). The most common method in archaeological studies is Inverse Distance Weight (IDW) method, a deterministic interpolation. This method was effectively used to create an ⁸⁷Sr/⁸⁶Sr isoscape of southwest Anatolia based on 283 shallow-, medium-, and long-rooted plant and snail-shell samples from 87 locations (Wong et al., 2021; see Table 2). However, as observed by Bataille et al. (2018) in their work on West European ⁸⁷Sr/⁸⁶Sr isoscapes, “the accuracy of these empirical interpolation methods is generally low due to the non-normal distribution of ⁸⁷Sr/⁸⁶Sr data and the non-continuous scalar patterning of ⁸⁷Sr/⁸⁶Sr variability”. In other words, when modeling predicted isoscapes of strontium for the purposes of bioarchaeological investigations of mobility, it is necessary to account for the multiple variables that cause the bioavailable ⁸⁷Sr/⁸⁶Sr to diverge from the geological ⁸⁷Sr/⁸⁶Sr.

To overcome the lack of quantitative understanding of controls in bioavailable ⁸⁷Sr/⁸⁶Sr variability, Bataille et al. (2018) developed a machine learning algorithm combining process-based models and regression techniques to predictively model ⁸⁷Sr/⁸⁶Sr isoscapes at various scales. For the basis of this model, interpolation of empirical ⁸⁷Sr/⁸⁶Sr data (as explained above) was combined with a mechanistic bedrock ⁸⁷Sr/⁸⁶Sr model, which hinges on the radiogenic equation accounting for the decay of ⁸⁷Rb into ⁸⁷Sr. Distinct from empirical models that rely on deterministic or geostatistical formulas, however, multiple covariates that impact bioavailable ⁸⁷Sr/⁸⁶Sr variability were incorporated into the model and the predictive model was trained by testing various regression techniques to evaluate and optimize its predictions. In the developing stages of the model (Bataille et al., 2018), researchers first computed the covariates “parent rock, aerosols, biological processes, relief, climate, and water dynamics” (following Capo et al., 1998), and the later version of the model was improved by inclusion of additional covariates like “global nitrogen and phosphorus fertilization, surficial deposits, global mean annual temperature, and an updated raster of global sea salt aerosol deposition” (Bataille et al., 2020). By testing the model utilizing multiple regression techniques, it was noted that random forest regression (RFR) yielded results with the highest precision (see below for details). The R-script for this global ⁸⁷Sr/⁸⁶Sr isoscape is made available open-access and it can be regenerated and updated by uploading new data, which is a great service to many disciplines including our research community. Machine learning approaches to bioavailable Sr isoscape reconstruction in world regions where proxy baseline data is scarce is rapidly developing (e.g., Janzen et al., 2020; Barberena et al., 2021; Gigante et al., 2023).

It is of utmost importance for the purposes of this study to emphasize that the global model performs best in world regions where abundant ⁸⁷Sr/⁸⁶Sr data from various substrates is available (e.g., Europe), while in regions lacking bioavailable ⁸⁷Sr/⁸⁶Sr data (e.g., Madagascar) predictions do not characterize the isotopic heterogeneity accurately (Bataille et al., 2020). Because Türkiye is among relatively data-poor regions for bioavailable ⁸⁷Sr/⁸⁶Sr, the global model would be expected to perform with low accuracy in this region relative to other world regions. In fact, the previously published global model did not include any datapoints from Türkiye. Including empirical ⁸⁷Sr/⁸⁶Sr data from Anatolian sites into the global model would optimize the regional model for Türkiye. Furthermore, the inclusion of this dataset could also aid in increasing the resolution of the model at a global scale. Before focusing on the state of strontium isotope ratio research in the archaeology of Türkiye and the compilation of bioavailable ⁸⁷Sr/⁸⁶Sr data, we will briefly review available isotopic data compilation initiatives in neighboring regions to contextualize this effort.

With the rise of open data movement, scientific communities around the world have begun to embrace FAIR data principles (Findability, Accessibility, Interoperability, and Reusability, after Wilkinson et al., 2016). In addition to the long-standing open-access database platforms in (paleo)environmental research (e.g., Neotoma) and radiocarbon dating (e.g., CARD), global databases of stable isotope ratios (e.g., waterisotopes.org, IsoBank, Faunal Isotopes Database) and open-access platforms hosting regional isotope ratio repositories (e.g., IsoMemo network via Pandora) have developed in the last 20 years. In the last decade, there has been an upswing in the development of regional, open-access databases for bioarchaeological isotope ratio data in the circum-Mediterranean, as well. The IsoArcH database, which incorporates bioarchaeological isotopic data from the Greek and Roman world (Salesse et al., 2018), AfriArch repository that integrates bioarchaeological isotopic measurements from African archaeological sites (Goldstein et al., 2022), IsoMedIta: A stable isotope database for medieval Italy (Mantile et al., 2023), and SrIsoMed: An open access strontium isotopes database for the Mediterranean (Nikita et al., 2022) are notable recent initiatives in the regions just west and south of Anatolia (Türkiye). However, the “Greater Ancient Near East” from the Balkans to the south Caucasus, and from northern Sudan to the Persian Gulf and Iran, remains on the margin of these geographies covered by the compilation of isotopic data (also see Nafplioti, 2011, 2021; Frank et al., 2021a,b). This is mainly due to the historical divides between disciplinary fields of European, Graeco-Roman and Near Eastern archaeology, language barriers that set research communities apart, and the much-delayed arrival of isotopic analyses to the archaeology of southwest Asia.

To address this lacuna and to move the field of isotopic analysis forward in the greater ancient Near East, the project “BioIsoANE: An open-access repository of bioarchaeological isotopic analysis in the greater ancient Near East” was initiated in recent years (Irvine and Yazıcıoğlu Santamaria, 2021). Once completed, the BioIsoANE website will host a comprehensive bibliography of published bioarchaeological isotopic analyses and a database of all published bioarchaeological isotopic data georeferenced to an interactive map of the Greater Ancient Near East (as defined above), which can be queried by regions, site names, cultural and chronological periods, sample types, isotopes reported, and research themes. Additionally, the website will contain information about the theory, principles, and methodologies of isotope chemistry, isotopic analyses, and their applications, as well as guidelines for best practices in project development, sample collection and preparation, data analysis, interpretation, and presentation, to serve not only specialized researchers but also students as a pedagogic resource on the utility, power, and limitations of isotopic analysis in bioarchaeology (Yazıcıoğlu and Irvine, 2025). The review of bioarchaeological ⁸⁷Sr/⁸⁶Sr analysis conducted to date in Türkiye, presented below, and the compilation of data sourced from published research, which we have incorporated into the global map of ⁸⁷Sr/⁸⁶Sr isoscapes in this study, are made possible by the efforts in compiling the BioIsoANE database.

⁸⁷Sr/⁸⁶Sr analysis in Türkiye only began in the 2010s and has progressed slowly but steadily in recent years. Table 1 presents (A) a list of seventeen archaeological sites (ordered chronologically from distant to recent past), where ⁸⁷Sr/⁸⁶Sr analysis was conducted to investigate various research questions, followed by (B) a list of additional five sites where bioarchaeological ⁸⁷Sr/⁸⁶Sr analysis is ongoing. Among sites with published data, research themes pursued by ⁸⁷Sr/⁸⁶Sr analysis include human mobility (eleven sites), animal management and herding practices (seven sites), management of aquatic resources (one), timber provenance in shipwrecks (one), and glass provenance/recycling (one). Research into ⁸⁷Sr/⁸⁶Sr isoscapes is a relatively novel field in the region, pursued by the authors of this paper in central Anatolia and the study by Wong et al. (mentioned above) in southwestern Anatolia. Pioneering work in isotopic analysis in Anatolian archaeology is to be credited to excavation projects with a vision of long-term, multi-disciplinary landscape archaeology like Neolithic Çatalhöyük in central Anatolia and Classical period Sagalassos in the southwest (Figure 1). Early prehistoric sites of southeastern Anatolia and south-central Anatolia constitute hot spots for multi-isotopic analyses, including strontium. The impetus for the use of ⁸⁷Sr/⁸⁶Sr analysis in early prehistory is a fine-resolution understanding of human-environment interactions across the transition to sedentism and the incipient and advanced stages of domestication of plants and animals.

Beginning with the Bronze Ages, the rationale for the use of ⁸⁷Sr/⁸⁶Sr analysis shifts from human-environment interaction to the investigation of long-distance human mobility to shed light on the multi-ethnic composition of protohistoric Anatolian communities. Notably, in each case, ⁸⁷Sr/⁸⁶Sr data has revealed that human mobility in proto-urban and urban societies are far more complex than what is implied by the textual record. For example, at the coastal EBA cemetery of İkiztepe, where the site's excavators had hypothesized long-distance migrations via the Black Sea to explain the early dispersal of proto-Indo-European speakers into Anatolia, ⁸⁷Sr/⁸⁶Sr data revealed individuals engaged with transhumant pastoralism in addition to possible emigrants from distant regions (Welton, 2015). At Kültepe (ancient Kanesh) in south-central Anatolia, which is renowned for its large corpus of cuneiform texts belonging to an enclave of resident merchants from the city of Assur in northern Iraq, ⁸⁷Sr/⁸⁶Sr data helped identify immigrants from multiple other locations including north-central Anatolia and possibly the upper Euphrates or other Bronze Age cities in south-central Anatolia (Yazıcıoğlu Santamaria, 2017). At another cosmopolitan MBA-LBA city, Tell Atchana (ancient Alalakh), combined ⁸⁷Sr/⁸⁶Sr and aDNA analysis shed particularly interesting light on the period of societal collapse between the levels characterized by palaces with administrative archives (Ingman et al., 2021). And finally at Hierapolis and Ephesos, ⁸⁷Sr/⁸⁶Sr was used for identifying pilgrims (potentially from Europe) in the early Christian era (Wong et al., 2018). These examples demonstrate that ⁸⁷Sr/⁸⁶Sr analysis of human remains is a particularly relevant methodology in text-aided archaeology because it provides contextualized empirical data that help us conceptualize human mobility and ethnicity beyond the straightjacketing culture-area paradigm and “the tyranny of the text”.

The map of Anatolia in Figure 1 shows that the geographical distribution of sites with analyzed remains is far from even. The number of archaeological samples and baseline reference samples analyzed in each study according to sample type are summarized in Table 2, which demonstrates this uneven distribution across sites and regions in further detail. For example, currently there are tooth enamel samples from 495 human individuals buried at Anatolian sites across the periods, and of this total, 76 are from southwestern Anatolia, 124 are from south-central Anatolia, and 223 are from southeastern Anatolia, while only 72 samples are from north-central Anatolia. To a certain extent, this trend reflects the geographical distribution of archaeological sites with well-preserved and curated human remains unearthed by modern excavation methods. Archaeological fauna analyzed for documenting animal management and herding practices show a somewhat similar distribution: of the total 302 sampled terrestrial animals, 231 are from sites in southeastern Anatolia, 18 are from south-central Anatolia and 53 are from north-central Anatolia, while in southwestern Anatolia faunal samples are represented only by fish (14 samples). Analyzed archaeobotanical samples are minimally represented (5 Neolithic plant samples and 4 LBA wood samples). Reference samples analyzed for ⁸⁷Sr/⁸⁶Sr baseline determination in each study are diverse (Table 2). Sample types include archaeological and modern plants, small fauna, snail shells, fish, water, soil, rock, salt, and dung. In terms of sample type distribution, modern plants constitute the largest dataset (298 samples) with an additional 16 samples from archaeological plants, followed by 64 samples from modern and archaeological small fauna, 62 samples from modern and archaeological snail shells, 60 soil/rock samples, 30 modern and archaeological fish samples, and 13 water samples (Supplementary material 2).

To reiterate, while the temporal span of ⁸⁷Sr/⁸⁶Sr data is quite impressive, there remain significant geographical gaps in the northern half of the country. For the purposes of this study, the uneven distribution of analyzed baseline samples is significant because the number of georeferenced empirical datapoints will add to the weighted probability in the isoscape interpolation and it will influence the accuracy of the resultant ⁸⁷Sr/⁸⁶Sr isoscape map.

In archaeological mobility studies utilizing ⁸⁷Sr/⁸⁶Sr analysis at prehistoric Anatolian sites (e.g., Göbekli, Nevali Çori, Çatalhöyük, Domuztepe, İkiztepe), researchers have typically used small and sedentary archaeological fauna (rodents, fox, pig) from the researched site as proxy samples for the determination of the site's local food catchment area, while in more recent studies (e.g., Körtik Tepe), a broader local catchment area has been characterized by water and modern plant samples (Table 2). Because land-snail shells are relatively easy to collect and have a high concentration of Sr, this sample type has also been frequently used in research typically at Bronze Age and Classical period sites (e.g., Kültepe, Tell Atchana, Ortaköy, Hierapolis, Ephesos), where the food catchment area is expected to be wider than the immediate vicinity due to food provisioning strategies of state societies. At Tell Atchana and Kültepe, for example, baseline determination for the site's catchment area was carried out at the scale of the settlement plain and sampling spots were selected by targeting geologically distinct units using bedrock geology maps as reference.¹ At the highland site of Ortaköy, where the research question was concerned with the use of pasturelands for herding strategies, snail shells and rock samples were collected from distinct geological units for baseline determination within a much broader radius around the site (Pişkin et al., 2023). Baseline determination for the human mobility study at Kültepe was extended to a regional scale in accordance with its research questions concerning inter-site mobility between MBA cities involved in a trade network. To this end, snail shells and plant samples were collected with regular intervals along the watersheds of two major rivers (Kızılırma and Yeşilirmak) in central Anatolia with the objective of understanding the mixing effect of these rivers. Analysis of samples nearby known Bronze Age sites were prioritized given budget constraints (Yazıcıoğlu Santamaria, 2017); samples analyzed in later years are included in this study (Supplementary material 2) and analyses are ongoing.

The first systematic study for creating a local isoscape, which enables the ⁸⁷Sr/⁸⁶Sr characterization of a continuous surface rather than scattered spots on the map, was carried out around the prehistoric site of Köşk Höyük in south-central Anatolia (Meiggs et al., 2018). Modern grasses were collected from sampling areas that were chosen not only based on bedrock composition but also different geomorphological contexts where alluvial and colluvial processes contribute to the mixing of sediments. Using kernel smoothing and spline features in ArcGIS, a local ⁸⁷Sr/⁸⁶Sr isoscape with contour lines was interpolated, which was then used for mapping the grazing grounds of Chalcolithic sheep/goats based on results from incrementally sampled enamel. In recent years, a larger-scale systematic study was conducted for regional baseline determination in southwestern Anatolia, in which modern snail shells and short-, medium-, and long-rooted plant samples were collected from each sampling spot using a hand-held GPS, and an estimated ⁸⁷Sr/⁸⁶Sr isoscape was created using the inverse distance weight method (Wong et al., 2021), as described above.

Reviewing the empirical data layer included with the publication of the R-script made available open-access by Bataille et al., we realized that ⁸⁷Sr/⁸⁶Sr data generated by archaeological research projects conducted in Türkiye were not present in this dataset (Figure 2). Likewise, bioarchaeological and baseline ⁸⁷Sr/⁸⁶Sr data from Anatolian sites were largely absent from the SrIsoMed database online, in which predominantly geological ⁸⁷Sr/⁸⁶Sr data is compiled from Türkiye.² We compiled all available publications of archaeological mobility and provenance studies and dissertations utilizing ⁸⁷Sr/⁸⁶Sr in Anatolia from the BioIsoANE database (see above) and added unpublished baseline data from the authors' ongoing analyses. First, all ⁸⁷Sr/⁸⁶Sr data from archaeological and baseline samples were compiled into an excel sheet for comprehensive coverage (summarized in Table 2), and then this data was sorted and screened according to sample type to assess their relevance and applicability to the purposes of this study. Datapoints that were reported by authors as not having passed the quality control checks (e.g., ⁸⁸Sr/⁸⁵Rb or ⁸⁴Sr/⁸⁶Sr QC) and datapoints from replicate runs were eliminated. Next, all data was georeferenced by recording latitude and longitude information from original publications. In cases where coordinates were not presented in the original publication, latitude and longitude information was inferred from published maps showing the sampling locations, and accuracy level for the coordinates of each sample was recorded in the datasheet (Supplementary material 2). Data reported only in graphs and with no map for sampling location were not included in this study. After this data cleaning process was completed, further exclusion criteria were introduced to the dataset. Because the purpose of this study is to increase the accuracy of the existing predictive ⁸⁷Sr/⁸⁶Sr isoscape model by incorporating empirical and in situ bioavailable data, a conservative strategy was followed: all human ⁸⁷Sr/⁸⁶Sr data and ⁸⁷Sr/⁸⁶Sr data from all domesticated species that are grazed and herded were excluded from the dataset so as to eliminate any possibility of skewing the resulting map by data that may be influenced by anthropogenic factors. This resulted in an empirical dataset of 688 ⁸⁷Sr/⁸⁶Sr values from archaeological and modern flora, wild fauna, mollusk (snail shell), and soil samples that could be incorporated into the global ⁸⁷Sr/⁸⁶Sr isoscape model.

Prior to generating the random forest model for Türkiye, we assessed the extent to which the existing global model was accurate for Türkiye without integrating additional data, so that we could later evaluate its improvement. To test this, we reproduced Bataille's Sr isoscape for Türkiye in QGIS and populated it with measured sample data from the systematic study by Wong et al. (2021) to calculate an offset between the measured ⁸⁷Sr/⁸⁶Sr values at a given point as compared to the predicted ⁸⁷Sr/⁸⁶Sr value for that same GPS coordinate (Figure 3). We found that agreement was generally good, with plant remains performing slightly better than those for snails (cf. Britton et al., 2020). Both data types were only precise to the fourth decimal point, or ten-thousandths. A higher precision model (to five decimals) would allow for stronger interpretive power, as is desirable for interpretation of mobility in the archaeological record. Considering that no ⁸⁷Sr/⁸⁶Sr data from Türkiye were used in the generation of the original global model, the performance of the model in predicting values was encouraging and suggested that adding empirical data from Türkiye to the model should notably improve the predictive power of the model.

After identifying offsets between predicted and empirical ⁸⁷Sr/⁸⁶Sr measurements in Türkiye based on the general global model, we used the existing R-script and modified the global dataset including 688 new analytical values from baseline samples at and around Anatolian sites to create a new isoscape model for Türkiye (Figure 4). This model used a random forest decision tree based on a priori relationships between variables and multiple iterations to generate a robust model that can then be tested and validated with additional measured data. In order to create a raster of predicted strontium ratios, we began by obtaining geolocated data of Sr isotope analysis samples from the global model (Bataille et al., 2020) in addition to covariate rasters. This spreadsheet included more than 16,000 entries, while the rasters included geologic and climatic data related to the ratio of bioavailable ⁸⁷Sr/⁸⁶Sr. The original script was then trimmed and altered to only include code necessary to generate a raster of predicted ⁸⁷Sr/⁸⁶Sr for the region. After trimming the code, we ran the model several times with adjusted parameters. The final model cropped the extent of the model to the area surrounding Türkiye, before running the random forest regression on a dataset of 8,000+ georeferenced samples. This created a model with a relatively small range of ⁸⁷Sr/⁸⁶Sr values (~0.705–0.719), in which the predicted strontium isotope ratios matched the regional trends in ⁸⁷Sr/⁸⁶Sr identified by previous research. For example, mountainous areas in the northeastern parts of Anatolia had significantly higher ⁸⁷Sr/⁸⁶Sr ratios than other areas, consistent with the findings by Bentley and Knipper (2005) on horst formations of Palaeozoic granites and metamorphic rocks in Europe. After creating our raster of predicted ⁸⁷Sr/⁸⁶Sr ratios for Türkiye and surrounding areas and our accompanying raster of standard errors, we imported the prediction raster into ArcGIS to conduct geospatial analysis and create graphics. We symbolized the raster using a percent clip stretch symbology in order to best display the continuous nature of ⁸⁷Sr/⁸⁶Sr ratios. We then overlaid prominent archaeological sites, cities, and national borders over our raster in order to provide important contextual information to assist us in our later analysis. We constructed buffer zones around our archaeological sites of interest and conducted zonal statistics over these buffers to find the mean predicted ⁸⁷Sr/⁸⁶Sr ratios for each site, so as to compare the means to our established values for these sites to check for model accuracy.

To assess the accuracy of the model within Türkiye, we also computed the predicted standard error across the landscape. This entailed using additional code in order to create a second raster to represent the standard error for each cell of predicted ⁸⁷Sr/⁸⁶Sr. The code paralleled the original random forest regression, but it was adapted to predict standard error instead of ⁸⁷Sr/⁸⁶Sr using the same predictors to best match each cell. On the map presented in Figure 5, colors proceed from cooler to warmer as the standard error progresses from smaller to larger. As we would expect, standard error is relatively smaller in areas represented by a higher number of analyzed reference samples, for example, around Hierapolis, where the systematic isoscape study by Wong et al. (2021) was conducted.