Abstract

This paper presents results from the MayaSim model, an integrated agent-based, cellular automata, and network model representing the ancient Maya social-ecological system. The model represents the relationship between population growth, agricultural production, soil degradation, climate variability, primary productivity, hydrology, ecosystem services, forest succession, and the stability of trade networks. Agents representing settlements develop and expand within a spatial landscape that changes under climate variation and responds to anthropogenic impacts. The model is able to reproduce spatial patterns and timelines somewhat analogous to that of the ancient Maya, although this proof-of-concept model requires refinement and further archaeological data for calibration. This paper aims to identify candidate features of a resilient versus vulnerable social-ecological system, and employs computer simulation to explore this topic, using the ancient Maya as an example. Complex systems modelling identifies how interconnected variables behave, considering fast-moving variables such as land cover change and trade connections, meso-speed variables such as demographics and climate variability, as well as slow-moving variables such as soil degradation.

Keywords:

Social-Ecological System, Archaeology, Cellular Automata, Network Model, Trade Network, Agent-Based Model

Introduction

1.1

Few topics gain as much cross-disciplinary interest as the rise and fall of ancient civilisations. The story of development and demise in complex societies contains narratives of the human endeavour threatened by devastating droughts, greedy rulers, foreign imperialism, and overuse of natural resources, among others. Societies are, however, a set of interacting elements which as a whole express characteristic features, interpreted as emergent properties of underlying processes at multiple scales. Designing a holistic approach to understanding social-ecological systems requires methods which simultaneously observe patterns in many dimensions, a kind of observation for which van der Leeuw (2012) argues that traditional Western science is not very well equipped. An analogy is the example of solving a Rubik's Cube, in that one cannot get the cube 'in order' by dealing first with one side, then the next, and so forth. The only way to arrive at order is by looking at the patterns on all sides simultaneously, and not favouring any particular one at any time (van der Leeuw 2012). This paper presents a method to identify candidate features of a resilient versus vulnerable social-ecological system, and employs complex systems science, using computer simulation to explore this topic using the ancient Maya as an example.

1.2

A number of research questions are presented for exploration:

• What dynamics lead to the development of the densely populated and interconnected human geography of the ancient Maya?
• Is it possible to use computational social science to 'grow' the three Maya temporal periods of the Preclassic (1000 BC - AD 250), Classic (AD 250 - 900), and Postclassic (AD 900 - 1500)?
• How does the simulated social-ecological system develop and respond to changing conditions, and what modelled indicators warn of vulnerability?

In order to explore these research questions, a simulation model was designed and calibrated for the landscape of Central America. Model runs produce temporal and spatial patterns that can be understood through examining the underlying assumptions of the different integrated components of the model. MayaSim is a combined agent-based, cellular automata, and network model that represents the ancient Maya social-ecological system. Agents, cells, and networks are programmed to represent elements of the historical Maya civilisation, including demographics, trade, agriculture, soil degradation, provision of ecosystem services, climate variability, hydrology, primary productivity, and forest succession. Simulating these in combination allows patterns to emerge at the landscape level, effectively growing the social-ecological system from the bottom up. This approach constructs an artificial social-ecological laboratory where different theories can be tested and hypotheses proposed for how the system will perform under different configurations.

1.3

The model is able to reproduce spatial patterns and timelines somewhat analogous to that of the ancient Maya's history. This proof of concept model requires refinement and further archaeological data for calibration to improve results, although it is noted that there is little empirical evidence by which to validate such models, and such evidence is generally site-specific and discontinuous through time.

1.4

The purpose of the model is to better understand the complex dynamics of social-ecological systems and to test quantitative indicators of resilience as predictors of system sustainability. An integrated agent-based, cellular automata, and network model was constructed using the software Netlogo (Wilensky 1999). The full model, code and documentation is available in Heckbert (2012) via the www.openabm.org website, and further description of the model in the context of Maya archaeological literature is presented in Heckbert et al. (in press).

Methods

2.1

The MayaSim model represents settlements as agents located in a gridded landscape. The model is constructed using the software Netlogo (Wilensky 1999). The software interface, shown in Figure 1, presents the spatial view of the model with graphs tracking model data and a user interface for interacting with the model. The view can be changed to observe different spatial data layers within the model. The model operates at a spatial extent of 516,484 km² with a 20 km² resolution. Temporal extent is 650 times steps, each representing roughly 2 years.

2.2

Upon model initialisation, base GIS layers are imported using the Netlogo GIS extension. Static cell variables are set therein, dynamic variables are reset to default values set on the user interface, and settlement agents are randomly initialised in the spatial landscape. Imported spatial data include elevation and slope (Farr et al. 2007), soil productivity (FAO 2007), and temperature and precipitation (Hijmans et al. 2005). Data for soil productivity was slightly smoothed using the Netlogo function diffuse to make differences along polygon boundaries less stark, but all other data is unchanged from the cited sources. Data is resampled at a 20km² resolution using the Netlogo GIS extension (see Wilensky 1999). All GIS data is high resolution and deemed to be of high quality except possibly the soil productivity dataset. The latter dataset can be replaced with finer resolution regional or country-specific data, but because the region extends across several modern countries with inconsistent data sets, the best globally consistent dataset available was used. Future work will aim to replace the existing dataset with the Harmonized World Soil Database (FAO et al. 2009) for finer spatial resolution and consistency across country borders.

Figure 1. MayaSim model interface with interactive controls, spatial view, and graphs tracking model data. Agents operate on a cellular landscape and are connected by links within a network.

2.3

The simulation begins with calculations of biophysical variables for water flow and net primary productivity, and these are further used to calculate forest succession, agricultural production, and ecosystem services. Settlement agents interact with the spatial landscape to generate agricultural yield through cropping, derive benefit from local ecosystem services, and generate trade benefits within their local trade network. The combined benefits of agriculture, ecosystem services, and trade drives demographic growth including migration. Simulating the integrated system reveals how the social-ecological system functions through time. Additional agents include a 'migrant' agent who settle new locations, a 'raindrop' agent which routes hydrological surface flow, and a 'network search agent' who traverses links between connected settlement nodes to calculate network statistics. The following section describes biophysical and anthropogenic processes that are tracked in the model.

Biophysical functions: climate, hydrology, soil, primary productivity, forests and ecosystem services

2.4

Spatial data for precipitation and temperature (Hijmans et al. 2005) representing current conditions (1950 to present day) is adjusted within the model using assumptions that mimic the paleoclimatic variation presented in Prufer et al. (2011). Variation in precipitation is more pronounced towards the northwest of the case study area. This diagonal north-south / east-west variation in rainfall is created over the current-day GIS data for precipitation, using the function:

(1)

where R_j,T is precipitation [mm] for cell j at initial time step T, and CL_n is a localized rainfall effect due to the presence of cleared land on neighbouring cells n = 1 … 8 with weighting parameter δ determining the strength of this effect. DF_j is the distance [km] of each cell from the top northwest corner of the map and is the furthest distanced cell from this point. RC_t cycles from + 20% to -10% linearly over a 56 time step cycle, and t = 1 … 650. This function serves to reduce and increase rainfall cyclically, with a more pronounced effect further towards the northwest.

2.5

A surface-flow hydrological model was devised to recreate past conditions based on elevation and rainfall. These data are used to calculate surface flow and location of potential seasonal standing water, consistent with Reaney (2008). Each cell generates a mobile 'raindrop' agent containing precipitation volume R_j,T. The raindrops follow the elevation data, repeatedly moving to the adjacent cell with the lowest elevation (and considering the summed volume [mm] of raindrops already at that location). If raindrops cannot move (i.e., a location is flooded), the raindrops 'pool' and form river and lake patterns, as shown in Figure 2 depicting simulated water flow under climate variation assumptions. This flow pattern can be validated against hydrological processing functions in GIS software. The function serves to move water based on elevation, and can generate the spatial distribution and surface water flow as precipitation varies across the climate cycle. See equations in model code readily explored in Heckbert (2012). The result of these equations is a finely-resolved spatial pattern that is able to recreate the potential river, lake, and wetland systems under climate variation assumptions. An important caveat of the model assumptions, however, is the lack of heterogeneity in infiltration rates across the region. In reality, the surface water in many areas drains rapidly to the aquifer owing to the karst limestone bedrock, however spatial data to inform this process was not able to be sourced, and as such a constant rate of infiltration is assumed across the landscape, and is a topic for further research effort.

Figure 2. Simulated surface water flow [m3] (blue on white background, darker blue being greater flow) based on climate variation assumptions of a) -15%, b) 0%, and c) +15% change in annual precipitation in equation 1.

2.6

The base GIS layers for rainfall and temperature are used to calculate net primary productivity, which in turn is used to calculate forest succession. GIS soil data is used to calculate agricultural productivity, and all these combine to calculate provision of ecosystem services. While we recognize that these relationships are complex, simplifying assumptions allow initial representations as methodologies are further refined.

2.7

Forest succession operates as a cellular automata model, where the state of a cell is dependent on internal conditions and is influenced by the condition of neighbouring cells. In this proof-of-concept stage, cells take on one of three general forest states that represent climax forest, secondary regrowth, and cleared/cropped land, referred to as state 1, 2 and 3 respectively. The forest state is decremented for 3.5% of randomly selected cells, to represent natural disturbance. The disturbance rate is linearly amplified by population density of nearby settlements to represent local wood harvesting, to a maximum of 15%. Cells advance in their forest state based on the time since last disturbance and the relative net primary productivity of the cell. Once the time since last disturbance is above a threshold (40 * NPP_{max j, t} / NPP_{j, t} years for secondary regrowth and 100 * NPP_{max j, t} / NPP_{j, t} years for climax forest, to account for spatial variation in net primary productivity) the forest converts to the new state. For conversion to climax forest, a cellular automata function is applied that requires a number of neighbouring cells to also contain climax forest. This rule represents the need to have local vegetation for seed dispersal.

2.8

Net primary productivity NPP_{j, t} [gC m^{2 -1} yr^-1] is a function of precipitation and temperature, calculated based on the Miami model (Lieth 1972) as:

(2)

where R_j,t is precipitation [mm] and T_j is temperature [degrees C].

2.9

For each cell, agricultural productivity AG_j,t is calculated as:

(3)

where SP_j is soil productivity (FAO 2007) [index 1-100], S_j is slope [%], WF_j is water flow calculated as the sum volume of water agents traversing any given cell j, as depicted in Figure 2, and SD_j,t is soil degradation [% loss of productivity]. Soil degradation occurs at a constant rate of 1.5% per time step for each cropped cell. This constant rate of soil degradation is obviously a simplifying assumption, as soil management and maintenance of soil productivity is in itself a complex social-ecological system. The β parameters are weightings for calibration, and presented in Heckbert (2012).

2.10

Ecosystem services are modelled by quantifying the availability provisioning services (as defined in MA 2005; TEEB 2010) relating to water, food, and raw materials. This is a subset of ecosystem services and does not include a full set of indicators which would incorporate supporting services (for example erosion prevention), habitat services (such as maintenance of genetic diversity) or cultural services (such as inspiration for culture, art, and design). The current ecosystem services equation incorporates a subset of four important services provision based on arable soils, precipitation, access to available freshwater, and timber resources. Ecosystem servicesES_j are calculated as:

(4)

where AG_j,t is taken from equation 3, R_j,tis taken from equation 1, WF_j,tis the simulated water flow volume, and F_j,tis the forest state [1-3], ESD_j,tis a catch-all proxy variable for all other ecosystem services degradation [%] as a function of population density.

Anthropogenic functions: agriculture, trade, and demographics

2.11

Each settlement agent i maintains at least one cell j for generating agricultural yield. Settlements perform an agriculture benefit-cost assessment considering the costs of production, travel cost given the distance of the cell from the settlement site, and with larger settlements achieving economies of scale, modelled as;

(5)

where BCA_j,t is the total benefit provided from agriculture yield, κ_j, α, and φ, are crop yield and slope parameters, AG_j,t is again taken from equation 3, γ is the establishment cost of agriculture (annual variable costs), O_j is the agriculture travel cost as a function for distance from the city and a per km cost parameter, and P_j,t is population of the settlement.

2.12

The benefit-cost of agriculture function generates yields that are spatially distributed based on individual conditions of the cells and the location of settlements. Costs of production, including distance from settlements, results in adding cropped cells, generating yield and increasing population, which in turn add more cropped cells, but causes soil degradation. The system adjusts over time in response to the spatially-explicit agricultural benefit-cost.

2.13

A series of functions represent trade within a spatially connected network of agents. It is assumed that through the process of specialization, settlements that are connected to one another within a network will generate benefits from trade. It is assumed a larger network produces greater trade benefits, and also the more central a settlement is within the network, the greater the trade benefits for that individual settlement. To model these benefits, settlements are connected via a network of links that represent trade routes. As a simplifying assumption of how they connect together, it is assumed when a settlement reaches (or drops below) a certain size, they will add routes (or allow routes to degrade) to nearby settlements within a 40 km radius. At each time step, the size of the local network is calculated as well as each settlement's centrality within that local network, further discussed below.

2.14

Combining the functions for agriculture, ecosystem services, and trade benefit, total real income per capita RI_i is calculated as:

(6)

where N_j,t is the network size [# nodes], C_j,t is the centrality [degree] and TC_i is the travel cost, and ϑ parameters are prices for agriculture, ecosystem services, and trade, respectively. Benefits from agriculture are calculated only for cells under cropping production AGJ_i,t = 1 … n whereas ecosystem services are calculated encompassing the entire 'area of influence' of each settlement IJ_i,t = 1 … m which is based on the population size of the settlement, increasing linearly to a maximum of 40 km in diameter for the most populous settlements (those with populations greater than 15000 people), as taken from Heckbert et al. (in press) and interpreted from Chase and Chase (1998). Travel cost measures the relative 'friction' of different land cover types, and is represented as:

(7)

where S_j is slope and WF_j,t is simulated water flow volume, both described in previous equations, resulting in areas of higher slope being relatively more costly to travel through, mitigated by the presence of flowing water for canoe transport.

2.15

The terms NJ_i,t network size, and C_i,t centrality in equation 6 are calculated using a network search algorithm custom-built for this application. The network search algorithm is the most detailed procedure in the model, and uses a recursive search function using 'network walking agents' who travel every possible combination of routes along a network path. These report to the network nodes (settlements) the total size of the local cluster, and the position of the node within that cluster (degree). The network search algorithm is the lengthiest procedure in terms of amount of code written, is the most computationally expensive procedure in the model (more than doubles the total run time of a simulation of 650 time steps). The outcome of the procedure is depicted in Figure 3 for time step 250, with the network links coloured black and the range of values on the visual scale halved compared to following figures of the same indicator, for visual purposes.

Figure 3. Simulated 'trade strength' [last term of equation 6] (red on white background, darker red being greater trade strength) based on centrality, size of local network cluster, and travel cost of each settlement (temple icon), connected via trade routes (black links). A local cluster is identified for visual purposes with a blue oval over the same cluster, with a) and b) being zoomed-in views of the wider system depicted in c).

2.16

We see in Figure 3 an area highlighted by a blue oval, zoomed-in at three different scales. In Figure 3a), we see that this local cluster is not connected to adjacent clusters, and N_i,t network size [# of nodes in network] is small. In Figure 3b) we see that the trade strength (red colouring) is greater for larger-sized local clusters given the larger value of N_i,t. Lastly, in Figure 3c) we observe that central nodes in larger clusters (those located in the darkest red regions), have greater trade strength compared to nodes at the periphery of the same cluster.

2.17

The value of each contribution to RI_i,t is determined by the weighting ϑ, which is static for ϑ_AG and ϑ_ES; however, the value of trade ϑ_TRt is dynamic. Specifically, ϑ_TRt increases each time rainfall decreases, according to the climate variation assumptions. The assumption here is that settlements specialize production within an overall trade network to increase the value of trade goods relative to other commodities. The effect is to linearly increase the value of trade each time the climate cycle is in decline. This assumption is obviously rudimentary, and further formulations are explored.

2.18

After determining RI_i,t, settlement demographics account for births, deaths, and migration. The birth rate is assumed to remain constant at 15%, while death rate and out-migration decrease linearly with increased RI_i,t per capita, with a maximum out-migration rate of 15% and a maximum death rate of 25% per annum. Settlements with a population below a minimum number required to maintain subsistence agriculture are deleted. Settlements that register out-migration above a minimum threshold of the number of people required to maintain subsistence agriculture create a 'migrant agent'. The migrant agent uses a utility function (Heckbert et al. 2010) to select locations to create a new settlement. The migration utility function is calculated as:

(8)

where λ parameters are weightings for travel cost and ecosystem services, and ES_j,t is taken from equation 4, and D_j is the distance from the origin settlement to the potential new settlement site.

2.19

Model outcomes are dependent on parameterization of equations 1-8. This stage of model development is offered for proof of concept, for transparency, and to receive feedback on model assumptions. Model code and parameterization is available for download via the http://www.openabm.org website in Heckbert (2012).

Results

3.1

This section describes results and sensitivity analysis from simulation runs. This proof-of-concept stage model reveals familiar patterns analogous to the Maya historical timeline and offers a quantitative basis for these assessments. The model was initialised according to default values (Heckbert 2012). Simulation runs were conducted with a temporal extent of 650 time steps with a spatial resolution of 20km². Figures 5-7 describe model results, depicting mean values for 20 runs (with confidence intervals [a=0.05] where relevant), and reporting on social-ecological indicators.

3.2

Figure 4 presents spatial outcomes for five indicators, at 200 time step intervals. Population density, forest condition, settlement trade strength, soil degradation, the condition of ecosystem services and forests each contain a narrative describing the development and reorganisation of the simulated social-ecological system. By time step 200, settlements have expanded into all regions, first occupying areas with greater ecosystem services and progressively growing with agriculture development. Population densities are higher in areas where settlements have clustered and formed local trade connections. By time step 400, as the value of trade increases, the population dramatically increases, extending local trade connections to 'global' connectivity. A notable fringe exists between the connected network and the periphery . The centre of the global trade network is approximately in the broad region where ancient Maya captials of Tikal and Caracol existed. The condition of the forest is markedly changed, with only small patches of climax forest remaining in agriculturally unsuitable areas, forming ecological refugia within the near-completely settled landscape. Dramatically, by time step 600 the trade network has disintegrated, the centre of the most densely populated areas is nearly entirely abandoned, leaving only a small number of locally connected settlements of any notable size in what was once the fringe of the globally connected network. Abandoned cropland and significantly decreased fuelwood harvesting allows broad-level secondary regrowth, and climax forest eventualy expands out from its refugia to an extent similar to pre population expantion levels. However, soils and ecosystem services remain severely impacted and large scale resettlement of degraded areas is not possible.

Figure 4. Population density, forest condition, settlement trade strength, soil degradation and the condition of ecosystem services are shown at 200 time step intervals. Darker colouring shows increased a) population density (blue), c) trade strength (red), d) soil degradation (red), and e) ecosystem services (green). Forest condition b) depicts three states of cleared / cropped cells (yellow), secondary regrowth (light green) and climax forest (dark green).

3.3

The model reports quantitiative indicators through time and can be used to drill down in order to explore the dynamcis of the development and reoganisation observed. Figure 5 presents the total population of all simulated settlements and contributions to real income by ecosystem services, agriculture and trade, respectively. In the first _ of the simulation, ecosystems services provide the majority of value. Ecosystem services values are superseded by agriculture by time step 150, and both are superceded by trade around time step 350. Neither recover as trade values decrease in the latter half of the simulation run, and population adjusts accordingly.

Figure 5. Total population [# people, primary axis] of all simulated settlements over time, and contributions to real income by ecosystem services, agriculture and trade [# proxy value units, secondary axis]. Ecosystem services values are superceded by agriculture by time step 150, and both are superceded by trade around time step 350. Neither recover as trade values decrease in the latter half of the simulation run.

3.4

The rapid change in the value of trade can be explained by examining Figure 6, which depitcs the total number of settlement nodes, the number of nodes within the largest cluster, and total natural capital, which is represented as the total sum of ecosystem services taken from equation 4. The network grows from local clusters to a near-globally connected system. Periodic perturbations give the clusters structure which inevitably form the 'skeleton' of the global structure. However, when natural capital has reached its lowest level, perturbations result in cascading failure in the network. Although natural capital recovers, the number of settlements does not due to the lack of trade network structure.

Figure 6. Total number of settlement nodes and number of nodes within the largest cluster [primary axis], and total natural capital [total sum of ecosystem services values, secondary axis]. The network grows from local clusters to a near-globally connected system through growth in link connections and periodic perturbations which give the clusters structure. However, when natural capital has reached its lowest level, perturbations result in cascading failure in the network.

3.5

Figure 7 depicts soil degradation and forest condition by three states of cleared / cropped cells (state 1), secondary regrowth (state 2) and climax forest (state 3). The initial period, corresponding to roughly the first third of the simulation run, is descibed by accelerated decline of climax forest and inhibited regrowth is a result of cropping and timber harvesting for construction and fuelwood as a result of increasing population levels. The following period, corresponding to roughly the second third of the simulation run, shows continued repressed regrowth, and increased areas of cleared/cropped land as population pressure results in marginal lands being put under agricultural production. As a result, soil degradation increases at its highest rate during this period. The last third of the simulation run shows a rapid decline in cleared/cropped land with corresponding large scale secondary regrowth, and eventual succession into climax forest which recovers to near pre-development levels.

Figure 7. Forest condition [% of total area, primary axis] by three states of cleared / cropped cells (state 1), secondary regrowth (state 2) and climax forest (state 3), and soil degradation [total proxy soil productivity units lost, secondary axis]. Accelerated decline of climax forest and inhibited regrowth is a result of cropping and timber harvesting. Increased cleared/cropped area results in higher soil degradation rates, with soils beginning to recover after climax forest again becomes dominant.

3.6

The timeline for each indicator presented in Figures 5-7 reveals the underlying dynamics of the modelled social-ecological system. Table 1 presents an integrated narrative of model indicators at 100 time step intervals, akin to the analogy of examining a Rubik's Cube by looking at the patterns on all sides simultaneously (van der Leeuw 2012). This complex systems perspective highlights that a reductionist method of identifying factors contributing to resilience and vulnerability cannot be used to identify a single cause (or a linear progression of single causes) of the unravelling of the social-ecological system.

Table 1: Narrative description of modelled indicators for population, forests, soils, and trade in 100 time step intervals
	Time Step
	100	200	300	400	500	600
Population	Number of settlements grows based on migration to unsettled areas with high ecosystem services	Population increases with agriculture development	Population increases with trade value	Population reaches height, significant perturbations with climate variability	Precipitous decline in population	Population does not recover to prior levels and settlements of any significant size are located at the former fringe of the most densely populated areas
Forests	Gradual decline of climax-dominated forest and gradual increase of secondary regrowth and cleared / cropped areas.	Roughly even distribution of three forest states.	Climax forest reduced to 10% of landscape area, timber harvesting and cropping suppresses secondary regrowth.	Cleared / cropped land accounts for 70% of land surface, climax forest reduced to refugia.	Broad-scale secondary regrowth.	Climax forest begins to emerge from refugia to nearly regain pre-development levels.
Soils and Agriculture	Agricultural expansion into prime soils.	Increase of cropped areas in broadly settled landscape.	Expansion of agriculture into marginal and more distant locations, soil degradation rate increases.	Broad-scale cropping albeit with diminishing returns relative to land area under production. Soil degradation reaches height.	Abandonment of crops and loss of economies of scale in agriculture production, legacy soil degradation inhibits re-sowing crops.	Soils remain unproductive and areas suitable for agriculture are predominantly in the former fringe of densely cropped areas.
Trade and network structure	Little connectivity of trade routes, with small network cluster sizes of <10 settlements.	Local connectivity of trade routes into well-defined clusters of <20 settlements.	Regional-scale connectivity of previous local clusters, with key settlements becoming critical network links.	Near-global network emerges and expands to fringes. Trade value supersedes contributions from agriculture or ecosystem services.	Global network fully connected through key locations, overall network structure vulnerable to perturbance.	Cascading network failure results in existence of only small local networks at the former fringe.

Discussion and conclusions

4.1

The MayaSim model presents a proof-of-concept agent-based model representing key elements of the ancient Maya social-ecological system. Agents, cells, and networks are programmed to represent demographics, trade, agriculture, soil degradation, ecosystem services, climate variability, hydrology, primary productivity, and forest succession. Simulating each of these in combination allows patterns to emerge at the landscape level, effectively growing the social-ecological system from the bottom up. This allows investigation of what conditions lead to the development of densely populated and highly interconnected human geography of the ancient Maya, and reveals dynamics of how past societies impacted their environment and vice versa.

4.2

The model verification process involved debugging and individually testing modules within the overall model. Through this process, it was observed that model outcomes were particularly sensitive to the relationship between soil degradation and the rate of increase in trade value, again within the range of non-extreme climate variation. The relationship lies in the fact that settlements with higher per capita real income, due to trade flows, increase in population and local marginal lands are put under production. Areas with better soils are able to maintain populations and trade connections over multiple cycles of disturbance and in the end become critical nodes in the trade network. This in turn increases local soil degradation. When soil degradation becomes severe, a critical node in the overall global trade network may not be able sustain sufficient local population to maintain trade connections. The demise of a globally significant node can then result in cascading failure in the network, with effects not necessarily located nearby, but everywhere on the network.

4.3

The verification process and manual debugging over thousands of runs produced an informal sensitivity analysis, which revealed that that the modelled Maya social-ecological system, under given assumptions, does not always reach a large population 'peak'. It is only under a narrow set of conditions that the population is able to grow large enough, and with the proper spatial distribution, to result in the globally connected trade network which in turn allows for further population growth. In order to get the peak which is analogous to the observed prehistorical record, a relative balance in the set of parameters relating to trade, ecosystem services, and agricultural real income must occur. In other words, the equations relating to demographics, soil degradation, and trade value must all be within a specific range for the population to rise substantially enough to create patterns that are recognisable as Preclassic, Classic and Post Classic temporal periods. The two primary parameters that seem to require balance are the increase in trade value and the rate of soil degradation.

4.4

Further research will focus on calibrating model inputs from agronomic models, inputting finer resolution soil and vegetation data, and running the model over a generated paleoclimate reconstruction. Rudimentary assumtions on how trade value changes can be replaced with other candiate assumtions. Future validation is possible using GIS-based archeological site data for settlement locations over time, comparison of simulated built areas to LIDAR data, comparing pollen core content to simulated catchment land cover, as well as network connectivity of settlements as evidenced through ceramic interactions over time.

4.5

Revisiting research questions in light of the model results, we can identify critical dynamics of how the Maya came to develop into a highly interconnected and densely populated social-ecological system. These processes can be interpreted into a temporal pattern of development and reorganisation, that in the simulated model produces outcomes roughly analogous to the archaeological description of the Maya Preclassic (< model time step 200), Classic (roughly time step 200-450), and Postclassic periods (> model time step 450). From a complex systems science perspective, there is no one cause of the dramatic re-organisation which occurs to mark the end of the Classic period, and notably in this model there is no 'drought', merely cyclical climate variability whose effect completely depends on the entirety of the rest of the state variables in the social-ecological system. Climate variability, soil degradation, deforestation, demographic pressure, and the physical configuration of trade networks are all factors which concurrently contribute to resilience or vulnerability, and cannot be 'ordered' or 'prioritised' into any given sequence of events which causes the re-organisation into the Postclassic period. Complex systems modelling can however identify how these variables interact, and provides a picture of how an interconnected system responds given its embedded fast-moving variables such as land cover change and trade connections, meso-speed variables such as demographics and climate variability, as well as slow-moving variables such as soil degradation.

Acknowledgements

The author would like to acknowledge the contributions of Christian Isendahl, Joel Gunn, Andrew Reeson, Simon Brewer, Tim Baynes, Vernon Scarborough, Arlen Chase, Diane Chase, Robert Costanza, John Murphy, Derek Robinson, Nicholas Dunning, Carsten Lemmen, Lael Parrot, Timothy Beach, Sheryl Luzzadder-Beach, David Lentz, Paul Sinclair, Carole Crumley and Sander van der Leeuw. This project was supported by Alberta Innovates Technology Futures, Portland State University, Arizona State University, Uppsala University, and University of Cincinnati.

References

CHASE, A.F. and Chase, D.Z. (1998). Scale and Intensity in Classic Period Maya Agriculture: Terracing and Settlement at the 'Garden City' of Caracol, Belize, Culture and Agriculture 20(2):60-77. [doi:10.1525/cag.1998.20.2-3.60]

FAO (2007). Soil Production Index. http://www.fao.org:80/geonetwork?uuid=f7a2b3c0-bdbf-11db-a0f6-000d939bc5d8

FAO / IIASA / ISRIC / ISSCAS / JRC, (2009). Harmonized World Soil Database (version 1.1). FAO, Rome, Italy and IIASA, Laxenburg, Austria.

FARR, T.G., Rosen, P.A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D. (2007), The Shuttle Radar Topography Mission, Rev. Geophys., 45, RG2004. [doi:10.1029/2005RG000183]

HECKBERT, S., Adamowicz, W., Boxall, P., & Hanneman, D. (2010). Cumulative Effects and Emergent Properties of Multiple-Use Natural Resources. Lecture Notes in Computer Science, 5683, 1-13. [doi:10.1007/978-3-642-13553-8_1]

HECKBERT, S. (2012). MayaSim: An agent-based model of the ancient Maya social-ecological system. http://www.openabm.org/model/3063/version/3

HECKBERT, S., Isendahl, C., Gunn, J., Brewer, S., Scarborough, V., Chase, A.F., Chase, D.Z., Costanza, R., Dunning, N., Beach, T., Luzzadder-Beach, S., Lentz, D., Sinclair, P. (2014, in press). Growing the ancient Maya social-ecological system from the bottom up. In: Isendahl, C., and Stump, D. (eds.), Applied Archaeology, Historical Ecology and the Useable Past. Oxford University Press.

HIJMANS, R.J., S.E. Cameron, J.L. Parra, P.G. Jones and A. Jarvis, (2005). Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965-1978. Available at: http://www.worldclim.org/ [doi:10.1002/joc.1276]

LIETH, H., 1975. Modeling the primary productivity of the world. In: Lieth, H., Whittaker, R.H. (Eds.), Primary Productivity of the Biosphere. Springer-Verlag, New York, pp. 237-263. [doi:10.1007/978-3-642-80913-2_12]

MA (2005). Millennium Ecosystem Assessment. Ecosystems and Human Well-being: Synthesis. Island Press. Washington, DC.

PRUFER, K., Scarborough, V.L. Chase, A.F., Chase, DZ., Cobos, R., Dunning, N., Gunn, J., Fedick, S., Fialko, V., Iannone, G., Lentz, D., Liendo, R., Lucero, L., Sabloff, J., Tainter, J., Valdez, F., and van der Leeuw, S. (2011). IHOPE Maya: Resilience and Rigidity in the Development and Disintegration of Complex Societies in the Tropical Lowlands of Mesoamerica, paper presented at Resilience 2011, Arizona State University, Tempe, Arizona.

REANEY, S.M. (2008). The use of agent based modelling techniques in hydrology: determining the spatial and temporal origin of channel flow in semi-arid catchments. Earth Surface Processes and Landforms. 33, pp. 317-327. [doi:10.1002/esp.1540]

TEEB (2010). The Economics of Ecosystems and Biodiversity. Mainstreaming the Economics of Nature: A synthesis of the approach, conclusions and recommendations of TEEB. Progress Press, Malta.

VAN DER LEEUW, S., Costanza, R., Aulenbach, S., Brewer, S., Burek, M., Cornell, S., Crumley, C., Dearing, J., Downy, C., Graumlich, L, Hegmon, M., Heckbert, S., Hibbard, K., Jackson, ST., Kubiszewski, I., Sinclair, P., Sörlin, S, and Steffen, W. (2011). Toward an Integrated History to Guide the Future. Ecology and Society. (16)4. [doi:10.5751/es-04341-160402]

WILENSKY, U. (1999). NetLogo. Center for Connected Learning and Computer-Based Modeling, Northwestern University. Evanston, IL. http://ccl.northwestern.edu/netlogo/.