Formal Tests for Resistance-Resilience in Archaeological Time Series

Riris, Philip; de Souza, Jonas Gregorio

doi:10.3389/fevo.2021.740629

ORIGINAL RESEARCH article

Front. Ecol. Evol., 13 December 2021

Sec. Interdisciplinary Climate Studies

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.740629

This article is part of the Research Topic Extreme Events in Human Evolution: From the Pliocene to the Anthropocene View all 11 articles

Formal Tests for Resistance-Resilience in Archaeological Time Series

$\r\nPhilip Riris*$ Philip Riris^1* $Jonas Gregorio de Souza\r\n$ Jonas Gregorio de Souza²

¹Institute for the Modelling of Socio-Environmental Transitions, Bournemouth University, Bournemouth, United Kingdom
²Culture and Socio-Ecological Dynamics Research Group (CaSEs), Universitat Pompeu Fabra, Barcelona, Spain

The study of resilience is a common pathway for scientific data to inform policy and practice towards impending climate change. Consequently, understanding the mechanisms and features that contribute towards building resilience is a key goal of much research on coupled socio-environmental systems. In parallel, archaeology has developed the ambition to contribute to this agenda through its unique focus on cultural dynamics that occur over the very long term. This paper argues that archaeological studies of resilience are limited in scope and potential impact by incomplete operational definitions of resilience, itself a multifaceted and contested concept. This lack of interdisciplinary engagement fundamentally limits archaeology’s ability to contribute meaningfully to understanding factors behind the emergence and maintenance of long-term societal resilience, a topic of significant interest that the field is in theory ideally positioned to address. Here, we introduce resilience metrics drawn from ecology and develop case studies to illustrate their potential utility for archaeological studies. We achieve this by extending methods for formally measuring resistance, the capacity of a system to absorb disturbances; and resilience, its capacity to recover from disturbances, with a novel significance test for palaeodemographic data. Building on statistical permutation and post-hoc tests available in the rcarbon package in the R statistical environment, we apply our adapted resilience-resistance framework to summed probability distributions of calibrated radiocarbon dates drawn from the Atlantic Forest of eastern Brazil. We deploy these methods to investigate cross-sectional trends across three recognised biogeographical zones of the Atlantic Forest domain, against the backdrop of prehistoric phases of heightened hydroclimatic variability. Our analysis uncovers novel centennial-scale spatial structure in the resilience of palaeodemographic growth rates. In addition to the case-specific findings, we suggest that adapting formal metrics can help archaeology create impact and engagement beyond relatively narrow disciplinary concerns. To this end, we supply code and data to replicate our palaeodemographic analyses to enable their use and adaptation to other archaeological problems.

Introduction

In this paper, we use radiocarbon frequency data from the Atlantic Forest of eastern Brazil to develop an empirically grounded archaeological test for resistance-resilience, focused on the period of the Medieval Climate Anomaly (MCA). The MCA is a multi-centennial period of anomalous warming between ∼1,000 and 700 cal BP (AD 950–1250), originally identified in northern hemisphere palaeotemperature proxies (Bradley et al., 2003; Mann et al., 2009). Narrowing data gaps across tropical South America have found support for its expression as a period of abrupt hydroclimate change over central and eastern Brazil too, with precipitation proxies suggesting depressed rainfall due to a northward shift of the Intertropical Convergence Zone (ITCZ) and lowered moisture delivery across the continent (Novello et al., 2018; Deininger et al., 2019; Lüning et al., 2019). Following the MCA, the Little Ice Age (LIA) manifests in the northern hemisphere as a period of cooling, which in South America shows a signal that varies from wet to dry, depending on the position within the path of the South American Monsoon System (SAMS). Of particular interest in the Atlantic Forest biome is a dry-wet dipole to the northeast and southwest of the South Atlantic Convergence Zone (SACZ), respectively (Novello et al., 2018; Utida et al., 2020). Records from the core of the SACZ, which lies astride the central Atlantic Forest, show a strong multidecadal to centennial-scale variability at around ∼900–500 cal BP, on the cusp of the transition between MCA and LIA (Novello et al., 2018). We expand the scope of recent palaeodemographic modelling work in southern Brazil (De Souza and Riris, 2021) and neighbouring regions of South America (Azevedo et al., 2019; De Souza et al., 2019) to encompass the entire Atlantic Forest biome across eastern Brazil and evaluate the resilience of pre-Columbian populations to known climate perturbations that affected this domain.

The Atlantic Forest, originally covering 1.5 million km² over a wide latitudinal and elevation range, is recognised as a global biodiversity hotspot, composed of evergreen forests, dry forests, mangroves, and other vegetation formations with high endemism. It is also one of the most endangered Neotropical biomes, with only ca. 12% of its former extent remaining at most (Ribeiro et al., 2011; Marques et al., 2021). In addition, the Atlantic Forest has a long history of human-environmental interaction. At the time of European contact, most of the coastal strip and hinterland of the La Plata basin was inhabited by Tupi-Guarani societies, who practised polyculture and management of secondary forests, as attested by ethnohistorical and archaeological evidence (Scheel-Ybert et al., 2014; Bonomo et al., 2015). This follows millennia of landscape modification by marine-adapted shell mound (sambaqui) societies, who practised a mixed economy with forest management and some cultivation (Scheel-Ybert and Boyadjian, 2020). The Tupi-Guarani presence on the coast also overlaps with the arrival of the Aratu tradition, who occupied large circular villages in central Brazil, but whose expansion towards and impact on the Atlantic Forest are still poorly understood (Robrahn-González, 1996). At the same time, the southern Brazilian highlands in particular are increasingly recognised as an anthropogenic landscape, where a peak in human occupation (by the Southern Jê/Taquara-Itararé Tradition) coincides with the expansion of Araucaria forests—one of the vegetation types in the Atlantic Forest biome (Iriarte and Behling, 2007; Robinson et al., 2018; Pereira Cruz et al., 2020)—and the emergence of a mixed agroforestry/horticultural subsistence pattern (Corteletti et al., 2015).

The growing alignment of archaeological and palaeoecological datasets in tropical South America (Mayle and Iriarte, 2014; Iriarte et al., 2020) has helped to underscore the importance of forests to Indigenous subsistence. In what has come to be termed “polyculture agroforestry systems,” a subsistence base made up of an extensive mix of cultivated, domesticated, and wild sources (Arroyo-Kalin, 2012; Clement et al., 2015; Roberts et al., 2017; Maezumi et al., 2018; De Souza et al., 2019; Iriarte et al., 2020) with a substantial input of aquatic resources (Albuquerque and Lucena, 1991; De Borges Franco, 1998; Wüst and Barreto, 1999; Hermenegildo et al., 2017; Azevedo et al., 2019; Prestes-Carneiro et al., 2019), pre-Columbian groups are inferred to have been more resilient when compared to contemporary intensive land use systems (De Souza et al., 2019). Furthermore, humid/arid cycles are a known driver of Atlantic Forest composition and extent (Scarano and Ceotto, 2015; Ledru et al., 2016; Barros Silveira et al., 2019). Over the MCA and early LIA period, archaeological data in Amazonia and the central Brazilian Cerrado evidence heightened levels of conflict and site abandonment (Azevedo et al., 2019; De Souza et al., 2019). While the heterogeneity of late prehistoric palaeoecological trends cautions against simple mechanistic models of culture change (Bush et al., 2021), it is worth noting that by 750 cal BP population growth in neighbouring Amazonia had likely, at best, stabilised (Arroyo-Kalin and Riris, 2021). We suggest that the reliability and availability of resources under polyculture agroforestry systems may also have been affected by hydroclimatic perturbations and propose this as a mechanism through which abrupt climate change impacted pre-Columbian population trajectories in the Atlantic Forest. Our working hypothesis is that the expression of the MCA across the Atlantic Forest adversely impacted the resource base of pre-Columbian Indigenous groups practicing polyculture agroforestry, which we anticipate is reflected in proxies for ancient population change over time.

To evaluate this hypothesis, we present a novel test of resistance-resilience (Nimmo et al., 2015) specifically developed for summed probability distributions of calibrated radiocarbon dates (SPDs). SPDs are one among several methods that assume the relative abundance of radiocarbon dates will reflect broad trends in “events” over time, as a proxy for relative population change in the past (Shennan et al., 2013; Timpson et al., 2014; Crema and Bevan, 2021). SPDs are, at present, one of the most common means of aggregating radiocarbon dates. Due to this, the study of population dynamics in the past—palaeodemography (French et al., 2021)—is by necessity performed within one or several statistical frameworks to account for the chronometric uncertainties, biases, and sampling problems that are inherent to such data. A major goal of these approaches has been to provide robust alternatives to visual interpretations of SPDs, either through alternative summation methods or the development of statistical models to control for diverse sources of error (Timpson et al., 2014; Crema et al., 2016; Bronk Ramsey, 2017; McLaughlin, 2019; Carleton and Groucutt, 2021). We build upon one of the most widespread such frameworks, available in the R package rcarbon (Crema and Bevan, 2021), to implement our approach towards resistance-resilience in the Atlantic Forest biome of eastern Brazil. It affords the ability to control for variation in sampling intensities through binning, mitigation of calibration curve effects on SPD shape, and direct statistical comparison of multiple SPDs concurrently through mark permutation tests (Crema and Bevan, 2021). Although other readily available approaches, such as kernel density estimates or event count modelling (Brown, 2017; Carleton, 2021) are also theoretically amenable to the treatment we propose here, they do not currently have the same range of model-testing functionality that rcarbon offers. We also scaffold our approach with Bayesian model-fitting methods (Crema and Shoda, 2021) to build a palaeodemographic context for the present study that iterates on our prior work in the region (De Souza and Riris, 2021).

Our aims are to: (1) develop and apply palaeodemographic tests for resistance-resilience, which we define below, to enable the systematisation of knowledge on this crucial topic, and (2) demonstrate their utility for statistically comparing responses to disturbances between different regions of the Atlantic Forest while evaluating their implications. Before presenting the methods, we draw on recent reviews of resilience in archaeological research (Bradtmöller et al., 2017; Fitzhugh et al., 2019; Degroot et al., 2021; Russo and Brainerd, 2021) to contextualise our method and outline the scope of palaeodemography’s potential contribution to the study of resilience in ancient socio-environmental systems. More broadly, we situate this work within the study of “extreme events” in the deep past, bearing in mind the absence of consensus definitions for extremity (Intergovernmental Panel on Climate Change [IPCC], 2012; Gregory et al., 2015), and turn our attention to a systems approach that is focused on empirical benchmarking of culture-climate links, through the lens of palaeodemography (Broska et al., 2020).

Resilience and Palaeodemography

Archaeological approaches to resilience are greatly influenced by the work of Gunderson and Holling (2002) on panarchy and adaptive cycling (Redman, 2005; Fisher et al., 2009; Kintigh et al., 2014). Recent reviews indicate that many, if not most, archaeological studies of ancient resilience draw on Resilience Theory (RT) in this sense (Holling, 1986; Folke, 2006; Bradtmöller et al., 2017; Fitzhugh et al., 2019). Within this body of work, adaptive cycling describes socio-environmental systems consisting of repeating phases of growth, conservation, release, and reorganisation in a domain of interest, which can be subsumed within (or themselves contain) cycles that operate across different spatio-temporal scales (Gunderson and Holling, 2002; Middleton, 2017). Due to the challenges associated with identifying appropriate scales and resolving different phases of adaptive cycling in the material record, studies of prehistoric resilience tend to adopt alternative strategies. Methods found within the archaeological resilience literature range from making qualitative point estimates or descriptions of conditions before and after a disturbance, to visual cross-referencing different proxies or time series (“wiggle matching,” cf. Bronk Ramsey et al., 2001; Blockley et al., 2018), as well as more rigorous quantification of relationships between different elements of a socio-ecological system (Kelly et al., 2013; Kintigh and Ingram, 2018; Carleton and Collard, 2020; DiNapoli et al., 2021).

In the vast majority of cases, resilience is employed largely as an interpretative device (Bradtmöller et al., 2017, p. 4), and consequently, the recognition of resilience (or vulnerability, stability, or fragility) in the archaeological record rests largely upon qualitative criteria (Degroot et al., 2021). Although this may offer internally consistent perspectives on how a socio-environmental system changed in response to disturbances, the findings of any single case study are rarely meaningfully comparable with those of another. Similarly, palaeodemographic studies tend to employ resilience heuristically (Leppard, 2015; Marsh, 2016; Bevan et al., 2017; Blockley et al., 2018; Riris and Arroyo-Kalin, 2019; Walsh et al., 2019; Palmisano et al., 2021), as opposed to a measurable property of a socio-ecological system like diversity, richness, or stability (Cantarello et al., 2017; Van Meerbeek et al., 2021), including in the Neotropics (De Souza et al., 2019; Ebert et al., 2019; Arroyo-Kalin and Riris, 2021). As the systematisation of archaeological knowledge been suggested to be key prerequisite for building coherent large-scale models of feedbacks between culture and climate (Kintigh et al., 2014; Perreault, 2019), we suggest there is a need to redefine resilience analytically and thereby overcome some of these limitations. At present, there are clear obstacles to this goal that are inherent to interpretative approaches.

To do so, we draw on approaches to the study of resilience from quantitative ecology, which has some noteworthy differences from the use of Resilience Theory in archaeology as described above. Although RT as high-level theory (sensu Gunderson and Holling, 2002; Folke, 2006) is similarly prominent in the ecological literature (Newton, 2021; Van Meerbeek et al., 2021), one area of vigorous research that has gone largely unappreciated in archaeological circles is its quantification, measurement, and comparison between different cases and in response to different (or common) disturbances. This drive to develop transferable metrics of system resilience (Orwin and Wardle, 2004; Pimm et al., 2019) does not have an identifiable counterpart in our field and stands in contrast to its almost exclusive use as an interpretative heuristic to frame archaeological problems (Redman, 2005; Bradtmöller et al., 2017; Table 1). Robust debate surrounds the definition of resilience metrics and the different elements of a socio-ecological system they purport to measure (Ingrisch and Bahn, 2018; Pimm et al., 2019; Van Meerbeek et al., 2021), however, it is clear that building understandings of resilience from concise operational definitions has afforded key insights into long-term biodiversity trends, informed conservation agendas, aided in the communication of scientific results, and provided a platform for ecologists to bridge the gap with stakeholders (Nimmo et al., 2015). We follow others in contending that research on prehistoric societies has a contribution to make to all of these research areas or closely analogous ones (Kintigh et al., 2014). While there is clearly an appetite on the part of archaeologists to adopt terminology across disciplines, the field has to date largely avoided the formalism necessary to obtain meaningfully comparable data on resilience (Olsson et al., 2015) that would enable such impacts to be made.

TABLE 1

Table 1. Summary of fitted Bayesian models.

We suggest that radiocarbon-based estimates of palaeodemographic dynamics are exceptionally well-positioned to provide solutions, as they provide an absolute chronological framework for assessing fluctuations that can be cross-referenced to known perturbations and couched within a hypothesis-testing framework. As noted, climate change often figures as an important, high-level driver of prehistoric population turnover and cultural transitions, whether as a source of unpredictability (high variance), extreme (pulse) events, gradual pressure, or other mechanism (Bevan et al., 2017; Blockley et al., 2018; De Souza et al., 2019; Ebert et al., 2019; Riris and Arroyo-Kalin, 2019; Walsh et al., 2019; Palmisano et al., 2021). We now investigate whether the transient, centennial-scale effects of the MCA, in the form of a weakened South American Summer Monsoon, are reflected in persistent changes to palaeodemographic growth rates (expressed as a percent rate of change) across the Atlantic Forest, as in parts of Amazonia (De Souza et al., 2019).

Materials and Methods

Archaeological Radiocarbon Dates and Hydroclimate Proxy

For our analysis, we use 1,347 georeferenced radiocarbon dates from the Atlantic Forest, which are part of ongoing efforts to compile all published anthropogenic ¹⁴C dates in lowland South America (De Souza, 2021). The South American Atlantic Forest is a megadiverse biome that spans approximately 25 degrees of latitude, from the northeast to the far south of Brazil, and as far inland as Paraguay and northern Argentina (Marques et al., 2021). Changes in the historical extent of these forests is typically compared to its expansion or contraction relative to more open biomes, notably grasslands and arid scrub, a pattern which is historically connected to variation in long-term climate patterns across eastern Brazil (Costa et al., 2018; Barros Silveira et al., 2019). Due to known north-south variability across Atlantic Forest ecoregions, and variation in the MCA across a north-south gradient (Novello et al., 2018), we anticipate that palaeodemographic fluctuations between 1,000–700 cal BP may be similarly spatially structured and/or asynchronous. In other words, if hydroclimatic change during our study period impacted environments enough to affect population dynamics, then we expect the extent and composition of forests to be a contributing factor to Indigenous responses. To capture potential differences in resistance-resilience between different ecoregions concurrently, we group the radiocarbon dates according to the phylogeographic regions identified by Ledru et al. (2016) on the basis of topographic, climatic, and edaphic criteria. In broad terms, these regions correspond to mixed evergreen forests in the southern area including Araucaria forests, semi-deciduous forest in the central area, and dry/open forest in the northern area (Marques et al., 2021). The proposed division also corresponds to the observed pattern of variation in the MCA signal in speleothem records (northeast, core and southwest of the SACZ) (Novello et al., 2018). To account for uncertainty regarding the historic distribution of the Atlantic Forest and whether dated archaeological sites fall fully within its modern boundaries in the strictest sense, we have chosen to include dates that are within 100 km of its consensus extent, using a spatial dataset produced by Global Forest Watch in collaboration with the Brazilian Ministry of Environment and the Brazilian Institute of Geography and Statistics. The database includes all anthropogenic dates since the Late Pleistocene peopling of South America across 543 sites in our study area, as well as dates that post-date the Conquest (Figure 1, top).

FIGURE 1

Figure 1. Archaeological and palaeoclimatic overview. Top: locations of radiocarbon dated archaeological sites (grey points) and Botuverá cave speleothem (yellow triangle in the Atlantic Forest (green). Bottom: Botuverá hydroclimatic proxy (Sr/Ca: black line) and smoothed fit (orange line) for 2,000–400 cal BP, with record mean and standard deviations. MCA period highlighted by shading.

We supplement this data with a high-resolution (sub-annual) geochemical record of SASM variability from Botuverá cave in Santa Catarina state, southern Brazil (Bernal et al., 2016; Novello et al., 2018). To aid the visualisation of this data and the period of interest in the overall context of Late Holocene hydroclimate across eastern Brazil, we fit a Generalised Additive Model (GAM) to the data using a thin-plate regression spline with the Restricted Maximum Likelihood method in the mgcv R package (Wood, 2011). An advantage of GAMs for trend summaries over loess regression is the automatic estimation of penalties and optimisation of smoothing parameter to reduce overfitting (Simpson, 2018). We present this fit for the period 2,000–400 cal BP alongside the record mean and standard deviations (Figure 1, bottom), indicating the heightened aridity apparent in the southern end of the Atlantic Forest during the MCA.

Bayesian Palaeodemographic Modelling

To provide an overarching context for our study of the Atlantic Forest ¹⁴C record, we analyse the entirety of the radiocarbon dates in our sample via summed probability distributions of calibrated radiocarbon dates (SPDs).

Our previous work in the eastern La Plata basin employed information criterion-based model comparison and Monte Carlo simulation to carry out null hypothesis testing of a fitted multipart palaeodemographic model (De Souza and Riris, 2021). While this approach succeeded in identifying a demographic transition to an exponential growth phase starting after ∼1,900 cal BP that lasted until European Conquest, conventional modelling of radiocarbon data such as this has well-known limitations. Most relevant for current purposes is the conflation of unaccounted-for sampling and calibration effects (Carleton and Groucutt, 2021; Timpson et al., 2021). Here, we adopt a Bayesian inferential framework (Crema and Shoda, 2021) to account for these issues, validate the overall trajectory of our Atlantic Forest SPD during this interval based on this prior knowledge (De Souza and Riris, 2021), and derive posterior samples of relevant demographic parameters.

Guided primarily by this preceding work (De Souza and Riris, 2021), we define a bounded exponential growth model for our Late Holocene radiocarbon dataset in the R package nimbleCarbon v0.1.0 (Crema, 2021). This package builds on the functionality of the nimble package (de Valpine et al., 2021), specifically using a Markov Chain Monte Carlo procedure (Metropolis-Hastings random walk algorithm with adaptive sampler) to calculate the probability mass of the SPD at each timestep and obtain posterior samples of model parameters. In this case, our parameters are the start (a) and end (b) points of our period of interest (1,900–500 cal BP) and the growth rate r with an exponentially distributed prior of λ = 1/0.0004 (Crema and Shoda, 2021). In the context of the archaeology of the Atlantic Forest, this period corresponds approximately to the transition between the shellmound tradition (sambaquis) and the expansion of major archaeological ceramic traditions in eastern Brazil, namely, the Tupiguarani, Una/Taquara-Itararé (Southern Jê), and Aratu traditions (Bonomo et al., 2015; Iriarte et al., 2017; De Souza et al., 2020). Because these traditions are generally agreed to have spread or intensified polyculture agroforestry practices, we expect growth rates to be broadly in line with other agrarian or semi-agrarian societies worldwide (Zahid et al., 2016). Following Arroyo-Kalin and Riris (2021) we also consider the possibility that populations had stabilised by or before Conquest, which we model via a logistic growth model with an additional carrying capacity parameter k, whose prior is set as an exponential distribution truncated between 0.001 and 0.2 after Crema and Shoda (2021). For both models we ran two chains of 30,000 MCMC iterations each, with 1,000 burn-in steps and a thinning interval of three. We use the inbuilt diagnostics to check for convergence and adequate mixing. Posterior predictive checks were achieved in an analogous fashion to the more common regression-based methods in rcarbon (Crema and Bevan, 2021) we have previously used (De Souza and Riris, 2021), namely by: repeatedly sampling random calendar dates from fitted models based on the posterior parameters, “back-calibrating” into ¹⁴C dates with randomly assigned errors, calibrating the datasets, and aggregating them into a simulated SPD envelope based on 500 iterations of this procedure (Timpson et al., 2014). We have not pruned our dataset, but bin radiocarbon dates from the same site that lie within 200 years of one another. The analysis is fully reproducible with the code provided in Supplementary Material.

A Test for Resistance-Resilience

In addition to the broad-scale perspective afforded by the Bayesian palaeodemographic modelling, we seek to recover cross-sectional differences in human responses to prehistoric climate change in the Atlantic Forest. To this end, we present a novel post-hoc significance test that adapts the functionality of “point-to-point” testing (Edinborough et al., 2017) to the established mark permutation test (Crema et al., 2016). The latter advantageously does not suffer from the same issues surrounding growth model definition as regression-based null hypothesis testing, while the former class of tests has demonstrable utility for evaluating the presence and significance of shorter-term (centennial-scale) “events” in SPD timeseries between two defined points in time (p₁ and p₂). Code to replicate our post-hoc test is provided in full in Supplementary Material. Per standard protocol in the aggregate analysis of radiocarbon dates, we bin dates originating from the same site that fall within 200 years of one another to minimise the overrepresentation of well-dated sites in our time series. We do not prune our dates. For the permutation tests, we calibrate dates using the using the ShCal20 curve (Hogg et al., 2020), except for those obtained on marine samples, in which cases we used the Marine20 curve (Heaton et al., 2020) and average ΔR offsets and errors based on the 10 nearest published values.¹ SPDs for the permutation tests are left unnormalised to minimise calibration curve artefacts in their shapes (Bevan et al., 2017).

For our test, we adapt ecological indices for resistance and resilience to the analysis of SPDs (Orwin and Wardle, 2004; Nimmo et al., 2015). Combining these system properties in a bivariate analysis enables the intuitive comparison of the behaviour of a proxy during a period of interest over a geographical region (Figure 2). In straightforward terms, resistance can be considered the capacity of a population to absorb and maintain its state in the face of disturbances. Analytically, resistance is defined as the deviation of a proxy from a baseline following a perturbation, standardised relative to the baseline (Orwin and Wardle, 2004):

FIGURE 2

Figure 2. Schematic outline of analytical approach. Cross-sectional data corresponding to the radiocarbon records of different geographical regions are rendered (via summed probability distributions) into longitudinal data. Measurements before (p1), after (p2), and during (p3) a disturbance D are taken from the palaeodemographic proxy and used to estimate resistance and resilience. Metrics provide an empirical basis for comparing spatiotemporal trends.

R e s i s t a n c e = 1 - \frac{2 \times | G_{b} - G_{x} |}{G_{b} + | G_{b} - G_{x} |} (1)

where G_b is the SPD rate of change at the baseline p₁, taken to be the start of a perturbation, and G_x is the index value after a defined period of time. Both metrics are bounded between −1 and 1, however, if using Equation 1 as-is, negative values of G_b (which are common in SPD-derived rate of change estimates) can lead to very small divisors when the difference between G_x and G_b is also small. In effect, this may cause Equation 1 to return values in excess of ±1. A simple and effective solution, which we advocate here, is to take the absolute value of G_bin the divisor too:

R e s i s t a n c e = 1 - \frac{2 \times | G_{b} - G_{x} |}{| G_{b} | + | G_{b} - G_{x} |} (2)

While any point in time may be chosen for G_x, the maximum deviation from baseline conditions is recommended for pulse perturbations, such as abrupt climate change (Van Meerbeek et al., 2021, p. 10). We therefore set this point (p₃) as the minimum rate of change experienced between p₁ and p₂. It should be underscored that our code (Supplementary Material) provides flexibility for any point in time between p₁ and p₂ to be chosen.

Resilience is defined as the extent to which the proxy is able to return to the reference condition following disturbance, standardised by the total relative change:

R e s i l i e n c e = \frac{2 \times | G_{b} - G_{x} |}{| G_{b} - G_{x} | + | G_{b} - G_{e} |} - 1 (3)

Here, G_b and G_x are as in Equation 1, while G_e is the condition of the proxy at the end of the perturbation (p₂).

A resistance value of 1 shows no effective change (maximal resistance to perturbations), zero resistance indicates a relative change of 100% of the baseline value, while negative values indicate the extent to which growth rates decline over the period of interest. For resilience, a value of 1 indicates full recovery at the end of a perturbation, 0 indicates no recovery from maximum impact. Negative values of resilience are returned when |G_b−G_e| > |G_b−G_x|, which can only occur when G_x is manually set at a point other than the maximum deviation (minimum growth rate) and is therefore not relevant here.

The expected ranges of Equations 2 and 3 are shown in Figure 3 for different parameter combinations. It is worth noting that because of their double exponential shapes, resistance and resilience are equal to zero for two different values of G_x: when G_x is itself equal to 0 and when G_x is two times G_b (Figures 3A,B). This means that when G_b is positive, zero resistance (100% change) is returned when the growth rate doubles (to 0.04, here). Similarly, when G_b is negative, zero resistance is returned when the growth rate decreases by the same proportion. The same is true of resilience when varying the endpoint (Figure 3C): if the start- and end-values (G_b and G_e) are the same, the formula yields 1, indicating that the system in question is fully resilient and successfully bounces back to its original condition. However, the equation does not distinguish directionality in the proxy. In other words, and as in the case of resistance, if a system bounces back and exceeds the original growth rate, the formula will return a value <1. This is counterintuitive to how archaeologists typically interpret resilience, as noted above, and demonstrates that these features of the equations must be kept in mind. Change to a system, even if “beneficial,” cannot necessarily be regarded as an instance of resilience. We return to this point in the discussion.

FIGURE 3

Figure 3. Resistance and resilience as functions of G_x and G_e. G_b is held constant at 0.02 (vertical dashed lines). (A) Resistance is equal to zero for two different values of G_x, when G_x = 0 and when G_x is two times G_b. (B) Resilience as a function of G_x. Note that different values of G_x may return the same Resilience index. (C) Resilience after Orwin and Wardle (2004), with a varying end point (G_e), and constant start (G_b) and minimum (G_x). Resilience equals 1 only when G_e = G_b. Lower and higher values of G_e yield Resilience <1.

Finally, we take advantage of the inbuilt functionality of the mark permutation test to compare the empirical resistance-resilience metrics with those obtained from envelopes of simulated SPDs. In effect, the permutation test generates a distribution of expected values of resistance and resilience for each region. Our analyses consist of determining whether and how observed values deviate from these expected values, which we highlight with a two-tailed test of significance of whether a given metric is significantly above or below the global aggregate. This is distinct from the p-value obtained through the mark permutation test itself, which we also report separately (Edinborough et al., 2017). We apply our point-to-point test for resistance-resilience to the “rate-of-change” version of the mark permutation test, as this allows easier comparison with the growth rate parameter (r) estimated through the Bayesian model fitting described above. Rather than binary states of resilience/vulnerability, our test indicates the relative location of empirical patterns along a bounded spectrum. To guard against our results being conditional solely on the choice of start- and end-dates, bearing in mind the likely existence of lead-lag relationships between our growth rate proxy and hydroclimate, we also perform sensitivity analyses with the point-to-point test. As the chosen frame of reference for the MCA is 300 years (1,000–700 cal BP), we also consider start- and end-dates within 30 years (10% of total time span) on either side of these dates in 1-year increments (1,030–670 cal BP) and report the ranges of the indices for these parameter combinations. Here we use 9,999 random mark permutations, with a running mean and backsight of 100 years. We recommend that the backsight, which is effectively an arbitrary parameter choice, be scaled to the duration of the phenomena under study. In this case, we believe a priori that centennial-scale change would make this choice acceptable. In principle, further sensitivity tests could also take changes in the backsight value into account as with the choice of date ranges.

Results

Bayesian Modelling and Parameter Estimation

The results of the Bayesian model fitting and diagnostic checks (Table 1 and Supplementary Figures 1–3) indicate adequate mixing and a sufficient effective sample size for both models. In good agreement with prior work on a subset of this radiocarbon dataset, the pre-Columbian population trajectory of Atlantic Forest conforms well to both exponential and logistic growth models between 1,900 cal BP and the start of Portuguese colonisation after ∼500 cal BP. In comparative perspective, we find that the exponential model performs better than the logistic (Table 2). Considering also that the k parameter converged to a lesser degree in the logistic model, we proceed with the more parsimonious exponential fit for subsequent posterior predictive checks, although in reality differences are minor (Supplementary Figures 2, 3). The 95% confidence envelope obtained from 500 simulated SPDs conditioned on posterior samples of the fitted model (Figure 4) shows that the empirical SPD is generally well within the predicted range, excepting close to the end of the analytical time range, when radiocarbon dates become sparser.

TABLE 2

Table 2. Model comparison summary.

FIGURE 4

Figure 4. Results of Bayesian radiocarbon analysis. Left: results of posterior predictive check, showing the empirical Atlantic Forest SPD (black line) within the prediction range of the fitted model (grey envelope), excepting close to the Conquest. Bottom right: marginal posterior distribution of the growth rate r over the interval 1,900–500 cal BP. Top right: exponential model fit.

Additionally, the results highlight the 95% range of the marginal posterior distribution of the exponential growth rate between ∼0.12 and 0.165% (Table 1), which is high relative to long-term estimates obtained elsewhere (Zahid et al., 2016) but approximately in line with recent palaeodemographic work (Crema and Shoda, 2021). The results may be, in part, an artefact of a constrained frame of reference, during a known period of significant expansion of pre-Columbian farming cultures (De Souza et al., 2020; De Souza and Riris, 2021). The Bayesian modelling does not directly impact the results of the resistance-resilience tests, but rather, sustains the findings of previous work in eastern Brazil, while also accounting appropriately for sampling error and calibration effects. These results inform our subsequent tests and invest them with the knowledge that, in aggregate, pre-Columbian cultures of the Atlantic Forest were undergoing largely uninterrupted exponential growth during the Common Era. Short-term perturbations, such as the MCA, must therefore be seen in the context of this overall trajectory.

Resistance-Resilience

The resistance-resilience test highlights key differences between the southern, central, and northern regions of the Atlantic Forest. Our sensitivity testing returned a relatively constrained range of values for both metrics within 30 years of the start and end of the MCA (Figure 5, bottom left), indicating that the choice of start date either side of 1,000 and 700 cal BP does not qualitatively alter the results. The obtained values for resistance across all three regions are statistically significant from their counterparts between 1,000–700 cal BP, while resilience returns no results that are distinguishable from the null hypothesis of no regional differences (Table 3). Before delving into the regional patterns, it is worth underscoring that the regions all have the same direction of change during the MCA. That is, the effects of the MCA disturbance on growth rates appear to be qualitatively similar, with generally negative trends observed, but differ in magnitude, timing, and ultimate impact. This underscores the importance of considering more than one index and grouping of dates concurrently; single metrics or settings are not as informative.

FIGURE 5

Figure 5. Results of point-to-point permutation test for resistance-resilience. Right column: Rate of change of SPDs in Northern, Central, and Southern subsets of the Atlantic Forest, with p-values for the underlying permutation test. Simulation envelope is based on 1,000 random mark permutations. Orange dots and shading show the defined start- and end-points of the MCA, purple dots are the minima of the time series between these points. Note that with a 100-year backsight, the time series starts at 1,900 cal BP. Bottom left: comparisons of resistance-resilience between regions.

TABLE 3

Table 3. Resilience-resistance metrics and significance test for Northern, Central, and Southern subsets of Atlantic Forest radiocarbon record, and time in years between start of disturbance and time series minima (lag).

The Central and Southern regions have resistance values on the same order of magnitude, reflecting the stagnation of growth rates close to zero in the wake of the MCA in both these settings (Figure 5). The timing of this phenomenon unfolds over two centuries in the Central Atlantic Forest, while in the Southern region, analogous impacts occur in less than a century. This implies that the “on-the-ground” experience of climate impacts may have been more acute in the southern Atlantic Forest compared to the centre and north. We return to the extremely low resistance in the Northern subset in the discussion, noting that it reflects a drop from an essentially stationary population to one contracting at a rate of −0.66% during the MCA.

We are cautious about overinterpreting the metrics for resilience, bearing in mind that there are no significant regional differences (Table 3). Qualitatively, all three rate of change time series have visual similarities during the 1,000–700 cal BP interval, i.e., a limited or non-existent return to pre-disturbance rates of change. Somewhat lower values for resilience in the Central and Northern areas may be related to their comparable lags of 20 and 83 years, respectively, as opposed to 209 years during which to recover in the south. In other words, the more abrupt, yet shorter, MCA impacts on populations in the south may have translated into a longer and more sustained recovery period. The Northern region stands out as the least resistant and least resilient region, with growth rates at 1,000 cal BP, as noted, close to zero that fall well below zero by 700 cal BP. It is, moreover, the only region to undergo a statistically significant negative deviation in the mark permutation test during the MCA, further marking out the narrow coastal strip of Atlantic Forest as unusual in the context of this study. The significant global p-value in the Central region is almost certainly largely because of the large negative deviation around ∼500 cal BP. The coastal strip here is coincidentally also the location of the first Portuguese colonies in Brazil, but we cannot conclusively attribute this pattern to the Conquest.

In summary, the results reveal low resistance to abrupt hydroclimatic change among our subdivisions of the Atlantic Forest following the start of the MCA, although apparently differing in abruptness of impact. There is also a consistently low ability to recover by the end of the MCA period, with resilience being less than half of its pre-disturbance value (<0.5) across all regions. Despite apparent commonalities and lack of statistically significant deviation, it is worth noting that by 700 cal BP, Northern growth rates are negative, Central growth is virtually stationary, and Southern rates are well above zero. As noted above, the interface between the MCA and LIA was a period of continuing climate change across tropical South America (Novello et al., 2018; Azevedo et al., 2019), which may be key to contextualising these qualitative differences, although we presently lack the resolution and sample size of radiocarbon dates to directly investigate further so close to European Conquest (Riris, 2019).

Our study of Atlantic Forest palaeodemography reveals correlations between abrupt climate change and key parameters, with a degree of spatiotemporal variability that suggests not all areas were equally affected. The results of our analysis introduce novel perspectives on current knowledge about Late Holocene palaeodemography in South America (Azevedo et al., 2019; Arroyo-Kalin and Riris, 2021; De Souza and Riris, 2021), complementing the findings of previous top-down demographic model-fitting on a broad scale with more focused interrogation of local records. It is plausible that impacts may have been differentially mediated in one or more ways. We now turn to exploring the implications of these findings for our case study specifically, as well as for the study of prehistoric resilience more broadly.

Discussion

Our findings affirm that a delayed demographic transition occurs in the Atlantic Forest relative to the Amazon basin following the adoption and development of sedentary agroforestry (De Souza and Riris, 2021). In contrast, in the Amazon, centres of local domestication and early adoption of exotic crops date to the early to mid-Holocene (Watling et al., 2018; Lombardo et al., 2020). Due to this, it is possible that the ability to reach carrying capacity would have been attained earlier in the Amazon than in the Atlantic Forest. Although we detected minimal differences between the logistic and exponential model fit (Supplementary Figures 2, 3), an exponential model of growth, besides being more parsimonious, is compatible with the cultural trajectories of the Atlantic Forest against those of the Amazon over the period of interest (Arroyo-Kalin and Riris, 2021). Before continuing, we underscore the persistence of Indigenous populations over our period of study. In the context of the last 1,500 years before the Conquest, Indigenous people were still experiencing unabated exponential growth (Figure 4). The relatively late occupation of the Atlantic Forest by farming populations must, however, be considered within its changing environmental context. The best example is the dispersal of the Tupi-Guarani, extremely rapid to the point where determining the precise chronology and direction of spread becomes unfeasible (Riris and Silva, 2021). Events of accelerated demic expansion, perhaps due to the migrant societies finding a new or relatively unexplored niche, may have contributed to the observed exponential growth, suggesting the Atlantic Forest societies were far from reaching carrying capacity at the eve of the Conquest. Furthermore, the role of late Holocene climate-driven forest expansion in the Atlantic Forest biome has to be considered (Smith and Mayle, 2017). As in other global cases of the expansion of agriculturists (Russell et al., 2014), the increase in the areas that could be exploited by polyculture agroforestry is likely to have facilitated the Tupi-Guarani dispersal (Iriarte et al., 2017). The expansion and contraction of certain niches would have influenced the dispersal and growth not only of the Tupi-Guarani, but also of Jê societies, whose subsistence relied on Araucaria forests in southern Brazil (Iriarte and Behling, 2007; Pereira Cruz et al., 2020).

Bearing these points in mind, we note that low resistance metrics for the Northern region, coupled with comparatively higher resilience observed in the centre and south, appear at first to correlate with the latitudinal gradient in the signal of the MCA observed in eastern Brazilian speleothems, which show a stronger deviation north of the SACZ (Novello et al., 2018). The speleothems at the core region of the SACZ, corresponding to our central sites, exhibit increased climate variability at the MCA-LIA transition (Novello et al., 2018), which may explain the decline in growth rates at the end of the MCA period (Figure 5). The northeast region of Brazil is expected to be the most sensitive to hydrological changes; climate simulations have confirmed its susceptibility to biome shifts when compared to the southern and central parts of the Atlantic Forest domain (Ledru et al., 2016; Costa et al., 2018). Relative biome (in)stability may therefore be an important factor in mediating the resistance and/or resilience of pre-Columbian populations to abrupt change, to a greater degree than absolute states such as “open” or “closed” vegetation (alternatively, forest or grassland). Quantifying tropical ecosystem stability—and hence the viability of different subsistence strategies in terms of resistance-resilience—against the backdrop of palaeodemographic dynamics is likely a useful avenue for future work. The main factor behind the observed north-south trend in resistance may be hydroclimatic regime variability, as the southern region receives precipitation from southerly, extra-tropical, oceanic sources in winter (Cruz et al., 2005; Campos et al., 2019). This is in agreement with the observation that, in the early to mid-Holocene transition (∼8,200 cal BP), the impacts of precipitation variability were lessened in more climatically stable regions, such as the southern part of South America (Riris and Arroyo-Kalin, 2019). Climate therefore seems a likely candidate for the main driver of the observed differences in resistance-resilience in our study area. Consequently, we reject the null hypothesis of no change to palaeodemographic growth rates across the Atlantic Forest during the MCA.

The different responses of Atlantic Forest societies during the MCA provides additional context for the expectations posed by De Souza et al. (2019), who argued in central-eastern Amazonia that polyculture agroforestry should provide a more resilient economic basis in the face of climate change than specialised, intensive land use systems. By separating palaeodemographic responses into the components of resistance and resilience, we are able to qualify this inference; the capacity to rebound following abrupt disturbances appear to be common across the agroforestry systems under study (resilience p > 0.05 in all cases), but responses to initial impacts (notably its timing and relative magnitude) appears to depend on local circumstances to a greater degree. In a topographically and biogeographically diverse environment, such as the Atlantic Forest, this is a relatively intuitive finding. Our findings parallel ecological resilience studies, where multiple metrics are found to be more capable of revealing different facets of the systems under study than single indices (Cantarello et al., 2017). Comparative studies will likely help to add further nuance, as well as better understandings of the components (species, cultivars, land use systems) that made up different forms of polyculture agroforestry in the past, which certainly had considerable internal variation. To this end, in addition to climatic/vegetation variability within the Atlantic Forest range, socio-cultural factors could be behind some of the observed differences in resistance-resilience. The available data suggest that the societies throughout the studied area practised similar forms of polyculture agroforestry, albeit in different niches (De Souza et al., 2020; Pereira Cruz et al., 2020), yet multiple distinct forms of social organisation existed prior to Conquest. Recent studies have emphasised pre-Columbian monumentality in southern Brazil as an indicator of territoriality and inter-group competition (De Souza et al., 2016), however, the emergence and maintenance of corporate architectural traditions can also be viewed as evidence and facilitators of strengthened intra-group cohesion (Iriarte et al., 2013). A possible consequence of more explicitly manifested and formalised social structures includes an ability to coordinate and manage resources above the kin level. Although we have focused on MCA climate as a driver of palaeodemographic impacts, we suggest that future work should consider quantitative approaches towards social and environmental factors as mediators of these impacts, and thus, investigate their possible role as facilitators of demographic and societal resilience.

We have focused here on the potentially deleterious effects of abrupt hydroclimate variability on population growth, highlighted previously elsewhere in lowland South America (Azevedo et al., 2019; Riris and Arroyo-Kalin, 2019). Our approach expects that the system(s) under study display resistance and resilience if baseline conditions are maintained or are returned to, respectively. There are, however, other scenarios for which our tests may be suitable, including inverse to those presented here. Returning to the indices as illustrated in Figure 3, one can also imagine “perturbations” to cultural systems such as the invention and spread of a technological innovation that raises growth rates over an interval of time. This hypothetical scenario would also result in values of <1 being returned in our metrics. Such a result is not because a given cultural system or group was (necessarily) adversely affected by rising growth rates, but because the nature of the system has measurably changed by the adoption of an innovation, reflected in its “low” resilience. The equations, and by extension our tests, are agnostic about whether change is a net positive or a net negative. We feel it is important to underscore this point, as it is contrary to commonplace understandings of resilience in the archaeological literature, which tends to interpret signs of “improvement” (i.e., increasing growth rates) as “high” resilience (Bradtmöller et al., 2017). Another, perhaps more intuitive, way to view resilience in this sense is that the prior cultural system was vulnerable to the introduction of an innovation, and so “collapsed” into a novel regime incorporating it. As always, caution must be exercised when deploying metrics; their broader societal and, in this case, palaeodemographic contexts are important for enabling considered interpretations of their meanings. The equations and indices make no distinction about the significance of low or high resistance-resilience values—careful consideration of the cultural system as a whole must accompany interpretation. Our aim is to encourage careful use in future applications, as the benefits of systematising the archaeological study of societal resilience are potentially great: comparability, communicability of results, rigour, and facilitating interdisciplinarity and impact.

Conclusion

This paper has sought to develop novel methods for the archaeological study of resilience in prehistory, a topic of significant concern (Bradtmöller et al., 2017; Russo and Brainerd, 2021). We have, both previously and separately, argued for the systematisation of the study of resilience in palaeodemography (De Souza et al., 2019; Riris and Arroyo-Kalin, 2019), in parallel to similar calls elsewhere (Bird et al., 2020). We are motivated by the need to make the archaeological study of resilience more rigorous and transparent, an essential part towards building interdisciplinary dialogue (Moser et al., 2019, p. 37). Here, we have drawn on mature research in quantitative ecology to propose a common framework for measurement of resilience in archaeology. Our method, which takes advantage of probabilistic, continuous time series data, can also be viewed as an answer to collapse-centric narratives on past societal responses to abrupt change (Middleton, 2017, p. 81; Koch et al., 2019; Degroot et al., 2021). We argue that the approach encourages attention to particular circumstances in archaeological and historical perspective, while also enabling the formalism that typifies the ecological resilience literature, and which we regard as a prerequisite to developing a comparative, deep time perspective on coupled socio-environmental changes and their consequences. We have demonstrated this by testing palaeodemographic responses to disturbances across and between biogeographically relevant zones of the South American Atlantic Forest during the Medieval Warm Period to explore the implications of hydroclimatic variability for the pre-Columbian societies of this biome. The resulting analysis of longitudinal and cross-sectional trends in palaeodemography reveal notable regional differences.

It remains hard to argue that there is a single “correct” way of viewing resilience (Middleton, 2017). Nonetheless, building a perspective that by its nature is rooted in probabilistic evidence of past population changes, itself is based on the ubiquitous use of archaeological radiocarbon dates, certainly has advantages that other proxies do not. Notably, our approach provides explicitly defined, quantitative, and reproducible alternatives to interpretative approaches, which accompany this paper for future use and refinement (Supplementary Material). We acknowledge that point estimates of resistance and/or resilience are likely imperfect reflections of a more complex picture in the past and that there will always be alternatives to the equations and tests we have implemented. Even so, we suggest that directly comparable data on prehistoric resistance-resilience provide a more useful point of entry into identifying what factors (behaviours, adaptations, institutions, ecologies) mediate the severity of disturbance impacts than incomparable ones, as well as whether trade-offs exist between different factors in terms of robustness and fragility. These issues already form the focus of substantive research (Peregrine, 2020), and the present work charts a path through them by focusing on measurable features of ancient population dynamics. Furthermore, our approach enables the archaeological study of resilience to locate commonalities of language with the ecological resilience literature.

Our method therefore paves the way for the systematic and rigorous evaluation of inferred disturbances to human societies and culture-climate links in general. We suggest that other well-known phases of abrupt climate change, for example the Younger Dryas/Bølling-Allerød, the Pleistocene/Holocene transition, the 8.2 and 4.2 ka events, and data permitting, the “Columbian Exchange” (Hamilton et al., 2021) would be ideal test beds for the approaches we advocate. While attention here has mainly been on climate as a driver, abrupt social change or reorganisation may be equally important, and interesting, targets for comparative analysis (such as periods of warfare, Edinborough et al., 2017). Archaeological palaeodemography is also uniquely positioned to offer enhanced estimates of effective population size and density, which are key to understanding the deep time earth systems impacts of humanity more comprehensively (Klein Goldewijk et al., 2017; French et al., 2021; Morrison et al., 2021). Finally, we anticipate that continuing analytical refinements in Bayesian inference (Carleton, 2021; Crema and Shoda, 2021) as well as hereto unexplored methods such as symbolic transfer entropy (Camacho et al., 2021) will open up the study of causal links between population resilience and induced disturbances.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

PR and JS: conceptualisation and writing. PR: analysis and figures. Both authors contributed to the article and approved the submitted version.

Funding

Publication was generously supported by the Bournemouth University Open Access Publication Fund. JS is supported by the MAHSA (Mapping Archaeological Heritage in South Asia) project, funded by the Arcadia Fund (UK).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Marc Vander Linden, Fiona Coward, Emma Jenkins, Kim Davies, Fabio Silva, and Adrian Newton for comments on earlier drafts of this manuscript. Three reviewers provided helpful and insightful feedback on the manuscript that improved it greatly.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.740629/full#supplementary-material

Supplementary Figure 1 | Results of Bayesian radiocarbon analysis with logistic prior. Left: results of posterior predictive check, showing the empirical Atlantic Forest SPD (black line) within the prediction range of the fitted model (grey envelope), excepting close to the Conquest. Bottom right: marginal posterior distribution of the growth rate r over the interval 1900–500 cal BP. Top right: marginal posterior distribution of the carrying capacity k over the same interval.

Supplementary Figure 2 | Exponential model trace plots and goodness-of-fit diagnostics.

Supplementary Figure 3 | Logistic model trace plots and goodness-of-fit diagnostics.

Footnotes

^ http://calib.org/marine/

References

Albuquerque, M., and Lucena, V. (1991). Agricultura tropical pré-histórica: um sistema de floresta úmida ou que integra o semi-árido? Ciência Trópico 19, 7–33.

Google Scholar

Arroyo-Kalin, M. (2012). Slash-burn-and-churn: landscape history and crop cultivation in pre-Columbian Amazonia. Quat. Int. 249, 4–18. doi: 10.1016/j.quaint.2011.08.004