Making the Ancestors: Materials, Manufacturing, and Modern Replicas of Recuay Monumental Stoneworks, Ancash Highlands, Peru PDF

MAKING THE ANCESTORS: MATERIALS, MANUFACTURING, AND MODERN REPLICAS OF RECUAY MONUMENTAL STONEWORKS, ANCASH HIGHLANDS, PERU by Melissa A. Litschi B.A., University of North Carolina, Chapel Hill, 2011 M.A., Southern Illinois University, Carbondale, 2012 A Dissertation Submitted in Partial Fulfillment...

MAKING THE ANCESTORS: MATERIALS, MANUFACTURING, AND MODERN REPLICAS OF RECUAY MONUMENTAL STONEWORKS, ANCASH HIGHLANDS, PERU by Melissa A. Litschi B.A., University of North Carolina, Chapel Hill, 2011 M.A., Southern Illinois University, Carbondale, 2012 A Dissertation Submitted in Partial Fulfillment of the Requirements for the Doctor of Philosophy Degree Department of Anthropology in the Graduate School Southern Illinois University Carbondale December 2022 Copyright by Melissa A. Litschi, 2022 All Rights Reserved CHAPTER 3 MEANINGFUL MATERIALS: SOURCING RAW MATERIAL SELECTED BY RECUAY STONE SCULPTORS Introduction, Background, and Preliminary Research The selection and procurement of raw materials are foundational to the manufacturing process. Not only do the properties of raw material limit what can and cannot be accomplished for a finished product, the materials and the locations and methods of procurement carry intangible meanings that are entangled with the finished product and the activities they participate in. These meanings can be read differently by various producers and users, and they can change through time. For example, Carrara marble was famously used by Italian Renaissance sculptors for its white color, fine texture, and the translucent quality it lent to finished pieces (Riccardelli et al. 2014). Today, the artists of the Don Bosco Academy import Carrara marble to their workshop high in the mountains of Peru both for its physical qualities and for embodiment of centuries old carving traditions and associations with religiously (Catholic) inspired Renaissance masterpieces (Canto y Estilo Huaraz 2017; Poma 2012), despite the costs and abundance of locally available travertine (see Chapter 5). At the same time, Carrara marble can be purchased throughout the world for the mundane purposes of countertops and floor tile, desired for its white color, name recognition, and affordability compared to other natural stone slabs, despite its poor durability and proneness to scratching, denting, staining, and potential cracking (Sima 2019). In the Andean world, materials have been demonstrated to be an important factor in imbuing objects with meaning. As discussed in Chapter 1, research in other Andean cultural periods has identified the importance of quarries and stone origins in symbolic communications and production of sacred spaces. These meanings are best understood during the Inka period 152 thanks to ethnohistorical documents from the conquest and colonial periods that have elucidated some of the meanings behind the patterned use of monumental stone observed archaeologically. During this period, it is clear that certain stones, be they mountains, outcrops, megalithic blocks, or small oddly shaped rocks, were recognized as powerful cognizant entities, who could intercede on behalf of their communities (Arriaga 1968 ; Guaman Poma de Ayala 1988 [ca. 1615]; Lau 2008; Ogburn 2004c; Salomon and Urioste 1991; Sarmiento de Gamboa 1942 ). Recent ethnohistorical and archaeological investigations have demonstrated the importance of certain material properties in Inka and Tiwanaku stone selection practices¸ especially for monumental, state-sponsored construction. For example, the Rumiqolqa quarry has been confirmed as the source of the tired stones in present-day Ecuador that were meant to become part of the Inka’s northern capital (Ogburn 2004b, 2004c). These stones not only illustrate the sacred nature and political power of the stones, but also their agency. Stone also played a major role in the restructuring of politically and religiously important built landscapes among the Tiwanaku between CE 500-700, where grey andesite was incorporated into prominent places in megalithic construction, replacing the red sandstone for the Late Formative period particularly between CE 250-500 (Janusek et al. 2013). Similar shifts were seen in megalithic carving practices. Janusek ties the sandstone to a mythological origin point and the more distant andesite to the emergence of Tiwanaku as a regionally prominent power and connections to astrological events. Further north, stone also played a major role in the construction and renovation of Chavín de Huántar, the principal pilgrimage site for the first pan-Andean religious cult and in the creation of megalithic sculptures depicting deities and shamanic rituals (Rick 2005). Early monumental structures were primarily constructed using vertical quartzite and grey sandstone 153 slabs, while the final monumental construction phase (BCE 900-780) used white granite and sandstone consistently juxtaposed against dark grey limestone (Rodriguez Kembel 2008). The consistent juxtaposition of “black” limestone and “white” granite and sandstone indicates the either color, texture, or stone origin possibly impacted the symbolic meaning of architectural features. The use of this contrasting material in Black and White Portal has been linked to iconographic themes expressing the duality between male and female (Roe 2008). Potential stone sources are identified within 12-15km of the site, all of which required crossing a river, however no interpretation of the sourcing practices is offered (Turner et al. 1999). Turning to Recuay practices, we again see indications that stone as a material has important symbolic dimensions. The Recuay used monumental stone in a very different way from preceding groups in the region, constructing temples with unique block-and-spall masonry around outcrops of “living rock” (bedrock) and venerating stone ancestors alongside ancestral mummies. This shift in monumental construction and carving traditions marks a dramatic regional alteration in the social, political, and ideological significance of stone. Sculpted stones depicted ancestors, chiefs, warriors, and feline or avian figures, rather than the mythological deities of the past. Material choice was an important step in the production process. Granites and green andesites appear in mortuary architecture, with contrasting white granite and black basalt used in chullpas at Huaullac in the Cordillera Blanca and bright green andesites that formed a red patina common in the chullpas at Chinchawas in the Cordillera Negra. This use of brightly colored stone in mortuary settings seems to continue among Recuay sculpture, with freestanding monoliths depicting mummified ancestors often carved from white, green, and blue stones while chiefly figures and puma appear more often on grey and beige stones used as architectural embellishments (Litschi 2018; Moretti and Litschi Forthcoming), see Figures 4.6-4.8. Chapter 4 154 explores several hypotheses on the associations between sculpture form and stone type in more depth, assessing the correlations between sculpture form, stone type, and intrinsic properties related to functional and aesthetic attributes. The increasing differentiation of stone selected for architectural blocks vs sculptures of sacred entities (Lau 2016; Moretti and Litschi Forthcoming; Schaedel 1952; Turner et al. 1999) between the Chavín and Recuay periods, combine with the prevalence of stone outcrops marking the center of sacred spaces (Lau 2010a, 2010b; Terada 1979), suggests a conceptual shift, during which stone transitioned from a raw material and medium for ideological expression to a sacred entity in and of itself. The first step to elucidating the role of monumental stones in Recuay society, is to identify their sourcing practices and characterize the important factors in material choice. As demonstrated by Janusek et al. (2013), locating quarries on the landscape through geochemical comparison and identification of physical quarries allows for more nuanced analyses of the role and meanings of materials than is possible through the identification of possible sources (e.g. Turner et al. 1999). Hypotheses have existed for decades regarding the origin of both building stone and sculpted stones utilized by the Recuay people. Often building stone is assumed to be sourced from on-site outcrops (Lau 2010a, 2010b, 2016, 2021) while Cerro Walun (Mejía Xesspe 1941) and the San Franciscan andesite quarry at Pongor (Soriano Infante 1950) have been suggested as the source of the green and light grey andesites respectively. This study seeks to test the hypotheses regarding sculpted stone and to more broadly identify whether any groups of sculptures within the MAA collection share a geologic source. Sourcing of carved Recuay stone sculptures is conducted through portable X-ray fluorescence spectroscopy (pXRF) on a set of 258 Recuay sculptures and 73 geologic samples. 155 Given the limited number and range of geological samples, it is not possible to conduct a thorough analysis of possible quarries. However, samples supplied by modern stone sculptors provide geologic data for Pongor’s San Franciscan Andesite quarry and for bedrock near the site of El Castillo. Geological samples were also assessed at Cerro Walun during site reconnaissance in 2020. These select geologic data will allow the two best-known hypotheses regarding stone provenance to be assessed. Furthermore, the identification of sculptures that share a common geologic origin will help to illuminate the number of quarries utilized throughout the entirety of the Recuay sculpture making tradition and facilitate more nuanced interpretations of these practices and related social contexts. The primary research question for this geochemical analysis are: - Do visually similar stones share a geologic source? - Is Cerro Walun the geologic source for the green andesites? - Is the San Franciscan andesite quarry at Pongor the geological source for the light grey andesites with speckle biotite and quartz inclusions? - Are geologic sources related to the intended use of the stone (architectural carving vs freestanding ancestral monolith)? Theoretical Concerns regarding Methodology Field Classification of Rock Types Rocks hold the deep history of the earth. Through the study of rock composition and formation processes we can learn about how particular landscapes emerged, from their initial formation to natural and anthropogenic transformations. By studying the lithology, minerology, and petrography of rock formations we can identify formation conditions and timescales of major lithic deposits. From an archaeological perspective, the classification of rock formations can help us to reconstruct past climates and landscapes, while the classification of utilized rock is 156 the first step in the analysis of entanglements between past peoples, landscapes, and their daily lives. Rock classification helps elucidate the uses, mechanical requirements, cultural necessities, aesthetic preferences, and other intrinsic properties of potential importance for tools, building stones, and other stone objects. Additionally, classification provides important information for the design of provenance studies, which elucidates the origin of lithic materials and informs interpretations of how and why certain materials were procured for specific purposes. A variety of classification systems exist for igneous, sedimentary, and metamorphic rock based on mineralogical, elemental, or petrographic composition. The most robust systems require destructive analysis to obtain accurate measures of percentages of specific minerals or elements within a sample. Two examples include the QAPF Modal Classification which classifies volcanic rock based on the relative percentages of the minerals quartz, alkali feldspar, plagioclase, and feldspathoid (Le Maitre et al. 2002; Streckeisen 1978) and the TAS (Total Alkali-Silica) classification system which uses the relative percent weight of silica dioxide vs sodium and potassium oxide when minerals cannot be distinguished in fine-grained rocks (Le Bas et al. 1986; Le Maitre et al. 2002). In cases where sample collection and destructive analysis are not possible, field classification schemes exist to allow for preliminary classification based on lithology (i.e., color, texture, and macroscopically visible minerology of rocks) (Adam 2002; Bureau of Reclamation 1998; Earth Science Reference Tables 2010; Le Maitre et al. 2002; Lutgens et al. 2014; Rapp 2009; Travis 1955). A modified version of the 1955 Travis classification system is used as the standard field classification method by the U.S. Bureau of Reclamation, who notes that the lithological classifications are broad enough to account for slight imprecisions in mineral identification as the differences do not affect the identification of mechanical properties of the rock that effect its 157 stability, which include texture, minerology, and planes of weakness (Bureau of Reclamation 1998). Thus, this method is appropriate for the descriptive characterization of architectural and sculptural stone in archaeological contexts, whose selection would have been influenced, in part, by similar mechanical properties. These criteria are also sufficient for the identification of visually similar bedrock sources to include in geochemical analysis of lithic materials. Consistency in classifications of the study sample were established in two ways. First, classifications were independently conducted by M. Litschi and M. Brito under the direction of a Huaraz-based geologist to ensure accuracy, and second, lithological classifications were compared to petrographic classifications made by M. Brito. See the subsection “Field Classification of Rock Type” under Materials and Methods below for additional information on classification methods utilized in this study. Figure 3.1 displays the field classification schema for igneous rocks used in this study and a summary of the classification criteria are summarized in Tables 3.1 and 3.2. These classification schema are best utilized on freshly broken rock, however many archaeological stones have an altered appearance due to weathering and patina formation. Most sculptures in this study had at least one freshly broken location to compare to the field classification guides, which showcased a limited range of stone types in the assemblage. In the rare cases where no fresh breaks existed, the patinaed surface was compared to other sculptures with fresh breaks for classification. In the case of the study region near Huaraz in the north-central Andes, the local geology surrounding Huaraz is highly varied. To the west, the Cordillera Negra is dominated by the Calipuy Formation of igneous andesites, along with substantial formations of ignimbrite, sandstones, shales, limestones, and marls. To the east, the Cordillera Blanca is dominated by an expansive intrusive flow of granodiorite, which interrupts the Chicama Formation of black shale 158 and quartzite (Cobbing et al. 1981). Figure 3.2 illustrates the geographic location of major Andean rock formations and Figure 3.3 illustrates the geologic timescale. Between these two parallel mountain ranges lies a high-altitude river valley, formed by the Río Santa and its tributaries, primarily following down from the glaciers that cap the Cordillera Blanca. Many Recuay sites, like Chinchawas and Yayno, are placed on ridgetops, close to stone outcrops. Within the river valley a number of sites, like Pumakayán, sit on rocky outcrops above the alluvial plain although it is likely that some sites have been buried by alluvial deposits over the past millennia. Previously, it has been assumed that Recuay groups used the local stones as building material for their megalithic masonry based on field observations of the macroscopically visible similarities in color and inclusions, and the observations of small-scale quarrying activity at Yayno (Lau 2010c, 2011, 2016). These assumptions have not been empirically tested through any mineralogical or elemental studies. While the elemental variation of classes of rock have not been defined for this region, studies in southern Peru and Bolivia have been able to utilize pXRF to differentiate sources based on trace elements (Janusek et al. 2013; Ogburn 2004b; Ogburn et al. 2013). 159 Table 3.1. Common Rock Types in the North-Central Andes Present in the Recuay Sculpture Sample. Summarized from Adam (2002) and Lutgens et al. (2014). Rock Name Andesite Dacite Rock Type Igneous, Extrusive Igneous, Extrusive Texture Often Porphyritic Aphanitic Rhyolite Tuff Granite Granodiorite Igneous, Extrusive Igneous, Pyroclastic Igneous, Intrusive Igneous, Intrusive Aphanitic Pyroclastic Phaneritic Phaneritic Sandstone Limestone Sedimentary Sedimentary Clastic Powdery Mineral Composition Plagioclase, pyroxene Intermediate, potassium & sodium feldspar, quartz, plagioclase, hornblende, biotite Potassium & sodium feldspar, 20-40% quartz Volcanic ash with mineral, glass, and lithic fragments Potassium & sodium feldspar, 20-40% quartz Intermediate, potassium & sodium feldspar, quartz, plagioclase, hornblende, biotite Quartz Calcite Figure 3.1. Field guide for the classification of igneous rock. Drawn by M. Litschi based on Lutgens et al. 2014. 160 Table 3.2. Common Minerals in Igneous Rock. Summarized from Adam (2002) and Lutgens et al. (2014). Mineral Felsic Mineral Plagioclase Feldspar Potassium Feldspar Quartz Mafic Minerals Biotite Hornblende Olivine Pyroxene Colors Luster Cleavage/Fracture Hardness Composition black, white, grey (translucent) pink, white (opaque) many vitreous 2 at 90, with striations 6 (Ca, Na)AlSi3O8 vitreous 2 at 90 6 KAlSi3O8 vitreous conchoidal fracture 7 SiO2 brown, black green, brown green-yellow dark green, black pearly vitreous vitreous vitreous one cleavage @ 60-120deg conchoidal fracture 2 at 89 or 91deg 2.5-3 5.5 7 5-7 K(Mg, Fe)3AlSi3O10(OH)2 Ca2(Fe, Mg)5Si8O22(OH)2 (Mg,Fe)2SiO4 (Mg,Fe,Ca,Na)AlSiO3 161 Figure 3.2. Structural Map of Ancash Lithology. Map prepared by M. Litschi, lithologic data from Schenk et al. 1999. 162 Figure 3.3. Selected Ancash Chronostratigraphy, information summarized from Cohen et al. 2013, Cobbing et al. 1981, and Pfiffner and Gonzalez 2013. Fm stands for formation. 163 Principles of Non-Destructive, Portable X-Ray Fluorescence Spectroscopy Portable x-ray fluorescence (pXRF) spectroscopy permits chemical characterization of inorganic elements in a sample through surface analysis of a single point on an object. Samples can be powdered or pressed into pellets as they are in traditional benchtop XRF. Alternatively, pXRF also allows samples to be analyzed as whole artifacts for a non-destructive solution. This non-destructive approach introduces different limitations compared to traditional benchtop XRF, including the production of less robust results. The benefit is the ability to conduct analysis in the field and on objects of cultural heritage that cannot be damaged during analysis. These limitations are discussed in more detail below. A full discussion of the differences between pXRF and traditional XRF can be found in Shackley 2011. During benchtop and portable XRF analysis, a sample is bombarded with x-rays generated by the x-ray tube in the XRF unit. When these primary x-rays have a higher energy than the atom they collide with, they displace electrons from the inner electron rings of atoms in the detection area. Electrons from outer shells will then fall into the innermost electron rings to replace the displaced electrons, in this process, they emit a secondary x-ray with a wavelength unique to each element (Figure 3.4). The detector in the instrument measures and records these secondary x-rays and produces a spectrum illustrating the intensity of x-rays detected at each energy level (Figure 3.5) (Potts et al. 1997; Schatzlein 2015; Shackley 2011; Williams-Thorpe et al. 1999). See Schatzlein 2015 for an accessible overview of XRF principles. The resulting spectrum has peaks at different energy levels corresponding to different elements. The size of the peak correlates with the concentration of the element in the sample. The net area under each elemental peak can then be calculated and used in qualitative and semiquantitative analysis or calibrated through comparisons to known standards to provide elemental 164 concentrations by weight percent or parts-per-million. This data allows comparisons of the chemical make-up of any inorganic material (Shackley 2011). For this study, focused analysis of a suite of trace elements, niobium (Nb), rubidium (Rb), strontium (Sr), yttrium (Y), and zirconium (Zr), is used to identify unique chemical signatures in stones that may correlate with stone provenance. These elements have successfully been used in provenance studies of obsidian, andesite, and sandstone in the Andes (Burger et al. 2006; Janusek 2015; Janusek et al. 2013; Ogburn 2004b; Ogburn et al. 2013; Tripcevich and Contreras 2013). Concentrations of trace elements can vary locally within a single stone formation based on differences in abundance during the formation process and therefore can be used to discriminate stone provenance of visually similar stone materials (Williams-Thorpe 2008). Additionally, these trace elements are less affected by weathering processes than bulk elements like iron (Fe) and have larger critical penetration depths which permits the analysis of a greater proportion of unweathered material (Potts et al. 1997; Williams-Thorpe 2008; Williams-Thorpe et al. 1999). For example, Potts et al. (1997) demonstrate a critical penetration depth (depth from which 99% of the secondary x-rays originate) of 0.17mm for iron, 1.1mm for rubidium, and 1.6mm for zirconium. 165 Figure 3.4. Diagram of an atom showing electron response to X-Ray bombardment (Redrawn from Schatzlein 2015) Figure 3.5. Geochemical Spectra View in Bruker's S1PXRF Software showing locations of peaks for Fe, Rb, Sr, Y, Zr, and Nb. 166 Limitations of Non-Destructive pXRF and Mitigation Procedures Traditional XRF spectroscopy measures a specimen’s chemical composition through bulk analysis (typically 1mg to 1g) of a pulverized, homogenized sample. Unlike traditional XRF spectroscopy, portable XRF (pXRF) has the ability to conduct non-destructive analysis through point analysis. The aperture for the instrument used in this study is a 10mm2 oval, with higher sensitivity in the center and diminished sensitivity towards the edges of the window (Bruker Elemental n.d.). Penetration depth of the x-rays depends on the density of the sample matrix and the energy at which an element fluoresces. For the heavier elements of interest in this study, critical penetration depth can be 1-2mm, improving the instrument’s sensitivity (Bruker Elemental n.d.). Surface irregularities can also impact pXRF analysis because air attenuates xrays, altering both the number of primary x-rays that reach the sample and the number of secondary x-rays that return to the detector. To mitigate air attenuation, sample points were carefully selected on each specimen to be relatively flat and smooth without large mineral inclusions (Forster et al. 2011). Three readings were taken across each sample and averaged together to account for heterogeneity arising from the mineral inclusions in each stone. Repeated studies have shown this method to produce an approximation of a bulk sample with ±5% precision for fine-grained rocks (Forster et al. 2011; Williams-Thorpe 2008). Finally, coatings and contaminants on the sample can attenuate and bias results (Forster et al. 2011). To account for this, sample points underwent low-impact cleaning with water and a toothbrush for 15 seconds to remove lichen and other contaminants and allowed to fully airdry before readings were taken (Ogburn et al. 2013). To validate the efficacy of this approach, one specimen was selected as a control during each collection period. A reading was taken in the same location on the control specimen each day the equipment was used for this study and error rates were calculated for each element of 167 interest to determine the impact of instrument performance variability. Access to the 2018 standard was difficult in 2020, so a second standard was selected of a different stone type. To account for inter-year precision, readings were repeated on an accessible sculpture that had been analyzed in 2018. Cluster Analysis and its Applications for Geochemical Analysis In order to use geochemical data to assess stone provenance, an analytical method is needed to group specimens based on the similarities of their geochemical signatures. Cluster Analysis and Principal Components Analysis are the two most common methods in this endeavor (Liritzis and Zacharias 2011; Papageorgiou 2020). Principal component analysis (PCA) is a data transformation technique commonly used to reduce the number of dimensions in multivariate datasets. In this transformation, new variables are created that summarize the covariance between original variables (Shennan 1997). Thus the 1st component explains the maximum amount of variance in the data, with each subsequent, uncorrelated variable explaining slightly less. Two of the new variables can then be plotted on a biplot to visualize the majority of the variance within the dataset in two dimensions (Baxter 1994; Papageorgiou 2020). PCA plots are particularly useful when working with geochemical data, as most datasets include three or more elemental concentrations or intensity ratios. Additionally, once individual observations are assigned to different groups (by visual separation or a cluster analysis), confidence ellipses can be drawn for each group within the dataset that project a summary of the Gaussian distribution for each group over the plotted data (Baxter 1994, 2015; Shennan 1997). As a visualization tool, PCA is an extremely powerful aid for multivariate datasets. However, PCA alone cannot make interpretations about the statistical significance of the separation of groups on the plot (Papageorgiou 2020). Therefore, PCA should 168 be reserved as an exploratory tool and a visualization tool for the results of a validated cluster analysis. Cluster analysis was first developed by biologists as a tool for taxonomic classification (Sokal and Sneath 1963) and is now one of the most used multivariate statistics by archaeologists (Baxter 1994, 2015). Cluster analysis groups individual observations based on the similarity (or dissimilarity) of all associated characteristics (a.k.a. variables) (Sokal and Sneath 1963; Tan et al. 2019). Several types of cluster analysis exist, including hierarchical, k-means, density-based, and fuzzy (see Baxter 1994 and Tan et al. 2019 for more in-depth discussions of these methods). Hierarchical clustering can either be agglomerative, where each observation begins as its own cluster and clusters are merged using a defined clustering criterion until all clusters converge, or divisive. The output for hierarchical agglomerative cluster analysis is a dendrogram, a tree that graphically represents the merging sequence and distances (height) between nested clusters (Baxter 1994; Tan et al. 2019). K-means clustering is a prototype-based technique that gives a one level clustering solution where a predefined number of cluster centroids are randomly placed, and observations are assigned to a cluster based on their proximity to a centroid. The process of placing centroids and assigning clusters continues until the cluster solutions do not change (Tan et al. 2019). Hierarchical clustering has been criticized as an archaeological tool because the subjects of archaeological analysis do not necessarily have hierarchical relationships the way species do, at some point archaeological specimens are no longer related and should not be clustered together (Baxter 1994). For this reason, K-means clustering is often preferred by archaeologists because it does not assume a hierarchical relationship. These observations are applicable for geochemical data as well; distinct rock formation events are unrelated from each other as is the relative abundance of trace elements during the formation event. Therefore, a 169 hierarchical relationship should not exist between specimens from different lithic sources. However, in this study, K-means cluster analysis proved inappropriate because it requires a predetermined number of clusters and because it assumes equal cluster sizes. While an optimal number of clusters could be determined using hierarchical analysis as an exploratory measure, preliminary tests showed that the assumption of equal cluster sizes was an insurmountable hurdle. As suggested by Baxter (1994), hierarchical clustering can still be used when no hierarchical relationship exists because we are still operating on the assumption that observations are more or less similar to other observations on the basis of their graphical distance. In this case, the primary outcome of interest is a validated clustering partition rather than the clustering relationships depicted on the dendrogram. While hierarchical clustering is an unsupervised technique, in that no external knowledge is used to build the clustering criterion, it does require decision making and interpretation by the analyst. The analyst must select the appropriate similarity measure and clustering criterion, together known as the clustering procedure, as both affect the output of the cluster analysis (Baxter 1994; Papageorgiou 2020). Similarity is defined by the distance between samples within multidimensional space. Unlike other multivariate methods which calculate statistics based on a matrix of variables associated with each observation, cluster analysis uses this variable matrix to calculate a distance matrix that summarizes the pairwise similarity between all observations with a single distance coefficient (Shennan 1997), which can either be an actual vector distance or a correlation score (Baxter 1994). Common distance measures, along with their strengths and limitations, are described in Table 3.3. The clustering criterion is the algorithm used to join clusters based on their similarity score, methods differ in the point from which the distance between clusters is calculated. Common clustering criteria are described in Table 3.4. 170 Table 3.3. Common distance measures utilized in Hierarchical Clustering. Equations from Charrad et al. (2014) adapted from Seber (1984), descriptions from Everitt et al. (2011). Distance Measure Notes Equation* Euclidean: Usually the default 𝑑 The square root of the Large differences in one 2 𝑑(𝑥,𝑦) = √∑(𝑥𝑗 − 𝑦𝑗 ) sum of squared distances variable magnified, thus 𝑗=1 between vectors of susceptible to outliers variables Manhattan: Less influenced by large 𝑑 The sum of pairwise differences in one variable absolute difference 𝑑(𝑥,𝑦) = ∑|𝑥𝑗 − 𝑦𝑗 | between vectors of 𝑗=1 variables Maximum (Chebychev): 𝑑(𝑥,𝑦) = 𝑠𝑢𝑝 |𝑥𝑗 − 𝑦𝑗 | 1≤𝑗≤𝑑 Maximum absolute differences between vectors of variables 𝑑 Squared Euclidean: Preferred for Ward’s 2 The sum of squared method 𝑑(𝑥,𝑦) = ∑(𝑥𝑗 − 𝑦𝑗 ) distances between 𝑗=1 vectors of variables Pearson Correlation Most common Distance: Sensitive to outliers ∑𝑑𝑗=1(𝑥𝑗 − 𝑥̅ )(𝑦𝑗 − 𝑦̅) 𝑑(𝑥,𝑦) = 1 − The degree of linear Similarity of features rather 2 2 relationship between two than proximity of values, √∑𝑑𝑗=1(𝑥𝑗 − 𝑥̅ ) ∑𝑑𝑗=1(𝑦𝑗 − 𝑦̅) vectors of variables this may not be ideal for geochemical data *Where x and y are vectors of d number of variables for two different observations and j is the value of each variable 171 Table 3.4. Common clustering criterion utilized in Hierarchical Clustering. Information summarized from Everitt et al. (2011). Clustering Criterion Description Visualization Single Linkage Clusters merged based on the smallest distance between the nearest neighbors of all clusters, can lead to chaining Complete Linkage Clusters merged based on the smallest distance between furthest neighbors of all clusters Average Linkage Clusters merged based on the smallest average distance between all data points in existing clusters Centroid Linkage Clusters merged based on the smallest distance between centroids of all clusters Ward’s Method Clusters merged based on the minimum increase in the error sum of squares 172 Different clustering procedures can produce very different dendrograms (and thus clustering solutions) and it is therefore important to a) consider the most appropriate approaches for the data, and b) to repeat the cluster analysis using multiple clustering procedures. Various internal quality measures can then be employed to assess the internal structure of the resulting clusters, to select the most appropriate methodology, and to determine the optimal number of clusters from the results. These internal validation measures can help to mitigate bias from the analyst based on visual cues from the dendrogram and PCA plot. This also allows the analyst to mitigate one of the major limitations with cluster analysis, the ability for hierarchical clustering to create clusters that do not actually exist in the data (Baxter 1994). This happens either when clusters are combine/divided beyond their similar characteristics or when the clustering criterion misrepresents the structure of the data due to chaining or forcing equal cluster sizes (Baxter 1994). This study utilizes three validation measures to evaluate the most appropriate clustering procedure, cophenetic correlation, Adjusted Rand index, and Fowlkes-Mallows index. The first statistical tool to accomplish this is cophenetic correlation (Sokal and Rohlf 1962). The cophenetic correlation coefficient measures the similarity between the distance matrix and the cophenetic coefficient, defined as the dendrogramatic height at which two individuals or clusters merged (Baxter 1994; Saraçli et al. 2013; Sokal and Rohlf 1962). Thus, the cophenetic correlation coefficient can be used to assess how well the modeled clustering solution adheres to the original data. Alternatively, the cophenetic coefficients from two different clustering procedures can be compared to assess the agreement between the procedures (Baxter 1994; Saraçli et al. 2013). Baxter (1994) notes that there is some concern when utilizing cophenetic correlation to assess the clustering of archaeological data because it was built for biological data 173 and evolutionary relationships which are inherently hierarchical, an assumption that does not always hold in archaeological cases. Since cophenetic correlation considers the entire dendrogram rather than unique clustering solutions, the cophenetic correlation coefficient may be artificially lowered due to the clustering that happens above the last data-supported clustering solution (Baxter 1994). The next two validation measures used in this study to assess clustering procedures are both derived from external validation approaches, the Adjusted Rand index and the FowlkesMallows index. These two indices were developed to validate clustering procedures based on known cluster affiliations (Papageorgiou 2020), thus, the results of a clustering procedure are classified as correctly or incorrectly clustered in the same manner as errors are calculated in hypothesis testing statistics. In the case of archaeological and unprovenienced geochemical data, we do not have external data to validate our clustering analysis. Instead, we can use these two indices to conduct a pairwise comparison of multiple clustering procedures to assess overall agreement. The clustering solution that has the most agreement between clustering procedures on these two indices has the most likelihood of being true (Papageorgiou 2020). The Adjusted Rand index is an adjusted for chance index that quantifies the percentage of clustering decisions that are correct based on Type I and Type II error rates. The Rand Index provides the percentage of decisions that are correct (or the same as another solution) in a pairwise comparison of observation cluster assignment across two clustering procedures. The resulting ratio compares the number of correct classifications to the total number of classifications made (1). 𝑅𝐼 = 𝑇𝑃 + 𝑇𝑁 𝑇𝑃 + 𝑇𝑁 + 𝐹𝑃 + 𝐹𝑁 (1) 174 Using one clustering procedure as the “correct” solution, cluster assignments are assessed pairwise for all observations and are classified as true positive, true negative, false positive, and false negative. True positive (TP) occurs when all relevant observations were retrieved for a cluster assignment, true negative (TN) occurs when irrelevant data was omitted for a cluster assignment, a false positive (FP) occurs when irrelevant observations were retrieved for a cluster assignment, and a false negative (FN) occurs when relevant data is omitted from a cluster (Everitt et al. 2011; Rand 1971). The Fowlkes-Mallows index is a ratio of the positive predictive rate and the true positive rate (2). The positive predictive rate represents precision, and the true positive rate represents accuracy. Similar to the Adjusted Rand index, the Fowlkes-Mallows index uses one clustering procedure as “correct” and conducts a pairwise comparison of cluster placement for each observation, classifying each as true positive (TP), true negative (TN), false positive (FP), or false negative (FN) (Everitt et al. 2011; Fowlkes and Mallows 1983). 𝑇𝑃 𝑇𝑃 𝐹𝑀𝐼 = √𝑃𝑃𝑉 × 𝑇𝑅𝑃 = √ × 𝑇𝑃 + 𝐹𝑃 𝑇𝑃 + 𝐹𝑁 (2) A host of other internal validation approaches have been developed to help identify the optimal number of clusters from a cluster analysis. As stated above, these validation measures attempt to identify the optimal clustering solution based on maximizing the similarity of cases within clusters and maximizing the differences between clusters. Helpful measures in this endeavor include cluster cohesion (compactness) which measures the distance between cases within clusters, separation which measures distances between clusters, and connectivity which identifies the frequency with which cases are placed in the same cluster as their nearest neighbor (Handl et al. 2005; Liu et al. 2010). Charrad et al. (2014) have developed a package for RStudio (NbClust) that automates cluster validation across 30 different indices and suggests an optimal 175 cluster solution based on the most frequently returned optimal solution for each validation index (see Charrad et al. 2014 for details about the included indices). Two useful indices included in the NbClust package are the Dunn index and the Silhouette index (Charrad et al. 2014), which will be described in detail given their importance for evaluating and presenting the results from this study. All relevant equations can be found in Charrad et al. 2014. The Dunn index describes the ratio between minimum distance between clusters and the maximum distance between cases within a cluster for each clustering solution produced by a clustering procedure (except when K = 1). D= minimum separation maximum diameter (3) A simplified equation for the Dunn index (D) is given in (3), where minimum separation is the smallest pairwise distance between clusters in the solution and maximum diameter is the largest pairwise distance between cases in a shared cluster for the entire clustering solution. The larger the Dunn index value is, the more compact and well-separated clusters are in the solution. The Dunn index values can then be plotted on a line graph to easily compare the validity of different clustering solutions within a clustering procedure (Dunn 1974). The Silhouette index describes how well each case is clustered by a clustering procedure. For each clustering solution, a silhouette width is calculated for every case in the analysis. Silhouette width is defined as a ratio between the average distance to all other cases within the same cluster and the average distance to all cases in other clusters. Si = (bi -ai ) max(ai , bi ) (4) A simplified equation for the Silhouette index (Si) is given in (4), where ai is defined as the average dissimilarity between an observation i and all other observations in its cluster and bi is defined as the minimum value of the average dissimilarity between an observation i and all 176 observation in another cluster. This produces a width between 1 to -1, with values approaching 1 representing a well-clustered observation, values approaching 0 representing an observation equidistant between clusters, and values approaching -1 representing a poorly clustered observation. In general, any cases with negative values are in the wrong cluster. Individual Silhouette widths can therefore be used to assess how well classified observations are, while the average of all silhouette widths in a specific clustering solution can be compared to other clustering solutions to identify the optimal number of clusters or the appropriate clustering algorithm (Rousseeuw 1987). Materials and Methods Study Sample This study focuses on the Recuay sculpture collection housed at the Ancash Archaeological Museum (MAA) in Huaraz. These primarily unprovenienced sculptures were recovered from secondary and tertiary locations around the city of Huaraz, a project initiated by Augusto Soriano Infante in the early 1900s (Ravines 1989). Despite the analytical limitations presented by a collection with only regional provenience, this large sample of Recuay sculptures, and the availability of the sensitive and precise compositional analytical equipment, can help illuminate patterns of material choice among Recuay artisans in the Huaraz area. A total of 251 Recuay sculptures were selected for geochemical analysis, along with 3 Recuay gameboards, 4 preforms, and 3 replicas that are also housed at the museum. In addition, 73 geologic samples were analyzed to characterize regional geology and to provide the future possibility of conducting destructive analysis in order to generate a regionally specific empirical calibration for non-destructive geochemical analysis of stone artifacts. Recuay sculptures were strategically sampled to include all major sculptural forms (monoliths, lintels, slabs, and tenons) and all predominant stone types. Sampling was also impacted by the accessibility of the 177 specimens as we did not want to move any of the megalithic sculptures and needed sufficient room to safely access and analyze each with the instrument. For this reason, all sculptures displayed in the Parque Lítico and Exhibition Hall were sampled, as were easily accessible sculptures in the Repository (Deposito)48. Along with geochemical analysis, a preliminary identification of stone type was conducted through macroscopic visual examination of the mineral composition of each stone. Geologic samples include stones opportunistically selected from construction debris around the city of Huaraz, debitage from contemporary stone carving workshops, stones from the banks of the Río Quillcay and a dry riverbed (Río Seco) south of Huaraz, as well as a number of fieldstones and bedrock samples from the area around the proposed Recuay quarry at Cerro Walun. Geologic samples were selected to account for the variability of raw materials in the MAA Recuay sculpture collection and to serve as possible samples for destructive analysis in future research. For this reason, some samples including the destruction debris and river/fieldstones were selected even though they do not have any provenience. In the present analysis, these unprovenienced geologic samples could help establish overall geochemical variability throughout the region. See Table 3.5 in the results section for a summary of the sculpture forms and raw material types included in this pXRF analysis. Field Classification of Rock Type Given the inability to conduct any destructive testing in this research, rock classification for this study was conducted using field classification schemes that combine analysis of color, texture, and minerology (Le Maitre et al. 2002; Lutgens et al. 2014; Rapp 2009). See Figure 4.1 and Tables 4.1 and 4.2. Rocks were first classified by igneous (volcanic), igneous (plutonic), 48 Stone sculptures in the repository were stored on lower storage shelves, along walls, and in between rows of shelves, all of which were easily accessible for analysis. Additionally, sculptures were stored in rows in an open area of the repository. The central portion of this storage area was formed by 6 tightly clustered rows, with sculptures in the center (primarily monoliths of multiple stone types) inaccessible for pXRF or visual analysis. As monoliths were the most abundantly represented sculpture type in the sample, this should not have a negative impact on the sampling strategy. 178 sedimentary, and metamorphic. Volcanic, plutonic and sedimentary rocks were then further classified by rock name. All sculptures were examined macroscopically under dry conditions. Lighting conditions were difficult to control; sculptures displayed in the Parque Lítico were examined under full-sun or cloudy conditions depending on the day’s weather while sculptures displayed in exhibit halls and the repository were examined under dim fluorescent light. Color distortion was a concern under the fluorescent lights, so a high-powered flashlight was used to better illuminate the color of the stone matrix and mineral inclusions, if stone classification was still questionable a microscopic image was taken using a DinoLite with its own light source. Photos were taken of each sculpture using the macro setting, although colors were not always captured accurately. Microscopic photos were also taken of a selection of sculptures exhibiting common rock types by Mirko Brito using a DinoLite. Initial classifications were assigned to the Parque Lítico sculptures in 2016 (Moretti and Litschi Forthcoming) and refined in 2018 based on consultations with Mirko Brito49 (Brito 2018; Litschi 2018). A final classification of the entire sample (Parque Lítico, repository, and exhibit halls) was then conducted in 2020 to ensure there was no inter-year variability in classification decisions. The 2020 classifications were spot checked by Mirko Brito, who employed microscopic analysis at 28.6x magnification to classify mineralogical composition50 (Brito 2020; Brito and Litschi 2022), using a DinoLite AM4115ZT with polarized light, magnification range from 10x-220x, and 1280 x 1024 pixel resolution (Dunwell Tech n.d.). During this macroscopic analysis rhyolitic and dacitic rocks were difficult to differentiate so all were classified as rhyodacitic. 49 Mirko worked with a geologist at the local university (Universidad Nacional Santiago Antúnez de Mayolo UNASAM) to advise him on the volcanic and sedimentary rock types present within the collection. 50 Following methodology outlined by Roduit (2007), see Brito and Litschi (2022) for more detail. 179 Equipment and Analytical Settings Portable X-Ray Fluorescence (pXRF) spectroscopy was carried out at the Ancash Archaeological Museum during two field seasons (September 2018 and February 2020) using a Bruker Tracer III-SD. The Tracer III-SD comes equipped with a rhodium target x-ray tube and a 10mm2 X-Flash Silicone Drift Detector with a reported “resolution better than 150 eV at Mn KA1, count rate >180,000cps (Bruker Elemental n.d.).” Typically, the x-ray tube is encased in oil to shield against stray radiation, however, the particular unit used in this study was modified by the manufacturer by encasing the x-ray tube in epoxy. Replacing the oil with epoxy allowed the unit to operate at higher altitudes than the standard Tracer III-SD but it does slightly diminish the sensitivity of the instrument (Lee Drake, personal communication August 29, 2018)51. A comparison of a small subset of rock samples was run between the modified Tracer unit and an unaltered Tracer III-SD, both instruments produced raw and valid count rates in the same range. All data presented here was collected with the battery-powered unit in a handheld configuration over the manufacturer recommended 90 seconds52 using the “green” filter (6mil Cu/1mil Ti/1mil Al) to target elements from iron (Fe) to molybdenum (Mo), with a voltage of 40 KeV, current of 30 μAmps, and dry air conditions. These settings allowed the unit to excite elements with higher atomic weights, including the trace elements of interest (Rb, Sr, Y, Zr, Nb). The unit was connected to a laptop computer running Bruker’s proprietary S1PXRF software to control collection conditions and monitor spectra collection. Bruker’s proprietary X-Ray Ops software was used to specify collection conditions available in the S1PXRF controls. 51 Operating conditions in regard to altitude are not reported by Bruker for the Tracer III-SD, however, discussions with multiple Bruker representatives indicate that the Tracer III’s inability to operate at high altitude was a wellknown issue, arising in part from the reaction of the argon gas to the pressure changes at high altitude. 52 The longer the bombardment window, the more x-rays are sent and received, improving the sensitivity of the reading. This sensitivity needs to be balanced with concerns of sampling enough sculptures. A 90 second interval is suggested by the manufacturer as it produces sufficient sensitivity to reliably reflect relative elemental concentrations (Bruker Elemental, n.d.). 180 As discussed above under the “Limitations of Non-Destructive pXRF and Mitigation Procedures” subsection, the non-destructive nature of pXRF introduces three primary concerns that do not exist in destructive XRF procedures, these are 1) scatter and attenuation of x-rays due to irregular surfaces, 2) chemical weathering or contamination of the sample surface, 3) heterogeneity of the material. Effects due to irregular surfaces were controlled by always taking pXRF readings on a flat surface of the sculpture, so that asperities in the sculpture surface treatment would be the only factor impacting x-ray scatter and attenuation (Forster et al. 2011; Shackley 2008, 2012). As discussed in Chapter 4, two surface treatments were consistently utilized throughout the MAA sculpture assemblage, including a rough polish for plutonic and porphyritic andesitic rocks and a dimpled pecked surface for all other stone types, meaning that within lithological classes of rock, surface irregularities were consistent across all samples. Surficial changes in relative elemental concentrations can arise through contamination or mechanical and chemical weathering of the rock which produces a patina. For these reasons, it is recommended to take reading over freshly broken surfaces without any patina. This is not possible in non-destructive studies like the present analysis. Investigations of critical penetration depths show that the penetration of x-rays into a sample matrix is dependent on the wavelength of the x-rays. As the heavier elements of interest in this study fluoresce with higher energy photons, the critical penetration depth for these elements is deeper (1-2mm below the surface) (Bruker Elemental n.d.; Ogburn et al. 2013). Based on recent breaks observed in the Recuay sculpture assemblage, the average patina was less than 0.5mm, which should not impact the suite of trace elements used for geochemical sourcing. Issues with contamination and vegetal growth were addressed via a 15-second low-impact cleaning of all sample locations with water and a toothbrush (Forster et al. 2011; Ogburn et al. 2013). Finally, sample heterogeneity was address 181 by avoiding large mineral inclusions and averaging together the results from three readings taken across each sample to account for heterogeneity arising from the mineral inclusions in each stone (Forster et al. 2011; Williams-Thorpe 2008). Analytical Approach This study of geochemical variability within the MAA Recuay lithic sculpture collection can help to elucidate sourcing practices for Recuay stone sculpture by identifying patterns in raw material selection. Geochemical analysis via pXRF spectroscopy is used 1) to identify groups of sculptures that use the same raw materials and 2) to compare the MAA sculptures to a select sample of possible lithic sources in the Huaraz region. Typically, these questions would be tested using a quantitative approach where the elemental spectra are converted into concentration data for each element as a weight-percent for bulk elements or parts-per-million (ppm) for trace elements. Statistical analysis of concentration data can then reveal groups of objects with the same elemental signatures, indicating a shared provenance. However, quantitative analysis requires regionally specific calibration data to convert x-ray intensities into concentrations. The absence of such a calibration for this region prohibited quantitative analysis. Instead, a semiquantitative approach was utilized that employed statistical analysis of intensity ratios instead of element concentrations. As discussed above, the ratio between elements within the same pXRF reading should remain fixed, allowing provenance to be assessed (Forster et al. 2011). The sequential steps of this analysis, including initial data processing, element and ratio selection, precision validation, and statistical analysis are outlined below. Initial Data Processing. Qualitative analysis of the spectra is immediately possible using the S1PXRF software to identify the presence or absence of elements and to visually compare the relative amounts of elements between two spectra. However, to conduct semi-quantitative or quantitative analysis, processing of the data is required. A semi-quantitative approach was 182 selected for this analysis for the reason described above. In this semi-quantitative approach, Bruker’s proprietary software ARTAX was used to calculate the area under the curve for each element peak using Bayesian Deconvolution. These calculations were output into an excel file that included the total area under the curve of each element present in each observation. The area under the curve can be understood as the secondary x-ray intensity for each element. Small inconsistencies in primary x-ray counts and matrix effects (attenuation of both primary and secondary x-rays, granulometric effects from crystals in each rock sample, and differences in surface morphology where each reading is taken) between different readings make direct comparisons of x-ray intensities unreliable. However, we can assume that the primary x-ray counts and matrix effects affect each element analyzed simultaneously in the same manner (Forster et al. 2011). With this assumption, we can calculate intensity ratios of elements of interest and use those intensity ratios as the basis for discriminant analysis (De Francesco et al. 2018; De Francesco et al. 2008; Serra et al. 2017). These ratios should only be compared to other reading taken with the same pXRF unit and thus are not comparable to other studies. Element and Ratio Selection. A preliminary examination of net intensities was conducted and any reading with incongruously low numbers for every element (e.g., all element raw intensities under 100) was removed from the study as these numbers were likely the result of an error during the initial pXRF analysis, such as the equipment shutting off early. Intensity ratios for all combinations of trace elements were then calculated in Excel for each pXRF reading. These ratios were calculated independently for each stone type, placing the element with larger intensity counts on average across the sample as the numerator. This procedure allowed differences in relative concentrations to be magnified. Then, the mean intensity ratio and coefficient of variation was calculated for each elemental ratio for each sculpture to account for 183 heterogeneity in the geochemical composition of the raw material. These mean ratios were used in all statistical analyses. Samples identified as outliers on the basis the coefficient of variation were further evaluated, as a high CV could indicate a non-representative reading was included in the average. A non-representative reading could be caused by a reading of less than 90 seconds, movement of the equipment during a reading, a reading over an irregular surface, or a reading over an inclusion. The first three issues should cause high CV values for all intensity counts and ratios, while the final issue may have irregular effects across elements. If a sample has only 1 CV reading that is an outlier and the sample does not plot as an outlier in a PCA plot, it is assumed that the outlier value represents heterogeneity within the sample. Precision Validation. Next, the intensity counts and intensity ratios for the three sculptures selected as standards during analysis were analyzed to validate the precision of the instrument and to select which intensity ratios were appropriate to use in cluster analysis. The mean, standard deviation, coefficient of variation, standard error, and relative standard error were calculated in Excel using pre-programed functions. Intensity ratios with relative standard errors above 10% were deemed too imprecise to be included in the provenance study. This happened commonly with elements whose concentration was below the detection level of the instrument. Since the inter-year standard did not have repeated readings from the same location on the sculpture, a t-test was run to compare the mean counts and ratios for the three readings taken each year. Statistical Analysis. Finally, statistical analysis was conducted in RStudio using hierarchical agglomerative clustering to differentiate possible stone source groups for four classes of stone (rhyodacitic, unaltered andesitic, hydrothermally altered andesitic, porphyritic andesitic). For each stone type, the selected intensity ratios were input as variables for statistical 184 analysis. In preparation for the cluster analysis a base 10 logarithmic transformation was applied to all selected intensity ratios. This normalizes the distributions and creates equal variance across all intensity ratios (Bieber et al. 1976; Bishop and Neff 1989), preventing intensity ratios with large variances from obscuring the compositional importance of intensity ratios with smaller variances. A PCA analysis was then conducted as a preliminary visualization of the data structure. Cluster analysis was then conducted using 20 different hierarchical agglomerative clustering procedures. A dialectic process was utilized to select the best clustering procedure and optimal number of clusters. All combinations of five distance measures (Euclidean, squared Euclidean, Manhattan, Maximum, and Pearson’s correlation) and four clustering criteria (average linkage, complete linkage, centroid linkage, and Ward’s method) were assessed using cophenetic correlation. Eight clustering procedures were then selected for further analysis based on overall highest cophenetic correlation coefficient for each distance measure and clustering criterion. Optimal cluster number was then assessed for the eight clustering procedures using the NbClust package for RStudio. Optimal cluster number for each stone type was selected by the analyst based on an evaluation of the most frequently returned optimal cluster number by the NbClust package as well as the Dunn index, the Silhouette index. The Adjusted Rand index and Fowlkes-Mallows index were utilized to evaluate the agreement between different clustering solutions for the eight best clustering procedures. Finally, the optimal partitions were visualized on a PCA plot using an example from each agreeing set of clustering procedures and individual silhouette widths were assessed. This information was utilized by the analyst to select a final clustering partition to present and interpret. A future analysis of provenienced geologic samples should be used to further evaluate the solutions from this study. 185