Hydroclimate Variability Influenced Social Interaction in the Prehistoric American Southwest

When droughts and floods struck ancient agrarian societies, complex networks of exchange and interaction channeled resources into affected settlements and migrant flows away from them. Did these networks evolve in part to connect populations living in differing climate regimes? Here, I examine this relationship with a long-term archaeological case study in the pre-Hispanic North American Southwest, analyzing 4.3 million artifacts from a 250-year period at nearly 500 archaeological sites. I use these artifacts to estimate how the flow of social information changed over time, and to measure how the intensity of social interaction between sites varied as a function of distance and several regional drought patterns. Social interaction decayed with distance, but ties between sites in differing oceanic and continental climate regimes were often stronger than expected by distance alone. Accounting for these different regional drivers of local climate variability will be crucial for understanding the social impacts of droughts and floods in the past and present.


INTRODUCTION
Exchange networks are part of the broad toolkit of social and physical infrastructure humans use to manage environmental risk in social-ecological systems (Anderies, 2015). The environment can structure these exchange networks by influencing the costs and benefits of social interaction. Recent theoretical and empirical work highlights how spatial, social, and environmental factors interact with networks of exchange and interaction (Fafchamps and Gubert, 2007;Bloch et al., 2008;Nolin, 2010;Verdery et al., 2012;Freeman et al., 2014;Koster and Leckie, 2014;Hao et al., 2015;Schnegg, 2015). Distance is a key factor in such systems, making it difficult to monitor conditions in potential migration destinations (Anderies and Hegmon, 2011) and know the resources and reputation of potential interaction partners (Fafchamps and Gubert, 2007), as well as increasing the metabolic costs of transport (Drennan, 1984). For agricultural societies in water-limited environments, hydroclimate variability-changes in the balance of precipitation and evapotranspiration-may be another key factor. The benefits of interacting with others in different drought regimes can outweigh the costs of traveling longer distances. As a consequence, we might expect a greater "investment of social energy in the maintenance of social ties" between populations experiencing poorly or negatively correlated climate variability (Rautman, 1993). Norms and institutions that maintain ties between different climate regimes are likely to evolve (Durante, 2009). This process is difficult to measure in the present day due to the mismatch between the generational time scale on which cultural evolution occurs and the limited time horizons available to contemporary social sciences. Instead, we can turn to the archaeological record.
Archaeology focuses on the material correlates of human behavior and is unique in addressing how social and physical infrastructure modulate human interactions with the environment over long time spans. Not only do archaeologists catalog the remains of field systems, road networks, canals, and other components of hard infrastructure directly, but also the ceramics, raw materials, and luxury goods that are the material correlates of past networks of exchange and interaction. A powerful idea in archaeology is that, because of the interaction between societies and their biophysical environments, the spatial and temporal patterns of environmental variability can be used to predict "ideal" cultural responses and compared to archaeological observations (Halstead and O'Shea, 1989). Yet in practice it is often difficult to find archaeological data fit for purpose due to the incomplete nature of the archaeological record and the paucity of detailed paleoclimate data at the scales most relevant to human populations.
The North American Southwest is an exception. The climate of this region has been intensively studied by paleoclimatologists and climate modelers (Cook et al., 1999;Sheppard et al., 2002;McCabe et al., 2004;Herweijer et al., 2007;Cook et al., 2011;Bocinsky and Kohler, 2014;Coats et al., 2015;Routson et al., 2016;Ault et al., 2018). Its aridity aids archaeological site preservation and recovery, and nearly two centuries of survey and excavation have yielded extensive, high quality settlement pattern data (Hill et al., 2004). Detailed inventories of material culture at hundreds of archaeological sites provide an unparalleled view of the structure and dynamics of past social networks. This archaeological record attests to extensive exchange networks of durable goods such as ceramics and obsidian (Malville, 2001;Taliaferro et al., 2010;Mills et al., 2013a), and there is evidence for the long-distance transport of limited quantities of maize to the large regional center at Chaco Canyon (Benson et al., 2009;Benson, 2010). The populations of the Southwest also underwent massive social transformation, migration, and population decline in the late 13th century contemporaneous with one of the worst droughts in the last 1,000 years (Hill et al., 2004). Past work has suggested a relationship between the intensity of social interaction and patterns of drought variability, but has been limited by small sample sizes or sparse climate data (Rautman, 1993;Johnson, 1990;Cordell et al., 2007). The question is returning to the fore with the advent of high resolution climate observations and reconstructions, facilitating more detailed accounting of the spatial patterns of drought in the North American Southwest (Strawhacker et al., 2017;Strawhacker et al., 2020), and more detailed archaeological datasets . Simulations suggest that the precise nature of environmental variability is critical for exchange dynamics (Freeman et al., 2014). With these advances in our ability to map droughts in space and time comes the need to more precisely define what patterns of climate variability are actually important.
Here, I draw on hydroclimate data from the past and present to isolate specific reoccurring climate patterns, or modes of variability, in the American Southwest. I then compare these patterns to prehistoric social networks, inferred from a dataset of 4.3 million ceramic artifacts from nearly 500 archaeological sites, to examine the relationship between hydroclimate variability, distance, and social interaction over a 250-year span.

Archaeological Interaction Networks
I analyzed data from nearly 500 archaeological sites in the Southwest Social Networks (SWSN) database, a compendium of material-culture data from well-dated sites west of the Continental Divide in Arizona and New Mexico (Mills et al., 2013a;Mills et al., 2013b;Peeples and Haas, 2013;Borck et al., 2015;Hill et al., 2015;Mills et al., 2015). Version 1.0 of the SWSN database cataloged nearly 4.3 million ceramic artifacts and nearly 5,000 obsidian artifacts, providing quantitative estimates of the topology of the region-wide social network during five 50-year time steps spanning the period 1200-1450 CE (Mills et al., 2013a;Mills et al., 2015). Raw site-level ceramic counts were allocated to each time step according to an apportioning procedure that combined the occupation span of each site and the production span of each ware type with a parametric assumption of the wares' changing popularity through time (Mills et al., 2013a;Peeples and Haas, 2013).
I aggregated the point-based SWSN data into 10 km grid cells ( Figure 1A) so that the network estimates were less sensitive to local settlement dispersal or aggregation as reflected in the assemblages at individual sites (Paliou and Bevan, 2016). The choice of 10 km grid cells reflects a day's round-trip travel, bounding the area for farming and raw material collection around a site, so the procedure effectively smooths over the approximate area of each site's resource catchment (Varien, 1999;Hill et al., 2015). Then I calculated a modified version of the Jensen-Shannon divergence (Masucci et al., 2011) between the empirical frequency distributions of 15 decorated ceramic wares at each of the grid cells as where D ij is the divergence between the empirical frequency distributions of ceramic wares at sites i and j, P i is a vector of the proportions of ceramic ware type k in the assemblage at site i, and H(P) − k p k ln p k is the Shannon entropy of P measured in nats. This equation measures the similarity of two sites by the distributions of the ceramic types shared by both sites and the types exclusive to each. Analogous to the use of divergence measures in population genetics, divergence here is a (inverse) proxy for information flow. For visualization and analysis, I convert this divergence measure into a similarity metric by taking The resulting cultural similarity network ( Figure 1B) is similar to that resulting from other similarity measures such as the Brainerd-Robinson index (Mills et al., 2013a) save for different behavior in the tails of the distribution, but the Jensen-Shannon index provides a more natural interpretation as a measure of information flow. By focusing on a general measure of information flow, aggregate patterns of social interaction can be inferred regardless of the precise mechanisms of that interaction (e.g., trade, migration, shared history or raw materials). The index can thus be loosely interpreted as a probability of interaction between two sites, with identical patterns of ceramic discard at two sites indicating a higher probability of interaction than between sites that share no ware types in common.

Hydroclimate Variability
To estimate large-scale patterns of interannual drought and flood variability, I analyzed a 122-year record of the Standardized Precipitation-Evapotranspiration Index (SPEI) from across the US states of Arizona, New Mexico, Colorado, Utah, and California (Abatzoglou et al., 2017). The large spatial domain sampled variability across the western US climate zone, ensuring that estimated spatio-temporal climate patterns were not sensitive to the exact location and dimension of the archaeological study area. SPEI is the normalized deviation from the average climatic water balance for a given month on varying time scales (Vicente-Serrano et al., 2010). I focused on the 12-month SPEI calculated in the August of each year, which captures the water balance over the year leading up to each summer growing season. This index was calculated from 4 km temperature and precipitation grids interpolated from weatherstation observations using the topographically-sensitive PRISM algorithm (Daly et al., 1997).
Weather can vary for many reasons across space and time, so it is important to separate climatic signal from random noise. Principal Components Analysis (PCA) of spatiotemporal data is a common tool for extracting "modes of variability" in the climate sciences (Lorenz, 1956;Hannachi et al., 2007), but its use for this purpose is rare in archaeology (though see (Weiss, 1982;Cordell et al., 2007)). PCA of a dataset's space-time covariance matrix decomposes it into sets of orthogonal time series (principal components) and spatial patterns (eigenvectors), arranged by their contributions to the total observed variance (eigenvalues). The resulting modes of variability are an efficient means of representing a complex spatiotemporal field by a limited set of patterns.
After multiplying each grid cell by the cosine of its latitude to weight for differences cell area, I performed PCA on the stack of 122 annual SPEI maps via singular value decomposition. Then, I selected the leading modes of variability using both a scree test and North's rule of thumb, which is a method to detect degenerate patterns caused by temporal autocorrelation in the observed data (North et al., 1982). I rotated the leading PC modes using a varimax rotation in order to relax the spatial orthogonality constraints of the PCA analysis and reveal coherent, physically meaningful patterns (Richman, 1986). The resulting eigenvectors were then multiplied by the square root of the corresponding eigenvalues to yield correlation coefficients and were mapped in space. I refer to these resulting spatial patterns as empirical orthogonal functions (EOFs), and their FIGURE 1 | The Southwest Social Networks dataset, version 1.0 (Mills et al., 2013a). (A) Site locations for all time periods in the dataset, aggregated into 10 km patches to reduce biases from local settlement aggregation and dispersal. Site sizes are estimated from room counts. (B) Social networks reconstructed from ceramic assemblage similarity between each pair of sites. A similarity coefficient of one means that the sites share the same decorated ceramic wares in the exact same proportions, and a coefficient of zero means there is no overlap in the ceramic assemblages. The similarity network can be interpreted as the degree of social interaction and cultural transmission between sites, whether via migration, trade, or copying. Networks are shown over successive 50-year time spans, starting at 1200 CE.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 620856 associated time series as principal components (PCs). The PC amplitude time series were then compared to the observational record, and the signs of the eigenvalues and vectors were reversed to match the historical record (so that a positive time series value corresponded to a positive SPEI and vice versa). To determine whether these patterns were robust over time, the observed EOFs were compared to the EOFs of a SPEI reconstruction over the past millennium (Steiger et al., 2018) (see SI for details).

Least-Cost Networks
Distance ultimately constrains social interaction, as the further one travels to interact with a partner the greater will be the cost in time, energy, and other resources. In order to control for the effect of distance on social interaction, I calculated the least-cost network between all sites in the SWSN network. The topography of the study area was represented using a 90 m SRTM DEM, resampled to 250 m to reduce computation time and smooth fine-scale topographic noise. A cost matrix was calculated containing, for each DEM cell, the amount of time in seconds it would take a foot traveler to move to each of the 16 neighboring cells. Time costs were calculated using a version of Tobler's hiking function, which estimates walking speed from terrain slope. The function was modified to make it isotropic (i.e., averaging the uphill and downhill walking speeds) and adding an extra penalty to very steep slopes consistent with human cognitive biases (Pingel, 2010). This cost matrix (time) was then inverted to represent conductance (speed), facilitating a sparse matrix representation and estimation of least cost paths using efficient graph theoretic algorithms (van Etten, 2017). The resulting transition matrix was used to calculate all pairwise isotropic least cost paths between the centroids of each pair of 10 km grid cells containing archaeological materials.

Spatial Interaction Models
Spatial interaction models are used across the social and natural sciences (Wilson, 1971;Fotheringham and O'Kelly, 1989;Sen and Smith, 1995;Bavaud, 2008;Murphy et al., 2010;Head and Mayer, 2015). In a regression context, a spatial interaction model estimates the pairwise flow-of resources, migrants, or information-among entities as a multiplicative function of predictors influencing the production and attraction of flows as well as measures of their mutual separation or other generalized costs of moving. Archaeologists have used statistical spatial interaction models sparingly (Tobler and Wineburg, 1971;Hodder, 1974;Johnson, 1990) because of the rarity of archaeological data on social interaction strength, although the method is common in simulation studies where data quality is less of a restriction (Bevan and Wilson, 2013;Evans et al., 2011;Davies et al., 2014;Paliou and Bevan, 2016). The conceptual justification for the use of spatial interaction models on archaeological networks is similar to that used in molecular ecology (Murphy et al., 2010), with information flows among a spatially-structured metapopulation measured by the divergence of those populations (Mesoudi, 2018). Data of this type have three features that make traditional statistical spatial interaction modeling difficult. These are: 1) the data are bounded between zero and one; 2) the measures are pairwise symmetric; 3) we have no exact functional expectations for the specific terms in the spatial interaction model because empirical work on this scale and type is rare. To address these issues, I used a generalized additive model (Wood, 2006), a semiparametric extension to generalized linear models useful for more complex spatial interaction models (Lebacher et al., 2018).
Specifically, I fit models of the form where the logit function maps the data from [0, 1] to [−∞, +∞], t is the time step, f () is an arbitrary function estimated during model fitting using penalized cubic regression splines, τ i and τ j are time-varying random effects for the nodes incident on each edge, and ϵ is Gaussian error. This model assumes only that information flows are at equilibrium with settlement population, not that the populations themselves are at equilibrium (Wilson, 2008). The τ terms account for the non-independence of edges that share a node, and were estimated using a maximum likelihood population effects correlation structure appropriate for pairwise data (Clarke et al., 2002). I compared the AIC, BIC, and R 2 of models fit using maximum likelihood with and without the EOF terms, and refit the best performing model using restricted maximum likelihood (Clarke et al., 2002;Shirk et al., 2018).

Six Drought Patterns Explain 83% of Observed Drought Variability in the American Southwest
I used PCA to decompose the 122-year gridded observational record of western US summer moisture availability into a reduced set of spatio-temporal patterns. The leading six principal component time series together explain 83% of the variance in the observational record. The PCs represent time series that are maximally representative of the entire data set ( Figure 2). I rotated the six PCs before mapping, in order to capture more physically meaningful patterns and minimize statistical artifacts. PCs beyond the leading six were not retained for rotation and mapping, as they represent spatially and temporally incoherent variability and spurious correlations introduced by sampling error in the observational record. The same PC time series are also present in the coarse-resolution SPEI reconstruction, where they explain 96% of the reconstructed variance over the last millennium ( Supplementary Figures S5, S12, S13).

Different Drought Patterns Are Associated With Different Zones of Oceanic or Continental Influence
To reveal the latent spatial structures associated with the temporal modes of variability, I mapped the spatial patterns associated with each of the leading six PCs (Figure 3). The results are robust, recurring patterns of spatially-coherent variability, and can be interpreted as the degree to which the 122-year record at each grid cell correlates with the associated rotated PC time series. These spatial patterns are known as the (rotated) empirical orthogonal functions (EOFs). The patterns are consistent across observations and reconstructions (Supplementary Figure S9) and regardless of the exact SPEI time scale used to calculate them, which supports their overall robustness. The spatial and temporal patterns associated with the leading six PCs allows us to trace the sources of each mode of variability back to the global climate system. The origins of each drought pattern can be determined by examining the EOF maps, along with the correlations of the PCs to global sea surface temperatures and the timing of extreme dry and wet years. EOF1 reflects southwesterly flow from the tropical Pacific, bringing moisture across the low desert zones of California and Arizona. The pattern attenuates with elevation and as distance from the ocean increases. PC1 shows a broad drying trend to the present day, possibly related to increased evaporative demand due to recent warming, although the spatial pattern in the associated EOF is not itself anthropogenic. EOF2 similarly represents southeasterly flow from the Gulf of Mexico, centered on eastern New Mexico. As with EOF1, the pattern attenuates with increasing elevation and distance from the ocean due to orography and continentality, respectively. It represents cyclonic storms coming from the Gulf of Mexico, in turn influenced by variability in Atlantic sea surface temperatures. PC2 shows a major dry period centered on the Texas/New Mexico drought of 1956. EOF3 represents northerly flow associated with polar continental cold fronts, and its associated PC shows a wet peak in the 1983 Salt Lake City floods. EOF4 represents the influence of westerly flow off the Pacific Ocean and the orographic effect of the Sierra Nevada mountains intercepting this flow, and is associated with events such as the 1924 drought in California. EOF5 is centered over the great plains and attenuates across the Rocky Mountains, and was most strongly expressed during the Dust Bowl of the 1920s. EOF6 is centered on the Colorado Plateau, likely reflecting local circulation of hot continental air masses.

The Intensity of Social Interaction Decays Nonlinearly With Distance
I calculated the cost of moving between each pair of archaeological sites as the shortest amount of time it would take a foot traveler to move between them. I then used a nonlinear regression model to estimate the functional relationship between distance and interaction. The null model for the statistical network analysis was that distance alone explained the intensity of social interaction as measured by the similarity in the decorated ceramic assemblages at each pair of sites. This null model was sufficient to explain 37.8% of the variance in the ceramic similarity data. The empirical distance deterrence function estimated on all time periods using a penalized regression spline predicts a falloff in interaction at distances of more than 100 h ( Figure 4). As expected, the resulting distance-based network predicts many strong interactions at close distances, and the residuals of the model show long-distance transitive ties. FIGURE 2 | Time series associated with the leading six PCs for the observational period, after varimax rotation. The y-axis corresponds to the 12-month Standardized Precipitation-Evapotranspiration Index (SPEI), the normalized deviation from the average climatic water balance for a given month on 12-month time scale. SPEI values can be interpreted as z-scores in a normal distribution (i.e., a value of 1 is one standard deviation wetter than average for that location, −1 is one standard deviation drier). 10-year moving averages superimposed over raw annual values.

Hydroclimate Variability Explains a Moderate But Clear Proportion of the Intensity of Social Interaction
A model predicting information flow using distance and climatic dissimilarity, measured as the absolute difference between the EOF loadings of a pair of sites, explains 42.5% of the variance in the ceramic similarity data. The increase over the distance-only null model is moderate but statistically significant, and the EOF model is superior in all measures of parsimony and goodness-of-fit. This difference changes over time, and refitting each model on data from each time step individually reveals that the improvement in the explanatory power of the EOF model over the distanceonly null is most pronounced at and after 1300 CE (see SI). The improvement in explanatory power over the null model is quite small in the 1200 and 1250 CE time steps. This pattern suggests that ties shaped in part by the EOF patterns are more common during and after the period of regional relocation around 1300 CE.
The smooth functions estimated in the EOF model are all close to piecewise linear on the scale of the linear predictor, but the intensity of these functional relationships varies smoothly over time and across EOFs ( Figure 5). Increasing distance along a particular EOF sometimes increases the intensity of social interaction, as was expected ahead of time, but some EOFs (3, 6) appear to slightly reduce social interaction at larger differences. The smoothness penalty also selects some EOFs out of the model entirely by estimating functions close to a horizontal line, and almost all the functions are flat when the climate differences are less than 0.2. Surprisingly, the fluctuations in the effect size of a particular EOF have no clear association with the sign of the associated PC amplitude time series reconstructed for each period (Supplementary Figure S14). This suggests that additional dynamic processes, such as cultural memory or institutional growth and decay, are in effect on time scales longer than a single generation. These processes may explain why different sets of EOF patterns appear to influence social interaction before and after the period of drought and interregional migration ca. 1275-1300 CE.  Figure 2, derived from gridded PRISM climate data (Daly et al., 2008;Abatzoglou et al., 2017). These regions represent different oceanic and atmospheric influences; people living in the same EOF will often experience dry and wet years at the same time as one another.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 620856 FIGURE 4 | Empirical distance deterrence function estimated with a generalized additive model, describing how the intensity of social interaction, defined as the information flow between two settlements and measured by the similarity in their decorated ceramic assemblages, decreases as a function of distance. Shaded area indicates the 95% confidence interval for the smooth function. Dashed lines indicate key thresholds in the function at 10 and 100 hours.
FIGURE 5 | Estimated smooth functions describing how the intensity of social interaction increases or decreases with increasing distance along each of six spatial drought patterns from Figure 3, compared over five time steps. As above, information flow is inferred from the similarity of the decorated ceramic assemblages at each pair of sites. Climatic difference is defined as the absolute difference between the EOF loadings of each pair of sites. Shaded regions correspond to the 95% confidence intervals for the smooth functions.

DISCUSSION AND CONCLUSION
The six spatial patterns of hydroclimate variability isolated here are consistent with the general mechanistic understanding of hydroclimate variability in the American Southwest. These patterns represent different zones of moisture transport, reflecting the influence of topography and marine or continental moisture sources (Liu et al., 2010;Hu et al., 2011). These spatial and temporal drought patterns, and their hypothesized drivers in the global climate system, are largely consistent with those from other studies using varied observational data and time windows (Comrie and Glenn, 1999;Cook et al., 1999;McCabe and Dettinger, 1999;McCabe et al., 2004;Ryu et al., 2010;Seager and Hoerling, 2014;Herrmann et al., 2016). These patterns from the observational period also appear in drought reconstructions spanning the past millennium, emphasizing the fact that these are robust, time invariant spatial modes. Objective measures of hydroclimate variability, as opposed to point-to-point sample correlations, help isolate the most important drivers of that variability. Droughts and pluvials associated with tropical Pacific and Atlantic influences seem to have been most important for structuring social interaction, with ties connecting these regions greater than expected by chance and distance alone. Tropical Pacific sea surface temperatures are known to be the primary driver of variability in Southwest, with additional influences from moisture sources in the North Pacific and Atlantic (McCabe et al., 2004). A disruption in these patterns is thought to be one reason why droughts in this period led to such social transformation, as the networks of social infrastructure that had developed over previous centuries were unable to adapt fast enough to unusual conditions (Cordell et al., 2007).
In spite of the robustness of these spatial patterns, there remains considerable diversity in the functional responses of human social networks to these drought patterns. One possible explanation is that large-scale climate regimes influence the formation of ethnolinguistic groups. Quotidian interaction may have been biased toward groups of shared ethnolinguistic affiliation, as kinship and ethnicity both influence social exchanges (Nolin, 2010). At larger, regional scales, goods and information might also flow on sociopolitical hierarchies (Crumley, 1979). Populations in the late pre-Hispanic Southwest were also out of equilibrium (Hill et al., 2004). Strong social networks take time to form and effort to maintain and monitor. Free-riders who avoid that effort can damage this critical social infrastructure when it is most needed (Kohler and West, 1996). A simulation approach could better capture these processes and more clearly resolve social responses to interannual climate variability. Dynamic, as opposed to statistical, spatial interaction models can explicitly trace the coevolution of social and physical infrastructure networks (Bevan and Wilson, 2013). Simulations can also explore the biases that static archaeological data introduce in representing dynamic social processes (Crema et al., 2014(Crema et al., , 2016. The residuals of the statistical network models retain unexplained structure. These structures appear to represent cultural clusters, a common feature in social networks that is not accounted for by either distance or drought variability. One source of this error may be irrigation-dependent groups who relied on streamflow driven by remote precipitation and evapotranspiration and may thus have had more complex, indirect dependencies on the large-scale climate patterns isolated here (Nelson et al., 2010). At finer scales, the model residuals also display evidence of transitivity and triad closure, with more closed triangle structures than would be expected by chance. Although this feature is common in human social networks, it is also to be expected because the measure of cultural similarity is a metric subject to the triangle inequality. Statistical methods specifically designed for such structures will be useful in future work (Stillman et al., 2017). In addition, the archaeological data are not spatially extensive enough to sample the full range of hydroclimate variability. Given the relative spatial scales of the environmental and cultural data, there is a risk that many different correlated climate patterns will be indistinguishable at the scale of the cultural data. Correlations between competing hypotheses are a source of error in model selection using information criteria (Shirk et al., 2018). In spite of these concerns, these results highlight two key points: the need to use objective and physically meaningful measures of hydroclimate to assess the social impacts of climate variability in the past and the need to capture social dynamics out of equilibrium with the biophysical environment.
These results refine our understanding of the geography of human adaptation to climate and climate change, and emphasize the role of social interaction in increasing the robustness of human populations to environmental variability. Much of the world's food is still grown on small farms, and these farmers rely on complex spatial networks of formal and informal arrangements in much the same way as their predecessors have for thousands of years. Prehistoric exchange infrastructure evolved in part in response to robust, timeinvariant spatial climate patterns. But social adaptions to one mode of variability are fragile to changes in the nature of that variability (Janssen et al., 2007). Large-scale patterns of hydroclimate variability act as a dynamic selective environment in which societies evolve new norms and institutions for regulating social interaction. Tracing the flows of information and energy within these complex social-ecological systems is essential for understanding their long-term behavior, and leveraging our archaeological understanding of why such systems succeed or fail will be critical to anticipating the impact of impending climate changes on farming communities in the developing world.

DATA AVAILABILITY STATEMENT
All original contributions presented in the study are included in the article or Supplementary Material in a reproducible R Markdown document. A preprocessed version of the SWSN v1.0 data, aggregated to 10 km grid cells to mask individual site locations, is available in the Supplementary Material with links to additional input climate and topographic raster data. Further inquiries can be directed to the corresponding author.