• MayaSim: An Agent-Based Model of the Ancient Maya Social-Ecological System

  • 마야심: 고대 마야의 사회적-생태학적 시스템의 행위자 기반 모형

Image of 1996 article Image of 1996 article
  • MayaSim: An Agent-Based Model of the Ancient Maya Social-Ecological System
  • Journal of Artificial Societies and Social Simulation 16 (4) 11
  • Received: 20-Aug-2012    Accepted: 02-Jun-2013    Published: 31-Oct-2013
  • Article Tools
  • 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

  • 서론

  • 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.
  • 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.
  • 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.
  • 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).
  • 이 모델에는 사회적-생태학적 시스템의 복잡한 역학관계를 좀 더 이해하려는 목적이 있으며, 시스템의 지속가능성을 예측 변수로서 탄력성의 정량적인 지표들을 시험하려는 의도가 있다. 통합적인 행위자 기반, 세포 자동자와 네트워크 모델은 넷로고 소프트웨어(Wilensky 1999)로 작성되었다. 모델, 코드, 문서 전체는 www.openabm.org 웹사이트에서 Heckbert (2012)로 공개되어 있다. 나아가 마야 고고학 문헌 맥락에서 모델의 상세 설명은 Heckbert 외 (인쇄중)에서 확인할 수 있다.
  • Methods

  • 연구 방법

  • 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 km2 with a 20 km2 resolution. Temporal extent is 650 times steps, each representing roughly 2 years.
  • 마야심 모델은 행위자가 격자 지형에 위치하는 것으로 정착을 표현한다. 이 모델은 넷로고 소프트웨어(Wilensky 1999)를 사용해서 작성되었다. 이 소프트웨어 인터페이스(그림1)는 모델과 상호작용하기 위한 사용자 인터페이스와 모델 데이터를 추적하는 그래프를 통해 모델의 공간적 뷰를 표현한다. 이 뷰는 모델 나의 서로 다른 공간 데이터 레이어를 보기 위해 변경될 수 있다. 모델은 20km2 해상도로 516,484km2 정도의 공간에서 동작한다. 시간 정도는 650 회에 걸쳐 시뮬레이션되는데, 각 횟수는 대략 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 20km2 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.
  • 모델을 초기화할 때, 넷로고 GIS 확장기능을 사용해서 기본 GIS (지리정보시스템) 레이어가 들어간다. 정적 변수들은 셀에 내재되고, 동적 변수들은 사용자 인터페이스의 기본 값으로 재설정된다. 정착 행위자들은 공간에 무작위로 초기화된다. 삽입된 공간 데이터에는 고도, 경사(Farr 외 2007), 토양 생산력(FAO 2007), 온도, 강수량(Hijmans 외 2005)가 포함된다. 토양 생산력에 대한 데이터는 넷로고의 디퓨즈(확산) 함수를 사용해서 조금 수정되었다. 폴리곤 경계를 따라 생기는 차이를 완화하기 위함이었다. 그러나 다른 모든 데이터는 인용한 원본 그대로이다. 넷로고 GIS 확장기능을 사용해서(Wilensky 2009 를 보라) 데이터를 20km2 해상도로 다시 샘플링했다. 토양 생산력 데이터를 제외한 모든 GIS 데이터는 고해상도이고, 고화질로 간주해도 좋다. 토양 생산력 데이터는 좀 더 나은 해상도의 지역 데이터 혹은 특정 국가의 데이터로 교체할 수 있었지만, 연구 지역은 현재 여러 국가로 나눠져있어서 그러기에는 데이터에 일관성이 사라졌다. 그래서 가능한 전세계의 일관성있는 데이터를 사용하는 게 최선이었다. 좀 더 나은 공간적 해상도와 국경에 상관없이 일관적인 데이터를 위해 향후 연구에서는 '세계 통일 토양 데이터베이스'(FAO 외 2009)의 데이터로 교체하려고 한다.
Image of 1996 article
  • 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.
  • 그림1. 마야심 모델은 쌍방향 제어, 공간 뷰, 모델 데이터를 추적하는 그래프를 인터페이스로 가진다. 행위자들은 셀 지형에 있으며 네트워크에서 링크로 연결된다.
  • 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

  • 시뮬레이션은 생물물리학적 변수로 유수량과 순 일차 생산력을 계산하는 것으로 시작한다. 나아가 이들은 산림 천이, 농업 생산량, 생태계 서비스를 계산하는데 사용한다. 정착 행위자들은 작물로 농업 이익을 만들기 위해 공간 지형과 상호작용한다. 여기에서 지역 생태계 서비스로부터 이익을 얻고, 지역 교역 네트워크에서 교역 이익을 만들어낸다. 농업, 생태계 서비스, 교역의 복합 이익은 이주를 포함한 인구성장의 원인이 된다. 통합 시스템을 시뮬레이션하는 것으로 사회적-생태학적 시스템이 시간에 따라 기능하는 방법을 밝혀냈다. 추가의 행위자들은 새로운 장소에 정착한 '이주' 행위자, 수문학(水文學)상 표면 흐름을 따르는 '빗방울' 행위자, 서로 연결된 정착 노드를 가로지르는 '네트워크 검색' 행위자를 포함한다. 이들은 네트워크 통계 계산에 쓰인다. 이어지는 내용으로 모델에서 추척한 생물물리학적, 인류발생적 과정을 설명할 것이다. 생물물리학적 기능들: 기후, 수문학, 토양, 일차 생산력, 숲, 생태계 서비스 등을 살펴볼 것이다.

  • 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:
Image of 1996 article
  • (1)
  • where Rj,T is precipitation [mm] for cell j at initial time step T, and CLn 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. DFj is the distance [km] of each cell from the top northwest corner of the map and is the furthest distanced cell from this point. RCt 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.
  • 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 Rj,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.
Image of 1996 article
  • 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.
  • 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.
  • 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 * NPPmax j, t / NPPj, t years for secondary regrowth and 100 * NPPmax j, t / NPPj, 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.
  • Net primary productivity NPPj, t [gC m2 -1 yr-1] is a function of precipitation and temperature, calculated based on the Miami model (Lieth 1972) as:
Image of 1996 article
  • (2)
  • where Rj,t is precipitation [mm] and Tj is temperature [degrees C].
  • For each cell, agricultural productivity AGj,t is calculated as:
Image of 1996 article
  • (3)
  • where SPj is soil productivity (FAO 2007) [index 1-100], Sj is slope [%], WFj is water flow calculated as the sum volume of water agents traversing any given cell j, as depicted in Figure 2, and SDj,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).
  • 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 servicesESj are calculated as:
Image of 1996 article
  • (4)
  • where AGj,t is taken from equation 3, Rj,tis taken from equation 1, WFj,tis the simulated water flow volume, and Fj,tis the forest state [1-3], ESDj,tis a catch-all proxy variable for all other ecosystem services degradation [%] as a function of population density.
  • Anthropogenic functions: agriculture, trade, and demographics

  • 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;
Image of 1996 article
  • (5)
  • where BCAj,t is the total benefit provided from agriculture yield, κj, α, and φ, are crop yield and slope parameters, AGj,t is again taken from equation 3, γ is the establishment cost of agriculture (annual variable costs), Oj is the agriculture travel cost as a function for distance from the city and a per km cost parameter, and Pj,t is population of the settlement.
  • 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.
  • 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.
  • Combining the functions for agriculture, ecosystem services, and trade benefit, total real income per capita RIi is calculated as:
Image of 1996 article
  • (6)
  • where Nj,t is the network size [# nodes], Cj,t is the centrality [degree] and TCi 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 AGJi,t = 1 … n whereas ecosystem services are calculated encompassing the entire 'area of influence' of each settlement IJi,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:
Image of 1996 article
  • (7)
  • where Sj is slope and WFj,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.
  • The terms NJi,t network size, and Ci,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.
Image of 1996 article
  • 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).
  • 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 Ni,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 Ni,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.
  • The value of each contribution to RIi,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.
  • After determining RIi,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 RIi,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:
Image of 1996 article
  • (8)
  • where λ parameters are weightings for travel cost and ecosystem services, and ESj,t is taken from equation 4, and Dj is the distance from the origin settlement to the potential new settlement site.
  • 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

  • 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 20km2. 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.
  • 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.
Image of 1996 article
  • 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).
  • 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.
Image of 1996 article
  • 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.
  • 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.
Image of 1996 article
  • 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.
  • 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.
Image of 1996 article
  • 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.
  • 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
  • 100200300400500600
  • PopulationNumber of settlements grows based on migration to unsettled areas with high ecosystem servicesPopulation increases with agriculture developmentPopulation increases with trade valuePopulation reaches height, significant perturbations with climate variabilityPrecipitous decline in populationPopulation does not recover to prior levels and settlements of any significant size are located at the former fringe of the most densely populated areas
  • ForestsGradual 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 AgricultureAgricultural 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 structureLittle connectivity of trade routes, with small network cluster sizes of