Drinking Locally: A Water 87Sr/86Sr Isoscape for Geolocation of Archeological Samples in the Peruvian Andes

The analysis of 87Sr/86Sr has become a robust tool for identifying non-local individuals at archeological sites. The 87Sr/86Sr in human bioapatite reflects the geological signature of food and water consumed during tissue development. Modeling relationships between 87Sr/86Sr in human environments, food webs, and archeological human tissues is critical for moving from identifying non-locals to determining their likely provenience. In the Andes, obstacles to sample geolocation include overlapping 87Sr/86Sr of distant geographies and a poor understanding of mixed strontium sources in food and drink. Here, water is investigated as a proxy for bioavailable strontium in archeological human skeletal and dental tissues. This study develops a water 87Sr/86Sr isoscape from 262 samples (220 new and 42 published samples), testing the model with published archeological human skeletal 87Sr/86Sr trimmed of probable non-locals. Water 87Sr/86Sr and prediction error between the predicted and measured 87Sr/86Sr for the archeological test set are compared by elevation, underlying geology, and watershed size. Across the Peruvian Andes, water 87Sr/86Sr ranges from 0.7049 to 0.7227 (M = 0.7081, SD = 0.0027). Water 87Sr/86Sr is higher in the highlands, in areas overlying older bedrock, and in larger watersheds, characteristics which are geographically correlated. Spatial outliers identified are from canals, wells, and one stream, suggesting those sources could show non-representative 87Sr/86Sr. The best-fit water 87Sr/86Sr isoscape achieves prediction errors for archeological samples ranging from 0.0017 – 0.0031 (M = 0.0012, n = 493). The water isoscape explains only 7.0% of the variation in archeological skeletal 87Sr/86Sr (R2 = 0.07), but 90.0% of archeological skeleton 87Sr/86Sr fall within the site isoscape prediction ± site prediction standard error. Due to lower sampling density and higher geological variability in the highlands, the water 87Sr/86Sr isoscape is more useful for ruling out geographic origins for lowland dwellers than for highlanders. Baseline studies are especially needed in the highlands and poorly sampled regions. Because the results demonstrate that a geostatistical water model is insufficient for fully predicting human 87Sr/86Sr variation, future work will incorporate additional substrates like plants, fauna, soils, and dust, aiming to eventually generate a regression and process-based mixing model for the probabilistic geolocation of Andean samples.

critical for moving from identifying non-locals to determining their likely provenience. In the Andes, obstacles to sample geolocation include overlapping 87 Sr/ 86 Sr of distant geographies and a poor understanding of mixed strontium sources in food and drink. Here, water is investigated as a proxy for bioavailable strontium in archeological human skeletal and dental tissues. This study develops a water 87 Sr/ 86 Sr isoscape from 262 samples (220 new and 42 published samples), testing the model with published archeological human skeletal 87 Sr/ 86 Sr trimmed of probable non-locals. Water 87 Sr/ 86 Sr and prediction error between the predicted and measured 87 Sr/ 86 Sr for the archeological test set are compared by elevation, underlying geology, and watershed size. Across the Peruvian Andes, water 87 Sr/ 86 Sr ranges from 0.7049 to 0.7227 (M = 0.7081, SD = 0.0027). Water 87 Sr/ 86 Sr is higher in the highlands, in areas overlying older bedrock, and in larger watersheds, characteristics which are geographically correlated. Spatial outliers identified are from canals, wells, and one stream, suggesting those sources could show non-representative 87 Sr/ 86 Sr. The best-fit water 87 Sr/ 86 Sr isoscape achieves prediction errors for archeological samples ranging from 0.0017 -0.0031 (M = 0.0012, n = 493). The water isoscape explains only 7.0% of the variation in archeological skeletal 87 Sr/ 86 Sr (R 2 = 0.07), but 90.0% of archeological skeleton 87 Sr/ 86 Sr fall within the site isoscape prediction ± site prediction standard error. Due to lower sampling density and higher geological variability in the highlands, the water 87 Sr/ 86 Sr isoscape is more useful for ruling out geographic origins for lowland

INTRODUCTION
Over the past few decades, the analysis of radiogenic strontium isotopes ( 87 Sr/ 86 Sr) in skeletal and dental tissues has become a robust tool for identifying probable non-local individuals among archeological samples and for geo-locating them by matching sample 87 Sr/ 86 Sr to landscape 87 Sr/ 86 Sr. While massdependent isotopic fractionation does occur between dietary and water sources of strontium as they are processed by the human body, the effects of fractionation are negligible after laboratory normalization (Graustein, 1989, p. 492;Beard and Johnson, 2000;Knudson et al., 2010, p. 2;Lewis et al., 2017;Bataille et al., 2020). Due to this small magnitude of dietbody offset of strontium isotopes, the 87 Sr/ 86 Sr in animal tissues reflect the geological signature of underlying bedrock, which is incorporated into the mixed sources of food, water, and atmospheric elements in the food chain that are consumed and imbibed throughout the course of tissue development (Graustein, 1989;Bentley, 2006;Grimstead et al., 2017;Bataille et al., 2020). Strontium isotope analysis is typically employed by archeologists to examine questions of past migration and mobility: identifying first-generation migrants in a skeletal sample, constraining areas of possible provenience for those non-locals, and assessing the proportion of geologically local to non-local individuals in skeletal assemblages (e.g., Ericson, 1985;Price et al., 1994Price et al., , 2002Price et al., , 2015Evans et al., 2006;Knudson and Price, 2007;Slovak et al., 2009;Montgomery, 2010;Tung and Knudson, 2011;Frei and Price, 2012;Knudson et al., 2012b). Assuming mostly local foods were consumed, individuals with 87 Sr/ 86 Sr that differs from most other individuals in a sample, or from the measured 87 Sr/ 86 Sr of the local environmental baseline, were likely non-local during the period when the sampled tissue was developing.
Understanding how much difference in 87 Sr/ 86 Sr across the sample is required for a non-local designation and then reconstructing the likely provenience of those suspected non-locals depends on accurately characterizing the 87 Sr/ 86 Sr expected throughout the region(s) where they might have lived or traveled. Methods for defining the expected local 87 Sr/ 86 Sr range have shifted over the last few decades from statistically parsing human skeletal 87 Sr/ 86 Sr around the sample mean (e.g., Wright, 2005;Slovak et al., 2009;Price et al., 2015) to testing archeological and modern fauna with small home ranges (Price et al., 2002Evans and Tatham, 2004;Hedman et al., 2009Hedman et al., , 2018Knudson and Tung, 2011). More recently, many scholars have argued that the only reliable way to determine the local 87 Sr/ 86 Sr range is by analyzing the bioavailable strontium in local environmental reference or baseline materials like soils, plants, and local fauna, and then aggregating those reference 87 Sr/ 86 Sr to geological units, soil units, or statistical groupings of baseline and/or skeletal materials (Valentine et al., 2008;Evans et al., 2010;Maurer et al., 2012;Knudson et al., 2014;Kootker et al., 2016;Grimstead et al., 2017;Willmes et al., 2018;Barberena et al., 2019;Pacheco-Forés et al., 2020;Snoeck et al., 2020). Unfortunately, because underlying bedrock geology is highly variable even at a short distance, and because rock weathers unevenly throughout catchments, these methods do not enable the sourcing of a sample to any location beyond the immediate vicinity of the sampling location. This is in part because the geological units the isotopic catchments are aggregated to are themselves models of the underlying bedrock that do not account for uncertainty error. Furthermore, because mammalian tissues derive their strontium content from the mixing of plants, fauna, water, and atmospheric particles ingested from the diet, sampling any one reference material or overreliance on bedrock geology will inevitably fail to capture the mixing of strontium in the human diet or in the diets of other animal end-members (Bentley, 2006;Grupe et al., 2011;Burton and Price, 2013;Grimstead et al., 2017).
Because the 87 Sr/ 86 Sr of mixed-source isotopic catchments vary predictably across the landscape, models of baseline or reference landscape 87 Sr/ 86 Sr have become essential for reconstructing movements across the landscape throughout human history (Fisher et al., 2010;Frei and Frei, 2011;Laffoon et al., 2012;Knudson et al., 2014;Grupe et al., 2017a,b;Bataille et al., 2018;Willmes et al., 2018;Adams et al., 2019;Hoogewerff et al., 2019;Lengfelder et al., 2019;Scaffidi and Knudson, 2020). Isoscapes are spatially explicit models of isotope variability (Bowen, 2010;Ehleringer et al., 2010;Wunder, 2010;West et al., 2014;Pellegrini et al., 2016). Archeological and anthropological isoscapes have most commonly utilized a geostatistical kriging method to model isotopic variability across the landscape, which relies on underlying spatial auto-correlation in GPS-located isotope measurements to predict isotope values at unsampled locations (Copeland et al., 2016;Pellegrini et al., 2016;Wang et al., 2018;Willmes et al., 2018;Adams et al., 2019;Scaffidi and Knudson, 2020). Geostatistical models have been improved on in recent years by process-based isoscapes that rely on regression or machine learning methods to predict isotope variability as a function of environmental variables Bataille et al., 2014Bataille et al., , 2020Hoogewerff et al., 2019)these are rare in archeological and anthropological applications but are becoming more common Hoogewerff et al., 2019;Zimmer-Dauphinee et al., 2020). In either case, when predictive isoscapes are generated from samples reflecting the bioavailable strontium in the dietary catchment rather than bedrock outcroppings or soils containing weathered bedrock, they have the distinct advantage of capturing and predicting the mix of foods, water, and dust at every geographic location. This is important since studies have generally found that soil and bedrock 87 Sr/ 86 Sr alone are poor predictors for human tissue 87 Sr/ 86 Sr (Price et al., 2002, pp. 120-122;Knudson et al., 2014;Grimstead et al., 2017;Snoeck et al., 2020). Establishing accurate isoscapes of bioavailable 87 Sr/ 86 Sr variation in the natural environment and understanding what proxy materials best reflect human 87 Sr/ 86 Sr are prerequisites for accurately identifying non-locals in archeological assemblages and constraining their likely geographic origins to the most probable locations.
Isoscape modeling is notoriously tricky in mountainous landscapes due to their highly variable geology, hydrology, weather, and climate patterns (Bowen and Revenaugh, 2003;Yamanaka et al., 2015). In the Andes Mountains, seismic activity, droughts, flooding, landslides, and intensive mining and agricultural practices are additional factors that alter 87 Sr/ 86 Sr distributions around the landscape. This is expected to contribute to high heterogeneity of 87 Sr/ 86 Sr in higher-elevation zones. The study region of southern Peru is located in the Central Volcanic Zone of the Andes, where lava flows across a broad geographic area tap the same magma sources, homogenizing surface materials across regions within the volcanic arc (Scott et al., 2018). Similarly, recent findings by Serna et al. (2020, p. 10) confirm that in volcanic zones of Northern Patagonia are influenced away from bedrock 87 Sr/ 86 Sr by the deposition and mixing of "unradiogenic volcanic tephra". . . and "reworked loess mixing Quaternary and Jurassic volcaniclastic sediments." Serna et al. (2020) conclude that dust and sediment reworking processes dominate the geomorphology of bioavailable Sr in some volcanic zones-sampling in these areas is critical for correctly characterizing expected 87 Sr/ 86 Sr.
Another obstacle to accurate geolocation that is particularly pronounced in the Andes is the problem of equifinality (Price et al., 2002)-the fact that 87 Sr/ 86 Sr of geographically distant regions are often identical, so that samples with the same 87 Sr/ 86 Sr, recovered from the same location, may have actually originated in distinct places. For example, the recent metaanalysis of published 87 Sr/ 86 Sr from Andean archeological sites by Scaffidi and Knudson (2020) show that archeological human and animal tissue 87 Sr/ 86 Sr ranges from 0.7038 to 0.7239. The upper and lower ends of the human and animal range are slightly higher than the theoretical Andean range based on bedrock, where 87 Sr/ 86 Sr = 0.7050 to 0.7200 (Knudson et al., 2014). Variation follows an east-west gradient corresponding to increasing elevation, longitude, and distance from the coast, so that locally born individuals from similar longitudes might display identical 87 Sr/ 86 Sr even though they grew up thousands of miles apart latitudinally (Scaffidi and Knudson, 2020). Yet another obstacle is that the mixing of geological strontium sources in Andean food webs prior to their incorporation into human tissues-between soils, crops, fertilizers, animal protein, different water sources, dust, and sea spray-is poorly understood, both in the past and the present.
There is considerable debate over the use of relatively simplistic geostatistical kriging models for predicting 87 Sr/ 86 Sr variability across the landscape. Kriging assumes data are normally distributed, continuous, stationary, and spatially autocorrelated-nearer points should be more similar than farther points (Cressie, 1993;Webster, 2014, 2015). Because bedrock units are discrete blocks, Sr reservoirs are not continuous in 87 Sr/ 86 Sr across the landscape and 87 Sr/ 86 Sr are often nonnormally distributed. Therefore, geostatistical methods "can give ambiguous and inaccurate results due to the challenges of incorporating non-normal and skewed distributions" which limit the utility of kriged isoscapes for determining provenience (Bataille et al., 2014(Bataille et al., , 2020. Process and regression-based methods are ultimately necessary for developing high-resolution models of human 87 Sr/ 86 Sr variability (Hoogewerff et al., 2019;Bataille et al., 2020;Serna et al., 2020). Notwithstanding these limitations, kriging has been used to develop coarse, preliminary predictions of bioavailable 87 Sr/ 86 Sr for hydrological (e.g., Varouchakis et al., 2012;Johnson et al., 2018) and anthropological and archeological applications (Frei and Frei, 2011;Copeland et al., 2016;Wang et al., 2018;Willmes et al., 2018;Adams et al., 2019). Published bedrock and volcanic tephra 87 Sr/ 86 Sr are abundant throughout the Andes (Mamani et al., 2008;Scott et al., 2018), but these samples are limited in various parts of the Peruvian Andes, which limits the accuracy of existing bedrock 87 Sr/ 86 Sr models (e.g., Schenk et al., 1999). Therefore, this paper argues that modeling a geostatistical isoscape of bioavailable 87 Sr/ 86 Sr is appropriate as an initial step to characterizing 87 Sr/ 86 Sr variability in these sparsely sampled regions, an effort which will be improved on by subsequent sampling campaigns and future regression-based models.
This paper follows recent ecological and archeological studies to argue that because water averages strontium from throughout the dietary catchment (Frei and Frei, 2011;Frei and Frei, 2013;Zitek et al., 2015;Crowley et al., 2017), imbibed water from the dietary catchment may contribute to a portion of 87 Sr/ 86 Sr variation in archeological skeletal tissues. We report 220 new water 87 Sr/ 86 Sr from GPS-located water samples throughout the Peruvian Andes. We plot these together with 46 published Peruvian water 87 Sr/ 86 Sr (Flusche et al., 2005;Alaica et al., 2020), and generate a geostatistical water 87 Sr/ 86 Sr isoscape from this dataset trimmed of spatial outlier values (n = 249). We then test the validity of the model for geolocating previously compiled human archeological samples considered local to their burial location (from Scaffidi and Knudson, 2020) 1 . Ultimately, a mismatch (beyond the model error range) between the model predictions and a measured archeological sample will enable that individual to be identified as non-local during the time of tissue development and allow archeologists and other researchers in the Andes to constrain possible provenience matches along the modeled 87 Sr/ 86 Sr gradient.

Strontium Isotopes in the Andes: From Dietary Catchments to the Human Skeleton
Strontium enters the body through food and water consumed, and then substitutes for Ca ions in the bioapatite crystals of bone and enamel. While diet-body fractionation throughout the process of digestion and tissue incorporation is minimal, a trophic level effect called biopurification has been well documented throughout ecological systems. The Sr/Ca is reduced by approximately a factor of five in each trophic level, so that plants contain lower Sr/Ca ratios than soils, but higher Sr/Ca than meat (Comar et al., 1957a,b;Wasserman et al., 1958;Sillen and Kavanagh, 1982;Burton and Wright, 1995;Bentley, 2006). Because strontium follows calcium pathways into skeletal tissues, calcium-enriched foods like salt and seafood can be disproportionately reflected in tissue 87 Sr/ 86 Sr, which could override the signature of the isotopic catchment for individuals with significantly more consumption of salty or marine foods (Sillen and Kavanagh, 1982;Hodell et al., 2004;Slovak et al., 2009;Knudson et al., 2010). Thus, the Sr/Ca ratio and related 87 Sr/ 86 Sr in human body tissues are largely impacted by the proportion and types of foods in the diet. Assuming these foods were mostly locally procured, skeletal 87 Sr/ 86 Sr reflects the geological 87 Sr/ 86 Sr of the dietary catchment of residence during the period of tissue formation (Ericson, 1985;Price et al., 1994Price et al., , 2002Montgomery, 2010).
The 87 Sr/ 86 Sr in the food chain at a given sampling location is a mix of weathered bedrock eroded and mixed into soils, surface and groundwater carrying strontium from different weathered sources, atmospheric sources of strontium like sea spray, ash, and dust, and faunal and plant strontium reservoirs (Bentley, 2006;Bataille et al., 2020). The proportion of these different sources varies across the landscape, and is best understood by localized, multi-source mixing models (Montgomery et al., 2007;Bataille et al., , 2018Bataille et al., , 2020Willmes et al., 2018). In the Andes, 87 Sr/ 86 Sr baseline data from plants, fauna, and water are relatively limited compared to other world regions (but see Barberena et al., 2017Barberena et al., , 2019. Expectations for catchment 87 Sr/ 86 Sr are primarily derived from models based on bedrock outcroppings. Bedrock models illustrate that patterns in 87 Sr/ 86 Sr in the Peruvian Andes are influenced by the age and lithology of the two sub-ranges: the Cordillera Occidental (Western Cordillera), which rises sharply out of the coastal plain into a lower-elevation range, and the Cordillera Oriental (Eastern Cordillera), which is the highest elevation range to the east that gradually descends into the Amazonian basin (Figure 1). Coastal areas throughout the north and south of Peru are comprised by Quaternary sediments with contributions from various alluvial flows from the highlands (Bellido et al., 1956;Schenk et al., 1999), and thus, bioavailable 87 Sr/ 86 Sr is expected to be relatively heterogeneous at the country scale. The near-coastal Cordillera Occidental is comprised of late Cenozoic volcanic rocks, which are older in the southcentral Andes than in the north-central Andes. For example, bedrock outcropping 87 Sr/ 86 Sr in this front range varies from a mean of 0.7043 in Ecuador, to 0.7074 in Arequipa, to 0.7065 in northern Chile (see summary in Knudson et al., 2014, p. 206). In contrast, in the Cordillera Oriental, bedrock is dominated by Paleozoic andesites with more radiogenic (higher) 87 Sr/ 86 Sr, although bedrock has not been well-sampled in this region. In areas between the two ranges, like Lake Titicaca, 87 Sr/ 86 Sr should reflect a mix of alluvial deposits from both ranges, but again, bedrock samples are lacking (see summary in Knudson et al., 2014, pp. 206-207). Finally, the Amazon basin is comprised of Precambrian-Devonian sediments and alluvial contributions from the Cordillera Oriental (Bellido et al., 1956;Schenk et al., 1999). To the far east, Amazonian rivers show a broad range of 87 Sr/ 86 Sr, from 0.7056 to 0.7438, which reflect drainage of highly variable highland geologies in the eastern sub-range (Santos et al., 2015).
At present, it is unclear how bedrock 87 Sr/ 86 Sr variability is integrated into the human dietary catchment in the Andes. Generally, bedrock has broader 87 Sr/ 86 Sr ranges than the catchment-averaged bioavailable 87 Sr/ 86 Sr from biological substrates (Bataille et al., 2020), and of those substrates, omnivorous skeletons show broader ranges than soils, herbivorous fauna, and plants (Valentine et al., 2008;Ladegaard-Pedersen et al., 2020). However, in geologically active areas in the highlands like those around Cusco, Knudson et al. (2014) found wider 87 Sr/ 86 Sr variability in soils than in published bedrock. In coastal locations, soil 87 Sr/ 86 Sr matched bedrock 87 Sr/ 86 Sr but very few coastal bedrock samples were available for comparison. This study shows promise for establishing local bioavailable 87 Sr/ 86 Sr expectations through soil analysis but is limited by its geographic coverage-soils were tested from only 17 locations (n = 114), all but three of which were coastal. Knudson et al. (2014) also found that archeological human and faunal samples were consistently offset from soil 87 Sr/ 86 Sr in many locations. This makes sense since soils and rocks are generally not directly consumed, although rock contributions to the diet might be nominal in groups that process food with grinding stones (Johnson et al., 2019), practice geophagy, or consume geological materials as digestive aids-all were done in the ancient Andes (Johns, 1990). These inconsistencies between bedrock, soil, and skeletal 87 Sr/ 86 Sr could be due to differential weathering of bedrock throughout the steep Andes and/or the overall paucity of bedrock and soil samples from the region; nonetheless, this illustrates the need for a deeper examination of different substrates for the analysis of bioavailable Sr in the region.
Indeed, across the globe, researchers have debated which environmental baseline materials incorporated into the food web are the best proxy for human tissue 87 Sr/ 86 Sr. Baseline studies have accordingly varied in their reference material of focus according to the different landscapes of archeological study locations. Using small modern or archeological fauna to establish a baseline was the earliest approach, which remains one of the most common approaches. Because small animals tend to have small homes ranges, they should average many of the same elements of the isotopic catchment as human end members do (Price et al., 2002;Bentley and Knipper, 2005;Hedman et al., 2009;Kootker et al., 2016; Fernández-Crespo FIGURE 1 | Water sampling locations throughout Peru plotted in decimal degrees (WGS 1984 datum), relative to bedrock geology modeled by Schenk et al. (1999). Major cities and the general location of mountains in the Cordillera Occidental (Western Cordillera) and the Cordillera Oriental (Eastern Cordillera) are indicated in the legend. et al., 2020). However, fauna may not be local at all, and this is especially true for livestock (Knudson et al., 2012a;Grimstead et al., 2016). Furthermore, the home ranges of small animals can be larger than expected and may crosscut discrete 87 Sr/ 86 Sr catchments (Copeland et al., 2010). Or, faunal substrate animals may have consumed non-local food, especially when they fed at or near archeological settlements that relied to some extent on imported or fertilized foods. Furthermore, animal behavior and physiology can also lead to results that do not represent the dietary catchment well; for example, land snails burrow into discrete geological pockets (Thornton, 2011) and have an extraordinarily high rainwater contribution to their shell (Evans et al., 2009(Evans et al., , 2010, leading to abnormal 87 Sr/ 86 Sr compared to other fauna or plants from the same area. Plants are another common reference material for bioavailable strontium. When sampled from various root depths, they can provide a more representative mix of the dietary catchment (Åberg, 1995;Knipper et al., 2012;Grimstead et al., 2017), but caution still must be exercised when relying on plants for several reasons. First due to the process of biopurification, 87 Sr/ 86 Sr can vary throughout different plant organs and according to plant age (Bailey et al., 1996;Drouet et al., 2005;Drouet and Herbauts, 2008). Second, different species can also show distinct 87 Sr/ 86 Sr from the same location, due in part to their root systems accessing discrete pockets of strontium in the soil (Dijkstra, 2003;Poszwa et al., 2004). Third, differential contributions from sea spray and the acidity or other molecular differences of the microhabitat can cause plant 87 Sr/ 86 Sr to differ even when the same organ and species are sampled from the same location . A fourth complicating factor is the use of non-local fertilizers in either modern or ancient times, which can contribute a non-local geological signature to plant and soil samples (Bentley et al., 2004;Bentley, 2006;Evans et al., 2012;Thomsen and Andreasen, 2019). Sea bird guano (Julien et al., 1985;Szpak et al., 2012) and manure (Finucane, 2007) were commonly applied to agricultural fields throughout the pre-Hispanic Andes. Avoiding sampling from fertilized fields is therefore essential for accurately characterizing local bioavailable strontium throughout pre-Hispanic Andes (Knudson et al., 2014). However, in many parts of the Andes, unfertilized fields are rare to non-existent, so that systematic sampling of soils and plants to avoid contamination can be challenging if not impossible. Knudson et al.'s (2014) study suggests there is much more localized variability in bioavailable Sr that could be captured through systematic sampling of plants and other reference materials in the geologically and hydrologically complex dietary catchments of the Andes.
To that end, efforts are currently underway to develop local baseline models from environmental reference materials in discrete parts of the Andes, by a variety of research teams (e.g., Barberena et al., 2017;Barberena et al., 2019, inter alia;Serna et al., 2020). However, these efforts are largely localized to sub-regions of the Andes and constrained to sampling locations near targeted archeological sites. The only supra-regional analysis of Andean 87 Sr/ 86 Sr for archeological provenience published as of June 2020 is that of the Andean Paleomobility Unification (APU) project directed by Scaffidi and Knudson (2020). Arguing that humans provide the ultimate mixing model for 87 Sr/ 86 Sr in dietary catchments, this research shows that at a supraregional scale, archeological tissue 87 Sr/ 86 Sr alone can be used to distinguish provenience between nearby areas, provided that suspected mobility routes are between discrete isoscape areas, along an east-west or coastal-highland gradient. However, this archeological isoscape does not handle equifinality well, given that people living at a similar distance from the coast but at different latitudes show similar 87 Sr/ 86 Sr. Therefore, the ongoing aims of the APU project include developing broader regional coverage and exploring variability in different 87 Sr/ 86 Sr substrates to further refine expected local baselines in the Andes.
Skeletal tissue 87 Sr/ 86 Sr at a sampling location reflects the averaging of plant and animal foods I the diet. In prehistoric communities, these are assumed to have been mostly locally procured, with the plant portion of the diet dominating the 87 Sr/ 86 Sr in the tissues of omnivorous humans (Beard and Johnson, 2000;Montgomery, 2010). In the pre-Hispanic Andes, however, the local food assumption may not hold true, and may not hold true uniformly across different individual's diets within a community. As far back as the Preceramic, reciprocal exchange of maritime products and highland grains was practiced over robust long-distance trade networks throughout the Andes (Quilter and Stocker, 1983;Hastorf, 2003; see contributions in Dillehay, 2011;Cuéllar, 2013). Particularly with the rise of the first Andean states and empires, distinct diets developed between elites and subelites, so that access to imported foods varied between individuals in communities where social inequality was manifested through differential food consumption (Hastorf and Johannessen, 1993;Hastorf, 2003;Cuéllar, 2013). For communities or individuals with access to imported foods, diet-derived 87 Sr/ 86 Sr of human tissues might depart significantly from local plant, animal, and soil inputs. Therefore, it would not be surprising for archeological human 87 Sr/ 86 Sr in the Andes to deviate substantially from any mixed model of local soil, plant, and faunal 87 Sr/ 86 Sr. However, given that the human body is up to 60% water, with up to 31% of skeletal tissues comprised by water (Mitchell et al., 1945), the 87 Sr/ 86 Sr of imbibed water may also contribute significantly to skeletal tissues. If this is the case, local drinking water 87 Sr/ 86 Sr analysis could help diminish the effect of imported foods. Therefore, this paper examines the potential for water 87 Sr/ 86 Sr to predict 87 Sr/ 86 Sr in archeological human tissues, comparing water 87 Sr/ 86 Sr and measured vs. predicted 87 Sr/ 86 Sr between samples from different elevations, watershed sizes, and underlying geologies.

Assessing Water as a Proxy for Local Skeletal 87 Sr/ 86 Sr
Water is influential in the 87 Sr/ 86 Sr of a dietary catchment and varies somewhat predictably throughout watersheds and flowing surface water networks (Aubert et al., 2002;Frei and Frei, 2013;Bataille et al., 2014;Zitek et al., 2015;Crowley et al., 2017). Because the Sr/Ca entering food chains in the form of plants are directly related to the Sr/Ca of the water they consume (Sillen and Kavanagh, 1982), dissolved strontium in surface and groundwater heavily influences plant and faunal 87 Sr/ 86 Sr (Beard and Johnson, 2000;Bentley, 2006;. Fast-flowing rivers might be suspected to carry a greater proportion of strontium from distant rocks and soils relative to slower-moving rivers, but the 87 Sr/ 86 Sr of river water has repeatedly been observed as consistent regardless of flow rate (see summary in Bentley, 2006, p. 8), seasonal, or long-term temporal variation (Voerkelius et al., 2010). In geologically active, higher elevations or during seasons of low discharge, flowing water 87 Sr/ 86 Sr is higher (more radiogenic), reflecting faster weathering of higher 87 Sr/ 86 Sr rock from higher elevation water sources throughout the watershed's headwaters, depositing highaltitude sediments on floodplains (Bentley, 2006). The 87 Sr/ 86 Sr of these sediments is then mixed into the waters consumed throughout the dietary catchment at a given sampling location (see summary in Bataille et al., 2020, pp. 5-6). Hodell et al. (2004) report that water body depth also influences 87 Sr/ 86 Srlocations with shallower water tables reflect upstream mixing, while locations with deeper water tables reflect more of the underlying bedrock signature. Watershed size also influenced 87 Sr/ 86 Sr variability, with larger watersheds leading to more homogenized 87 Sr/ 86 Sr in rivers (Brennan et al., 2019).
There are several limitations to using water 87 Sr/ 86 Sr as a reference material. First, sampling locations low in a watershed might also have greater pollution levels from upstream agricultural runoff, so higher elevation waters more closely reflect the weathering of underlying bedrock (Buhl et al., 1991;Tricca et al., 1999;Eikenberg et al., 2001;Thomsen and Andreasen, 2019). However, a recent study of Danish surface waters  shows that much of agricultural lime is retained by soils, so agricultural contamination may not influence surface water 87 Sr/ 86 Sr as much as previously thought. Second, proximity to the ocean also impacts water 87 Sr/ 86 Sr. Because modern ocean water is uniform across the planet ( 87 Sr/ 86 Sr = 0.7092) (Veizer, 1989), sea spray and brackish water can drive up otherwise low 87 Sr/ 86 Sr in coastal plains (Whipkey et al., 2000;Bentley, 2006;Snoeck et al., 2020), particularly for water samples closest to the ocean. This can lead 87 Sr/ 86 Sr in coastal regions to deviate from expectations based on bedrock. Third, Sr concentrations are generally low in global surface waters (0.06 ppm on average) (Capo et al., 1998), so water samples can be difficult to measure. Finally, at a global scale, river water 87 Sr/ 86 Sr is less variable than other substrates (Peucker-Ehrenbrink and Fiske, 2019; Bataille et al., 2020), which may average out some of the 87 Sr/ 86 Sr variability needed for identifying non-locals in zones with homogeneous bedrock geology.
Notwithstanding these potential drawbacks, surficial and groundwaters have been used as 87 Sr/ 86 Sr reference materials for provenience human samples. The earliest approach involved testing 87 Sr/ 86 Sr of water sources and comparing those to people known or suspected to reside at those locations. Over four decades ago, Nixon and Helsby (1976) showed a relationship between strontium concentrations of tap water in 12 English towns and tooth enamel of their human residents. However, a more recent study found no relationship with modern tooth enamel 87 Sr/ 86 Sr and tap water 87 Sr/ 86 Sr throughout the Netherlands (Kootker et al., 2020). In another study, Tipple et al. (2018) found that modern human hair 87 Sr/ 86 Sr reflected tap water 87 Sr/ 86 Sr. They analyzed 41 student hair samples alongside 99 tap water samples collected from sites around three schools in Salt Lake City, Utah, and found that only 6 hair samples fell outside of the tap water ranges for those locations (Tipple et al., 2018, p. 4, inferred from Figure 4). However, all tap water and hair 87 Sr/ 86 Sr fell within 0.0070 of each other (inferred from Figure 4), and the three schools shared overlapping 87 Sr/ 86 Sr, so that without known provenience data it would be impossible to distinguish between origin sites for all but individuals with the highest 87 Sr/ 86 Sr. In contrast with the Nixon and Helsby (1976) enamel study, the Tipple et al. (2018) study shows that 87 Sr/ 86 Sr in hair reflects uptake from bathing waters rather than from imbibed water.
Recent studies suggest that dietary Sr dominates human skeletal tissue 87 Sr/ 86 Sr rather than drinking water or atmospheric contributions. For example, a recent porcine controlled feeding study (Lewis et al., 2017) shows that there is little correlation between drinking water 87 Sr/ 86 Sr and skeletal 87 Sr/ 86 Sr compared to high correlations between dietary 87 Sr/ 86 Sr and skeletal 87 Sr/ 86 Sr. However, this study is limited by a small sample size (n = 19 pigs) and the assumption that human water metabolism is like that of pigs. Similarly, a large study of French children (Glorennec et al., 2016) shows that Sr concentrations are dominated by dietary Sr contributions rather than drinking water or atmospheric dust. However, this study is limited in that tissues were not directly tested for 87 Sr/ 86 Srinstead, dietary Sr concentrations were estimated from known values of foods reported in national survey data. Furthermore, the sources of water were not controlled for. Bottled water was widely consumed by study participants, so the Sr concentrations may not reflect local water. Lewis et al. (2017) further clarify that drinking water can, in fact, contribute to skeletal 87 Sr/ 86 Sr in animals that drink frequently or in landscapes with Sr-rich waters. Additional testing of human tissues relative to known food and water sources is needed to further clarify how much water 87 Sr/ 86 Sr is integrated into human bioapatite relative to dietary sources in Andean landscapes.
The potential of using water 87 Sr/ 86 Sr for archeological skeleton provenience has varied according to study location and approach. Many baseline studies including water do not directly compare water and archeological human 87 Sr/ 86 Sr (Hodell et al., 2004;Bentley and Knipper, 2005;Frei, 2011, 2013;Lengfelder et al., 2019;Ladegaard-Pedersen et al., 2020;Pacheco-Forés et al., 2020), although some incorporate water into a mixed or predictive model and evaluate archeological data points against that model (e.g., Montgomery et al., 2007;Snoeck et al., 2020). Among studies in mountainous regions comparing skeletal 87 Sr/ 86 Sr directly to local water measurements, a relationship is reported between water elevation, sample variability, and fit of archeological samples to local water. In the Rhine Valley of the Alps, Hoogewerff et al. (2001) found that high elevation stream water 87 Sr/ 86 Sr in an area comprised of old, high Rb/Sr gneisses the river's headwaters matched historic skeletal 87 Sr/ 86 Sr better than waters overlying younger, lower Rb/Sr limestone. After excluding outlier archeological pig 87 Sr/ 86 Sr, Bentley and Knipper (2005) report a similar finding for Rhine Valley sitesarcheological samples overlapped with water 87 Sr/ 86 Sr from published studies, and showed more variability in the uplands than the lowlands. So far, direct comparisons of water and archeological 87 Sr/ 86 Sr have not been published for the Andes, but similar highland vs. lowland trends are to be expected for the similarly geologically complex river valleys of the Andes mountains.
The second approach to water 87 Sr/ 86 Sr sourcing involves testing water 87 Sr/ 86 Sr throughout a region, generating an isoscape predicting 87 Sr/ 86 Sr at unmeasured regions, and then validating that model with 87 Sr/ 86 Sr measurements from tissue samples. Testing how well a water strontium isoscape for the United States derived from bedrock weathering projections and river 87 Sr/ 86 Sr measurements  predicts mammalian body tissues, Crowley et al. (2017) found that that large-bodied mammalian tissues were more accurately predicted by the surface water model than vegetation and fish skeletons. They also report that the surface water model predicted mammalian tissues better for small, relatively noncomplex watersheds than for larger watersheds. Frei and Frei (2011) found that surface water 87 Sr/ 86 Sr predictions in Denmark were consistently 15% higher in modern snail shell and soil leachates than in water samples from the same location, but modern sheep wool coincided with water 87 Sr/ 86 Sr. Lengfelder et al. (2019) used water 87 Sr/ 86 Sr in the Alps to construct an environmental mixing model which predicted 87 Sr/ 86 Sr 87 Sr/ 86 Sr in 70% (80/115) zooarchaeological samples. Lengfelder et al. (2019) considered the prediction matching if it fell within the isoscape prediction ± (2 SD + mean offset between water and archeological fauna, which was ∼0.001 for their study area). These appear to be the only studies that validate predictive models based on water 87 Sr/ 86 Sr for modern mammalian or archeological tissues. It is important to note that none of these studies tested water 87 Sr/ 86 Sr models with human skeletal 87 Sr/ 86 Sr; in light of this and recent feeding studies critiqued above (Glorennec et al., 2016;Lewis et al., 2017), understanding how much water contributes to skeletal 87 Sr/ 86 Sr is an ongoing research question.
In the present study, water was selected as a potential baseline proxy for the bioavailable strontium in archeological human tissues for several reasons. First, water samples are relatively quick to process (as described in section "Sample Collection and Strontium Isotope Analysis Methods" below). Unlike for stable oxygen isotope analysis, which is impacted by evaporation and storage of water samples, there is no evidence that 87 Sr/ 86 Sr is impacted by water storage and evaporation, so old, unparafilmed samples can still be analyzed. Second, agricultural fields in many parts of the Andes are rarely planted up to the river line, to secure crops from seasonal flooding. Due to this, water samples from rivers are suspected to be less likely than plants or soils to be contaminated by non-local fertilizer (although this may not be the case for stationary water bodies). Third, unlike foods that were sometimes exchanged across great distances in the prehistoric Andes, water was difficult to transport and would have been largely imbibed from local sources. Particularly in dry coastal valleys of the Peruvian Andes, settlements were located within walking distance of rivers, which comprised some of the only sources of drinking water in many places 2 . While the composition and geographic origin of diets often varied within and between neighboring archeological sites due to social status differences, gender, occupation, or other social identities, drinking water should have been sourced mostly from the same nearby locations. Finally, surface water provides an inherently mixed model of local and upriver bioavailable strontium catchments that averages bedrock weathering, soil runoff, atmospheric, and upriver water contributions as it moves throughout the watershed. When local foods are consumed, the overall contribution of local water to the human diet may be further amplified-i.e., directly imbibed water is disproportionately represented in human tissues, and in the tissues of local plants and animals consumed. Based on the studies above, it is hypothesized that 87 Sr/ 86 Sr will vary according to elevation, watershed size, and underlying geological age. It is also hypothesized that homogenous water at lower elevations will better predict human archeological 87 Sr/ 86 Sr compared to geologically complex higher elevation zones.

Sample Collection and Strontium Isotope Analysis Methods
Water samples were collected by co-authors and collaborators 3 during the North American summer and early fall of the 2009-2019 field seasons that comprise the drier months of the Austral winter and spring (Supplementary Table S1). We focused on collecting surface water from Southern Peru, where the majority of the co-authors' archeological samples were excavated. Because lakes, rivers, and streams would have been the most commonly used drinking water sources for archeological communities, we focused our collection on flowing surface water near archeological sites, far from fertilized fields (following the prior methods of BSIRL and APU described in Zimmer-Dauphinee et al., 2020). Since much of the municipal water stores in Peru come from diverted surficial water, we also tested tap and fountain water samples to understand how accurately they reflect surface water 87 Sr/ 86 Sr, following Chesson et al.'s (2012) investigation of tap water in the United States. Between 15 and 50 mL of sediment-free water was collected from each location and each was categorized as surface water (rivers, streams/tributaries, lakes/ponds), ground or subsurface water (springs/wells), or man-made storage or diversion of surface water (reservoirs, canals, tap water/fountains). Non-acid washed collection vials were parafilmed to prevent contamination.
The geographic coordinates, water temperature, speed, location of sampling within the water body were collected by APU collaborators in a database developed by Scaffidi using the cloud-based and collaborative Collector app (ESRI) deployed on mobile devices equipped with Bad Elf precision GPS receivers accurate to within 2 m. Coordinates collected by the ACL and BSIRL were collected with consumer to technical grade GPS receivers ranging from 10 m to 1 m in locational accuracy, and limited hydrological data was collected from sampling sites.
Finally, a literature search for published water 87 Sr/ 86 Sr in the Andes was completed. Several publications list 87 Sr/ 86 Sr for water sources in and around the Andes (Coudrain et al., 2002;Grove et al., 2003;Santos et al., 2015), but only one study published sample collection coordinates necessary for addition to this study's database (Flusche et al., 2005). Coordinates were determined for published water 87 Sr/ 86 Sr points from the North Coast of Peru (Alaica et al., 2020). Sampling locations are illustrated relative to the bedrock geological units modeled by Schenk et al. (1999) in Figure 1.
Water samples were prepared for analysis in the Archaeological Chemistry Laboratory (Center for Bioarchaeological Research) at Arizona State University following established protocols (Pacheco-Forés et al., 2020). Water was strained into metal-free vials through a 2.5 µm pore size filter paper and acidified with 5% HCl to prevent microbial growth. In the Metals, Environmental, and Terrestrial Analytical Laboratory (METAL Lab) at Arizona State University, water was evaporated and precipitates were re-dissolved in concentrated nitric acid and diluted to a 2 M stock solution. From the stock, a 1 mL aliquot was measured for elemental concentrations on the Thermo Fisher Scientific iCAP quadrupole inductively coupled plasma mass spectrometer (Q-ICP-MS). Strontium was separated through ion chromatography using a Sr-Ca ion exchange resin either with manual columns using Eichrom brand Sr Spec resin (Horwitz et al., 1992) or TODGA resin on the prepFAST automated low pressure chromatography system (Pourmand and Dauphas, 2010;Romaniello et al., 2015). For samples separated on the prepFAST system, elemental concentrations were compared for sample aliquots taken before and after strontium separation, to verify consistent yields were achieved. After separation, the strontium-containing sample was completely evaporated in a Teflon beaker and refluxed with concentrated nitric acid and 30% hydrogen peroxide to oxidize organic complexes eluted from the resin. Samples were then evaporated and reconstituted with 0.32 M nitric acid. Finally, samples were diluted with 0.32 M nitric acid in the amount necessary to achieve a concentration of 50 ppb Sr, based on the elemental concentrations of each sample from analysis on the Q-ICP-MS.
Strontium isotope ratios were measured on the Thermo-Finnigan Neptune multicollector inductively coupled plasma mass spectrometer (MC-ICP-MS) using an Apex sample introduction system. Measurements were made after on-peak blank correction, rejection of voltage outliers greater than 2 standard deviations outside the mean and normalization of 87 Sr/ 86 Sr ratios to 86 Sr/ 88 Sr of 0.1194 (Romaniello et al., 2015). Prior to any sample analysis, a sequence of SRM 987 from five to 50 ppb was analyzed to ensure measured 87 Sr/ 86 Sr was accurate within 0.00002 over the entire range. The analytical sequence included bracketing concentration-matched SRM 987 standard every six samples, and frequent secondary check standards. The SRM 987 bracketing standard was measured as 0.710247 ± 0.000019 (2 SD, n = 63) with a standard error of 0.000001. Purified IAPSO seawater gave an accurate value of 0.709169 compared to the known value of 0.709187 (Romaniello et al., 2015). A secondary check standard was created to simulate a sample with poor chemical recovery, SRM 987 run at 50% of the concentration of the bracketing standard and doped with calcium ICP solution to a Ca/Sr ratio of 500. This provided both an accurate and precise value of 0.710231 ± 0.00004 (2 SD, n = 28) with a standard error of 0.000004. In addition to these analytical check standards, NIST 1400 was purified and measured in parallel to samples to evaluate user and laboratory accuracy and precision. Aliquots of NIST 1400 similar in amount to samples was processed with each batch of samples; over the course of this project, the measured value was 0.713118 ± 0.000053 (2 SD, n = 10) with a standard error of 0.000008. This was in good agreement with the GeoREM value of 0.71314 (Romaniello et al., 2015). 87 Sr/ 86 Sr was rounded to four decimal places before spatial and statistical analyses, as this has been found to be sufficient for discriminating between origins for human tissues (Bentley, 2006;Knudson et al., 2016). Then, point values were divided into subgroupings for statistical analysis. Following prior work (Scaffidi and Knudson, 2020), groups were assigned to geoecological zones of Pulgar Vidal (1981): coastal (0 -500 masl), yunga (500 -2300 masl) and highland (2300+ masl). Point values were also divided into subgroups by geological era (Cenozoic, Mesozoic, and Paleozoic). Finally, point values were divided into groups based on watershed size, which was arbitrarily defined as small (<500,000 ha), medium (500,000 -1 million ha), and large (>1 million ha). Because 87 Sr/ 86 Sr and the residuals were nonnormally distributed according to Shapiro-Wilks tests and could not be transformed into normal distributions (see results), the non-parametric Kruskal-Wallis test was used to compare central tendencies and variability between three subgroupings. Mann-Whitney U was used to compare central tendencies between two subgroupings.

Spatial Analysis and Geostatistical Methods
All spatial and geostatistical analysis was completed with the Geostatistical Analyst tools in ArcMap 10.6.2 following prior Andean isoscape modeling methods (Scaffidi and Knudson, 2020;Zimmer-Dauphinee et al., 2020). Sampling location points (n = 263) recorded in degrees, minutes, seconds, or UTM were projected to decimal degrees (WGS, 1984) as a common coordinate space. One sample was from the Maipo Valley in Chile, and due to its distance from the main geographic focus area, it was excluded from geostatistical analysis, leaving 262 Peruvian samples. Points were plotted in ArcMap 10.6.2 and projected to the Albers Equal Area Conic coordinate system, which is necessary for accurate modeling of metric distances with the least distortion possible at the regional scale in South America. The name, distance in meters, and type of the nearest surface water body (i.e., primary river, secondary river, or lake/pond) was determined using the near tool in ArcMap. Following Scaffidi and Knudson (2020, p. 6), elevation at each sampling location was extracted from the 250 m resolution Aster Digital Elevation Model (version 2) (Tachikawa et al., 2011) and coded as coastal (0-500 masl), yunga (midelevation, 500-2300 masl), or highland (2300 masl and higher).
Finally, each sampling location was spatially joined to DEMderived drainage basins/watershed boundaries created by the HydroSHEDS project (Lehner and Grill, 2013), primary rivers (Instituto Geográfico Nacional de Perú, 2015), and the bedrock geology model of Schenk et al. (1999).
Before generating the isoscape, an analysis for suspected aberrant 87 Sr/ 86 Sr values was carried out. This is important, as a high percentage of the highland Colca Valley samples were collected from puquios or wells made with non-local cements and/or quarried stone. Furthermore, some highland samples were collected from canals flowing from higher to lower field terraces, or from areas with a history of mining contamination. These factors contributed to suspicion of non-representative 87 Sr/ 86 Sr in some parts of the highlands, necessitating analytical scrutiny. In order to assess the dataset for the possibility of contamination, following prior archeological isoscape studies (Pellegrini et al., 2016;Scaffidi and Knudson, 2020), spatial outliers identified as high-low or low-high clusters were screened for with the Anselin Local Moran's I statistic in ArcMap and were trimmed from the dataset prior to generating the isoscape.
Because the trend analysis tool in ArcMap identified a strong east-west trend from less to more radiogenic, Universal Kriging with first order trend removal was used to detrend the dataset, model the residuals, and add the trend back in to the final model. Due to the non-normal distribution (see results), 87 Sr/ 86 Sr was modeled using trans-Gaussian kriging and first transformed with the arcsin transformation. The best fit model was selected based on the following criteria: lowest mean standard error between predicted and measured 87 Sr/ 86 Sr, a root-mean square and average standard error closest to each other, a mean standardized error closest to 0, and a root-mean-square standardized error closest to 1 (Oliver and Webster, 1990;Webster, 2014, 2015). Coincidental points were not averaged before kriging and Supplementary Table S2 presents the parameters used to achieve the best fit model. The best fit model was assessed in two ways. First, a random 20% subset of the water measurement points was removed, the remaining 80% training subset was modeled on the same parameters as the best fit model and measured vs. predicted 87 Sr/ 86 Sr was cross-validated for the 20% test subset. Second, the models were tested with published 87 Sr/ 86 Sr archeological human tissue measurements at known geographic coordinates from the archeological meta-analysis compiled by Scaffidi and Knudson (2020). The archeological test set was limited to 12 sites spanning the coast, yunga, and highlands, with samples of 10 or more individuals. Trophy heads and sacrificial victims were excluded, given the likelihood that they represented nonlocal individuals. The archeological dataset was also limited to sites with a low proportion of inter-quartile range (IQR) outliers relative to the sample size analyzed (see Scaffidi and Knudson, 2020); specifically, sites with an IQR outlier proportion more than 10.0% 4 were excluded (n = 493 archeological human 87 Sr/ 86 Sr). 4 Sites with high proportions of IQR outliers include places like the Inca site of Machu Picchu, where studies have shown high variation in 87 Sr/ 86 Sr because they were settled by people from diverse geographic origins (Turner et al., 2009; see discussion in Scaffidi and Knudson, 2020). Samples from settlements comprised of a high relative proportion of probable first-generation migrants would not be Each site's IQR outliers were trimmed from the archeological test dataset prior to plotting the points and extracting predicted 87 Sr/ 86 Sr from the kriged isoscape raster, elevation, geological units, and watershed size at each site.
Cross-validation and model tests produced two metrics that are emphasized here: prediction error, which is the difference between measured and predicted 87 Sr/ 86 Sr at a given location, and prediction standard error (PSE or uncertainty) at each sampling location or archeological site. Linear regressions were also used to assess the fit between measured water and predicted 87 Sr/ 86 Sr and between measured archeological 87 Sr/ 86 Sr from the Scaffidi and Knudson (2020) dataset and water isoscape predicted 87 Sr/ 86 Sr. The 87 Sr/ 86 Sr at each site was considered a match if it fell within the isoscape 87 Sr/ 86 Sr prediction ± prediction standard error at that location. Kruskal-Wallis tests were then used to compare 87 Sr/ 86 Sr prediction standard errors between site groupings based on elevation, geology, and watershed characteristics. Finally, Chi-squared was used to compare the proportions of archeological samples falling within and beyond the geostatistical prediction standard error at each site to understand whether higher prediction rates were associated with an elevational zone.

Andean Water 87 Sr/ 86 Sr Is Related to Elevation, Bedrock Geology, and Watershed Size
Plotting water sampling locations (n = 262) 5 demonstrates that water was collected from a wide range of watersheds and underlying geologies. Water was collected from 31 distinct watersheds, reflecting sampling coverage of approximately onequarter (24.4%) of the 127 watersheds within Peru (Figure 2). The geographic coordinates, elevation, and watershed characteristics of sampling locations are reported in Supplementary Tables S1, S4. Water was collected from locations corresponding to 10 bedrock geological age units (according to the Schenk et al., 1999 model), reflecting sampling coverage of approximately two-thirds (66.7%) of the 15 bedrock age units in Peru.
Elemental concentrations of 88 Sr are sufficient to measure 87 Sr/ 86 Sr in all samples, ranging from 0.0001 -0.0113 (M = 0.0007, SD = 0.0019, n = 220) (Supplementary Table S3). New 87 Sr/ 86 Sr from this study (n = 220) are reported in Supplementary Table S4. Throughout the entire sample water 87 Sr/ 86 Sr ranges from 0.7049 to 0.7227 (M = 0.7081, SD = 0.0027, n = 262). Water 87 Sr/ 86 Sr and residuals are nonnormally distributed according to Shapiro-Wilk, W(261) = 0.672, p < 0.001. Trend analysis in ArcMap shows a strong east-west trend where 87 Sr/ 86 Sr increases longitudinally with elevation and distance from the coast. Water 87 Sr/ 86 Sr is lowest near the coasts of Lima, Moquegua, and Nasca. Water 87 Sr/ 86 Sr increases to the highest (most radiogenic) around the eastern portion of Lake expected to match the local environmental baseline as communities comprised of mostly locally born individuals, and were thus excluded from the archeological test set.  Kruskal-Wallis H(2) = 50.47, p < 0.001] (Figure 4B). Watershed size (hectares) accounts for nearly 60% of the variability in water 87 Sr/ 86 Sr (Adjusted R 2 = 0.60); as watershed size increases, so does the standard deviation within the watershed. Watershed size accounts for 43.6% of 87 Sr/ 86 Sr variation (SD) across a watershed (Adjusted R 2 = 0.42). These differences in water 87 Sr/ 86 Sr also correspond to the geological age of the underlying bedrock geological units (Supplementary Table S5). Water overlying units deposited in the most distant past, the Paleozoic, are significantly more radiogenic (Median = 0.7087, SD = 0.0034, n = 31) compared to those from the Mesozoic (Median = 0.7076, SD = 0.0027, n = 134) and the most recent Anselin Local Moran's I test generates 13 spatial outlier 87 Sr/ 86 Sr that were excluded from the isoscape model, since they were suspected to reflect contamination from agricultural runoff or storage in non-native stone wells at those sampling locations. Of these outliers, nine are low 87 Sr/ 86 Sr surround by abnormally high 87 Sr/ 86 Sr. One of these is an abnormally low value from a high-elevation (3900 masl) stream in the Ayacucho Basin. The remaining eight abnormally low 87 Sr/ 86 Sr are from around the western and northern Lake Junín area collected by the Flusche et al. (2005) study. There are also four points with abnormally high 87 Sr/ 86 Sr surrounded by more typical low 87 Sr/ 86 Sr. Three of these are from the high-elevation Colca Valley of the Department of Arequipa-one is from a natural spring, and two are from man-made irrigation canals. The final high-low outlier is from a stream in Lucanas in the Ayacucho Basin. Interestingly, none of the tap water samples (n = 5) are outliers relative to surrounding sample 87 Sr/ 86 Sr, and none of the larger surface water bodiesrivers (n = 133) and lakes (n = 34)-produced spatial outliers.
Prediction errors for the model range from 87 Sr/ 86 Sr = 0.0004 to 0.0052 (Supplementary Table S6; prediction standard error is illustrated in Figure 5C). The geographic distribution shows a lower prediction standard error at highland sites in Cusco, Ayacucho, and Huancavelica, as well as around the Lima coast (Supplementary Table S6). The geographic distribution of prediction standard errors is partially a function of sampling coverage, with lowest errors in the well-sampled Department of Arequipa and neighboring departments in southern Peru ( Figure 5B). There are no statistically significant differences in the measured-predicted 87 Sr/ 86 Sr errors between groupings based on watershed size, elevational zone, or underlying bedrock geology age. Testing the isoscape with the archeological samples from 12 sites trimmed of IQR outliers (n = 493) shows mediocre performance of the surface water-only model, explaining about one-tenth (7.0%) of the variation in human skeletal tissue 87 Sr/ 86 Sr (Adjusted R 2 = 0.07) (Supplementary Table S7). Error between the isoscape predicted and measured 87 Sr/ 86 Sr for the archeological test set ranges from 0.0000 -0.0089 (M = 0.0012, n = 493, see Supplementary Table S7). Prediction standard error (uncertainty) at each site ranges from 0.0001 -0.0052.
The only site with a prediction standard error larger than 1 SD (SD = 0.0027) is the north coast Moche settlement of Huaca Colorada (PSE = 0.0031, n = 16).
The archeological 87 Sr/ 86 Sr (trimmed of their IQR outliers) are within the site's prediction standard error of the water 87 Sr/ 86 Sr isoscape prediction for 90.0% (451/493) of the test samples (Supplementary Table S8). More archeological samples from the coast (248/251 = 98.9%) fall within their site's prediction standard error range of the water isoscape than yunga (87/103 = 84.5%) or highland samples (117/139 = 84.2%), and this difference is statistically significant, X 2 (2) = 35.081, p < 0.001, n = 493. One highland site has a particularly low proportion of skeletal 87 Sr/ 86 Sr within the model uncertainty range-none of the nine samples from Turpo fall within the prediction of 0.7092 ± 0.0019 (they are all less radiogenic than expected based on the isoscape prediction) (Supplementary Table S8).

DISCUSSION
Overall, the water isoscape mirrors the east-west distribution of 87 Sr/ 86 Sr from less to more radiogenic illustrated by the Scaffidi and Knudson (2020) archeological isoscape, as well as the bedrock geology model of Schenk et al. (1999). Both this water model and the Scaffidi and Knudson (2020) archeological isoscape lack adequate sampling coverage in the north coast and highlands, all rainforest departments (e.g., Loreto, Ucayali, and most of the Madre de Dios), and within regions in between major cities in Peru. Nonetheless, the water 87 Sr/ 86 Sr isoscape illustrates important points about the heterogeneity of water in the highlands compared to coastal and yunga water, which has critical implications for understanding the suitability of sampling likely drinking water sources as a reference material for 87 Sr/ 86 Sr geolocation of archeological individuals.
Water 87 Sr/ 86 Sr Varies According to Elevation, Watershed Size, and Geology Measured water 87 Sr/ 86 Sr fits with expectations from previous water studies in mountainous zones that show more variable 87 Sr/ 86 Sr at higher elevations (Frei and Frei, 2011;Crowley et al., 2017). The isoscape predictions vary from less to more radiogenic from east to west, with a few areas of highly radiogenic predictions around the highlands of Junín, Cusco, the north coast of La Libertad, and around the northern jungles of the Madre de Dios. The lowest 87 Sr/ 86 Sr is predicted for the Lima, Nasca, and Moquegua desert coasts ( Figure 5A). Pockets of the highest and lowest projected 87 Sr/ 86 Sr are a function of poor sampling in some regions (e.g., the north coast, the Eastern slopes and the jungle, and the central and northern highlands), and should be groundtruthed by future sampling in those areas. Indeed, the isoscape is not reliable in areas where the prediction standard error is higher than the SD of the water samples (SD = 0.0027, Figure 5B). Furthermore, the fact that the water model only explains 30% of the variance in measured water 87 Sr/ 86 Sr on cross-validation demonstrates that sampling coverage should be further improved, and also that geostatistical methods may not produce accurate models in these hydrogeologically complex mountainous areas.
Notwithstanding the poor cross-validation results, important patterns emerge that can inform our understanding of bioavailable 87 Sr/ 86 Sr in Andean waters. Highland water samples displayed distinct 87 Sr/ 86 Sr than those from midelevation (yunga) and coastal zones, which demonstrates the homogenization of water strontium sources as rivers drain into the Pacific Ocean. Heterogeneity in water 87 Sr/ 86 Sr increases with watershed size, a finding that differs from the findings of Brennan et al. (2019). In the Andes, the largest watersheds are located along the eastern slopes of the Andes and the Amazon Rainforest (Figure 2). These watersheds drain some of the oldest and most variable Precambrian-Devonian geologies (Figure 1) from the Cordillera Oriental, which produce highly variable 87 Sr/ 86 Sr water. This is important for several reasons. First, according to this isoscape, 87 Sr/ 86 Sr in the Amazon Rainforest should not be homogenous, as has been often suggested by archeologists contemplating 87 Sr/ 86 Sr results for items traded into communities from Amazonian groups. In fact, and consistent with recent hydrological analyses that report distinct 87 Sr/ 86 Sr for Amazonian rivers (Bouchez et al., 2010;Santos et al., 2015), this study suggests that 87 Sr/ 86 Sr signatures are likely to be more heterogeneous on the eastern slopes of the Andes than on the western slopes. On the other hand, these water 87 Sr/ 86 Sr results suggest that within smaller watersheds, river and stream water averages strontium from fewer distinct headwaters, so that mixed, large water bodies of lower elevation watersheds should be better proxies for local 87 Sr/ 86 Sr than the higher elevation, large watersheds. Of course, it is most prudent to generate localized and mixed models for any region, but this study shows that care should be taken to systematically sample the hydrological system particularly, in the Andean highlands and in the largest watersheds.
Water type should also be considered carefully when devising baseline sampling campaigns. As other studies have argued, it is critical to sample the most likely drinking water sources for a given archeological site (Hodell et al., 2004;Montgomery et al., 2007;Pacheco-Forés et al., 2020;Snoeck et al., 2020). This study's examination of spatial outlier water 87 Sr/ 86 Sr generally agrees with this proposition. Of the 13 spatial outliers detected with possibly aberrant 87 Sr/ 86 Sr, none were large water bodies like rivers and lakes. Most were canals and one small stream, which may have been contaminated with fertilizer from agricultural runoff, as some studies have shown (Frei and Frei, 2013;Thomsen and Andreasen, 2019). However, a more recent study by Frei et al. (2020) shows that Danish waters were not contaminated by agricultural liming, and fertilizers were fully contained in the soils. Modern and ancient Peruvian fertilization practices and their impact on soil and water 87 Sr/ 86 Sr need to be studied in-depth to understand whether this impacts the local baseline. The other spatial outliers in the present study were from natural springs or man-made wells. Springs might reflect distinct groundwater 87 Sr/ 86 Sr, while stored waters from wells or stone-lined canals may take on the geological signature of their containers, which may themselves not be local. The fact that these outliers were all from highland locations could reflect the intensity of agriculture and the proximity of fields to water sources in the highlands that make them more likely to be contaminated. Mining contamination is certainly another factor that could contribute to abnormal 87 Sr/ 86 Sr around Lake Junín (Rodbell et al., 2014), although it does not explain the unusually high 87 Sr/ 86 Sr in the Colca Canyon. Spatial outlier analysis suggests that springs, wells, and canals could reflect nonrepresentative 87 Sr/ 86 Sr. However, it is possible that seemingly abnormal well and canal 87 Sr/ 86 Sr represents true variability of bedrock strontium contributions-these areas need to be more systematically sampled to understand the complexity of the local geohydrological system.
The fact that none of the five tap water samples were spatial outliers demonstrates that tap-derived drinking water 87 Sr/ 86 Sr is close to that of the source waters they are diverted from. This pattern may have been different if tap waters from very large cities like Lima, Arequipa, and Trujillo were sampled, but in Chiclayo, Chincha, and Cusco, tap samples are within the range of nearby river 87 Sr/ 86 Sr variation. Because there were no archeological samples in the archeological test dataset that were directly comparable at the tap water sampling sites, it remains unclear whether tap water could be useful for geolocation of human skeletons. A thorough understanding of Peruvian water management practices and water transportation materials would be necessary to conduct a large-scale tap water provenience study like those conducted in the US by Chesson et al. (2012) and Tipple et al. (2018). Another consideration is that wild animals would likely have consumed primarily water from rivers or lakes, and not the tap water, smaller streams, springs, and canals closer to human settlements. Therefore, it is suspected that water contributions of wild faunal protein sources in the human food web would be better reflected by preferentially sampling large water bodies outside of anthropogenically altered settlements. Further testing of different modern and ancient drinking water sources throughout the Andes is needed to assess the utility of these dietary substrates for provenience human skeletons.

Peruvian Water 87 Sr/ 86 Sr Isoscape Predicts One-Tenth of the Variation in Archeological Human Tissues
At a broad scale, the variation in water 87 Sr/ 86 Sr from the Peruvian Andes is similar to that observed in the meta-analysis of human archeological tissues across the Andes by Scaffidi and Knudson (2020). The median 87 Sr/ 86 Sr of the water samples from this study (Median = 0.7077, n = 262) is not statistically distinct from the archeological human 87 Sr/ 86 Sr from the Andes aggregated by Scaffidi and Knudson (2020) [Median = 0.7077, n = 1390, Mann-Whitney U W(1), p = 0.315]. The present study is limited as it does not model rainwater in areas with significant precipitation, or groundwater or fog-water 87 Sr/ 86 Sr, which may have been the source of a significant amount of imbibed water in arid coastal regions where rivers frequently run dry or underground. The water-only isoscape performs well overall at predicting archeological 87 Sr/ 86 Sr-archeological 87 Sr/ 86 Sr falls within the water isoscape prediction ± prediction standard error for 90% of the IQR-outlier-trimmed archeological samples. This 90% prediction accuracy exceeds the 70% reported by Lengfelder et al. (2019) for a mixing model in the Alps based on soil, dust, groundwater, drinking water, vegetation, and fauna. The high prediction accuracy of the water-only model makes sense given that the drinking water from riverine sources reflects a mix of upstream bedrock weathering, soil and alluvial runoff, and atmospheric strontium (Whipkey et al., 2000;Hodell et al., 2004;Bentley, 2006;Crowley et al., 2017). However, the poor correlation between the measured and predicted 87 Sr/ 86 Sr still leaves most of the variation in human archeological 87 Sr/ 86 Sr to be accounted for from dietary sources of strontium.
Since R 2 values rarely approach a one-to-one correlation when modeling complex human phenomena, it is informative that 87 Sr/ 86 Sr alone accounts for nearly one-tenth of the variation in the archeological human test set 87 Sr/ 86 Sr trimmed of IQR outliers. The poor test performance could be explained by the fact that the water model accounts for 30% of the variation in water 87 Sr/ 86 Sr across the landscape when cross-validated with random water points in the 20% validation set. These low R 2 values are undoubtedly influenced by shortcomings of the geostatistical approach for modeling bioavailable strontium, including its reliance on Euclidean distances to predict 87 Sr/ 86 Sr (Bataille et al., , 2020. However, it is important to note that the fit of isoscapes to archeological or human data is always assessed based on an interval from the predictionlinear fit statistics (R 2 ) are generally not reported for human studies. This is because, as Lengfelder et al. (2019, p. 246) discuss with respect to assessing the fit between isoscape predictions and archeological values, "isotopes in bioarchaeology will always remain proxy data and how much of this 'noise' can be tolerated depends on the specific hypothesis." Nonetheless, regression-based methods are necessary to improve higher accuracy predictions than are allowed by geostatistical models (Bataille et al., , 2020Hoogewerff et al., 2019). Predicting a continuous surface from linear features like rivers is also problematic, as spatial autocorrelation is constrained by distance to those rivers (Brennan et al., 2016). Particularly in the yunga and coastal Andes where settlements were constrained to river valleys, methods that constrain predictions to networks within specific watersheds (Bataille et al., 2014;Brennan et al., 2016Brennan et al., , 2019 would likely improve prediction accuracy for sites located within a walking-distance proximity to the river networks accessed by archeological human residents. The model will also be improved with the incorporation of bioavailable 87 Sr/ 86 Sr from faunal and plant sources in local food webs that were likely consumed by inhabitants at different sites, as well as collection of baseline data from undersampled areas. Care must be used in modeling animal and plant strontium sources, however, given that non-local plants and animals likely played an important role in many Andean groups with access to far-reaching and complex exchange networks. In particular, camelids like alpacas and llamas are common protein sources in the Andes-and many of these moved from the highlands to the coast and/or vice-versa throughout their lifetimes (Núñez and Dillehay, 1979;Browman, 1980;Thornton et al., 2011;deFrance et al., 2016;Tripcevich and Capriles, 2016;Mader et al., 2018), so that individuals with diets high in camelid protein might demonstrate correspondingly complex 87 Sr/ 86 Sr that does not necessarily echo the landscape 87 Sr/ 86 Sr of origin for human end-members. Similarly, marine proteins like dried fish were often exchanged over great distances into midelevations and highlands (Lynch, 1983;Quilter and Stocker, 1983;Hastorf, 2003), which could drive radiogenic isotope signatures upwards in communities with sustained and significant access to coastal exchange routes. Even in places with access to traded foods, though, community members may not have uniformly consumed those imported goods. Furthermore, communities may have altered their local production and import of different foods based on cultural or environmental conditions. In other words, faunal distributions and subsistence practices in many parts of the Andes varied at different temporal resolutions: seasonally, from year to year, or over many generations. At present, macrobotanical and stable isotopic reconstructions of local diets are limited from the 12 archeological test sites. To more accurately model the bioavailable Sr of local dietary catchments, it is critical that archeological evidence-based local food webs and water sources are identified, sampled, and integrated into multi-scalar models. In practice, this means that future efforts at 87 Sr/ 86 Sr mixing models in the Andes will best be achieved through local and collaborative, supra-regional scale process-based models, informed and contextualized by combining ethnohistoric and archeological information on diets with geological and hydrological data. While beyond the scope of this study, multiple stable isotope data points are essential for developing appropriately weighted 87 Sr/ 86 Sr models of past dietary catchments for addressing bioarchaeological and zooarchaeological questions about mobility practices.
The water 87 Sr/ 86 Sr isoscape improves equifinality in wellsampled regions that look homogeneous in the archeological isoscape, like within Arequipa, Cusco, and the Lake Titicaca areas (Scaffidi and Knudson, 2020, p. 13). In particular, the model has a low uncertainty ( Figure 5B) and a low prediction standard error (Supplementary Tables S7, S8) for samples at Beringa in the Majes Valley, and at Chokepukio in Cusco, where both the number of archeological samples and water samples are high. In these areas, water 87 Sr/ 86 Sr and the model predictions can be relied on to identify non-locals, but in areas with uncertainty higher than the water 87 Sr/ 86 Sr standard deviation (Figure 5B), the model cannot be relied on for this and further sampling is needed. The high uncertainty in most of the Andes reflects the absence of sampling throughout areas that are more difficult to access and beyond the research focus of this study's authors in southern Peru. Due to this high uncertainty, it is emphasized here that this preliminary model cannot be relied on for geolocation, and future models will seek to incorporate additional dietary proxies for bioavailable Sr, increased the geographic scope sampling, and incorporate regression-based methods to improve the model.
The fit between archeological samples and the water isoscape prediction is better in coastal and yunga regions than in the highlands. Yunga waters show the lowest prediction standard error (Supplementary Table S6) due to homogenized, wellmixed water and high sample numbers at these elevations. Coastal water is also well-mixed, but sea-spray contributions and the preponderance of salt-rich seafood in diets may have masked local drinking water signatures in several communities (Sillen and Kavanagh, 1982;Bentley, 2006). At the coastal settlements of Ancón and Huaca Colorada where a few individuals show 87 Sr/ 86 Sr not within the site prediction standard error range (Supplementary Tables S7, S8), reliance on seafood may have driven human 87 Sr/ 86 Sr higher than the local drinking water prediction. Indeed, carbon and nitrogen data from Ancón show they consumed primarily marine proteins (Slovak and Paytan, 2011). Carbon and nitrogen data are not yet published for Huaca Colorada, but given iconographic evidence for reliance of north coastal Moche people on seafood (Donnan and McClelland, 1999; see, e.g., Bourget and Jones, 2009), it is likely that Huaca Colorada's inhabitants relied heavily on marine products, just as other coastal people did (Knudson et al., , 2015Slovak and Paytan, 2011). In contrast, at the coastal site of Armatambo, all of the human samples fell within the isoscape prediction ± prediction standard error, even though they consumed a high proportion of marine foods (Marsteller et al., 2017). At Armatambo, measured 87 Sr/ 86 Sr have a mean error of 0.0012 more radiogenic than the drinking water prediction, which may reflect a uniform elevation of 87 Sr/ 86 Sr consistent with greater consumption of marine products with 87 Sr/ 86 Sr around 0.7092. Reliance on water storage and underground aquifers may also have contributed to higher 87 Sr/ 86 Sr at these sites. Since coastal sites are generally located along arid desert seashores with little year-round freshwater, coastal and desert-dwelling people often relied on stored water or used underground wells or aqueducts (Netherly, 1984;Denevan, 2001;Schreiber and Lancho Rojas, 2003;Lane, 2014;Scaffidi, 2019). Furthermore, the use of guano fertilizers at might have driven 87 Sr/ 86 Sr up. Differences in the proportion of marine foods, and/or guano, and/or variable amounts of sea spray could contribute to prediction error, and these factors should be examined in future models.
Highland waters have the highest prediction error rates (Supplementary Table S6), which likely reflects true heterogeneity in the 87 Sr/ 86 Sr of water and underlying geology. At Ch'isi on Lake Titicaca, the fact that all archeological samples were within the isoscape prediction ± prediction standard error could confirm that most people obtained drinking water from the lake. In contrast, the failure of the water isoscape to predict any of the archeological samples from the highland site of Turpo in Andahuaylas might reflect that these people consumed water primarily from a spring or small lake that deviates from the local river system. The poor fit is to be expected when there were no water samples from this area, but the limits of geostatistical modeling and the absence of dietary substrate materials in the model are also likely at play. Overall, the higher variability of highland water 87 Sr/ 86 Sr and the potential of high-87 Sr/ 86 Sr seafood masking to impact coastal 87 Sr/ 86 Sr demonstrates that localized conditions and diets must be modeled alongside of supra-regional patterns in bioavailable substrates. Given their greater variability, localized highland water models have immense potential to discriminate between different geographic origins if geo-hydrological systems, other isotopic proxies for diet are analyzed, and likely drinking water sources are understood and systematically sampled.

Sample Collection and 87 Sr/ 86 Sr Mixing Model Development in the Andes
Systematic sampling and coordinating data collected in different coordinate systems with varying techniques remain the most significant deficiencies in baseline isoscape modeling in the Andes. Particularly along Andean coastal valleys, it is often impossible to travel in a grid, sampling at every 500 km 2 , as recommended by prior studies (Garrett et al., 1990;Fordyce et al., 1993). Furthermore, modern roads do not always pass where ancient ones like the Qhapaq Ñan (the Inca and pre-Inca road system) did, which makes many archeological sites and their isotopic catchments difficult to access. Due to these factors, more accessible archeological sites around larger cities are disproportionately well-sampled compared to more remote sites (Figure 1). To advance regional and localized baseline models, a unified, cloud-based approach would maximize efficiency in the data collection and data aggregation phases and predictive accuracy of the isoscape products. One challenge to this project was aggregating spreadsheets from multiple field sampling excursions which collected different types of locational attributes in different coordinate systems and at different spatial resolutions-points not collected in the APU relational database required manual spreadsheet consolidation, which amplified opportunities for error. Another challenge is that environmental co-variates that could be analyzed for process-based modeling, such as temperature, soil type, or contamination levels, were not collected. Widespread adoption of a single cloud-based, collaborative database with environmental co-variates suitable for the region, such as the APU database generated by Scaffidi with the ArcCollector app (ESRI), is one excellent option for multi-user collection. With an inexpensive external GPS receiver, this database can be deployed on mobile IOS and Android devices, standardizing data collection and aggregation. Data is uploaded to the cloud whenever a wi-fi signal is detected, so users can plan routes around extant points. The 2019 APU team found that even if base maps were not downloaded, Peruvian cellular data enabled sufficient bandwidth for caching base map tiles. Data is stored in a web map at ArcGIS Online, which can be viewed, edited, and downloaded by users with institutional ESRI subscriptions. The APU water sampling map can be viewed at: https://arcg.is/1bj4CT; prospective users can be added to the web map by emailing Scaffidi. Collaborative databases like those developed in other world regions (Salesse et al., 2018;Willmes et al., 2018;Hermann et al., 2020) are absolutely critical for fully characterizing isotopic variability in this difficult-totraverse region.
This study shows that a water 87 Sr/ 86 Sr isoscape is useful at a coarse scale for discerning between samples from highland, mid-elevation, and coastal provenience in the Peruvian Andes-particularly when known trade and migration networks are in well-sampled areas, pass between elevational zones, and are implicated by the research question of a given study. The greater variability in highland locations and larger watersheds calls for proportionally more baseline sampling when investigating archeological populations from these geo-hydrologically complex high-elevation zones. This study demonstrates the need for a mixing model to explain the other nine-tenths of variation in archeological human 87 Sr/ 86 Sr, and the need for standardized, collaborative sampling and data sharing. It also demonstrates the importance of analyzing environmental materials alongside paleodietary evidence to understand dietary-driven deviance from the environmental baseline. Future research will analyze bioavailable strontium from plant and faunal baselines throughout the Peruvian Andes, with the eventual goal of incorporating these multiple bioavailable strontium sources into a new, process-based mixing model for the probabilistic geolocation of archeological human 87 Sr/ 86 Sr.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
BS funded (through grants to BS and to BS, TT, and KK) and analyzed the majority of the samples, collected the samples, and wrote the manuscript. TT collected the samples and revised the manuscript. GG assisted with isotopic analysis and contributed paragraphs on analytical and instrument methods. AA and LG funded (through a grant to AA) and analyzed the some samples. AA revised the manuscript. SM, AD, and ES collected the samples and revised the manuscript. All authors contributed to the article and approved the submitted version. University BSIRL and Natasha Vang facilitated a great deal of the sample collection, organization, and shipping to Scaffidi for analysis. In addition to BSIRL, ACL, and APU students, the following individuals also graciously assisted with water sampling: Manuel Mamaní Calloapaza, Jacob Bongers, Dave Reid, Jake Dean, Alex Menaker, Matt Biwer, Donna Nash, and Kirk Costion. Furthermore, the following individuals allowed us to use their archeological research projects as home bases for sample collection: Rebecca Bria, Véronique Bélisle, Willy Yépez Álvarez, Justin Jennings, and Donna Nash. We thank researchers at the METAL Laboratory at Arizona State University for assisting with instrumentation and data processing: Trevor Martin, Tyler Goepfert, and Natasha Zolotova. Andrew Zipkin of the Archaeological Chemistry Laboratory at Arizona State University was a limitless source of feedback about chemical and spatial analysis of strontium, and the SPATIAL short course (Isotopes in Spatial Systems) from the Inter-university Training for Continental-scale Ecology group at the University of Utah provided isoscape training for Scaffidi. Michael Scaffidi assisted with figure creation. Finally, we appreciate the invitation to participate in this special edition; we also thank the two reviewers and the editors for their feedback which greatly improved the manuscript.