Research ArticleANTHROPOLOGY

The demise of Angkor: Systemic vulnerability of urban infrastructure to climatic variations

See allHide authors and affiliations

Science Advances  17 Oct 2018:
Vol. 4, no. 10, eaau4029
DOI: 10.1126/sciadv.aau4029


Complex infrastructural networks provide critical services to cities but can be vulnerable to external stresses, including climatic variability. This vulnerability has also challenged past urban settlements, but its role in cases of historic urban demise has not been precisely documented. We transform archeological data from the medieval Cambodian city of Angkor into a numerical model that allows us to quantify topological damage to critical urban infrastructure resulting from climatic variability. Our model reveals unstable behavior in which extensive and cascading damage to infrastructure occurs in response to flooding within Angkor’s urban water management system. The likelihood and extent of the cascading failure abruptly grow with the magnitude of flooding relative to normal flows in the system. Our results support the hypothesis that systemic infrastructural vulnerability, coupled with abrupt climatic variation, contributed to the demise of the city. The factors behind Angkor’s demise are analogous to challenges faced by modern urban communities struggling with complex critical infrastructure.


The management of modern infrastructural networks that may be vulnerable to rare but catastrophic failure is of great concern (19). The destabilizing effect of extreme weather events, associated with planetary-scale climate change (10), creates new risks to the stability of the urban environments in which most humans now live. These risks have become more acute as urban conglomerations become larger, more complex, and trend toward low-density, dispersed settlement forms (11). Low-density urbanism is known to be stable over long time periods (12) but demonstrates systemic instability above a given scale and density (13). The occurrence of very large low-density urban settlements in the preindustrial period provides an opportunity to investigate the stability of urban networks over long time periods, particularly when stressed by climatic variation.

Understanding of and adaptation to high-impact, low-frequency weather events are hampered by limited long-term data. This is especially so with respect to the complex social response to infrequent climatic perturbation (14). To our knowledge, quantitative evidence of perturbation in networked urban infrastructure has not been demonstrated in the archaeological or historical record (15). In cases where infrastructure can be directly analyzed or credibly reconstructed, and where other aspects of the material culture or historical archive are of adequate fidelity, it should be possible to resolve the rate and extent of failures or outages in infrastructural networks and to identify the specific articulation between forcing processes, failing infrastructure, and social response. Here, we analyze the largest urban infrastructural network of the preindustrial world, at Angkor in modern Cambodia, Southeast Asia. We assess the topology of the network and its vulnerability to catastrophic failure in response to climatic perturbation, and consider the implications for the resilience of networked urban infrastructure.

System theory, complexity, and cascading failure

In a complex system with many interacting components, a catastrophic response may result from even a small local change or from a series of accumulated changes. Such a response may be amplified by positive feedbacks, leading to a tipping point (16): a threshold beyond which further change forces the system to rapidly, and often irreversibly, adopt a new structure and/or functions. For example, a tipping point in a forest ecosystem may be reached when an initial perturbation brought about by deforestation triggers a positive feedback between reduced regional precipitation, increased risk of wildfire, and accelerated forest dieback, leading to further reductions in precipitation, and so on (17, 18). A physics-based perspective on a tipping point encompasses a critical regime characterized by phase transition when some macroscopic features of the system (e.g., density) undergo a sudden transformation in response to a change in an underlying parameter (e.g., temperature) affecting the system composition and/or the strength of internal interactions. This criticality may be an emergent property of self-organizing systems (19).

One key symptom of systemic vulnerability in complex networks is the tendency for relatively small perturbations to cascade outward from a primary failure to trigger further failures in other dependent parts of the network that may be distal to the primary failure (8). This produces an important nonlinearity, where the consequences of perturbation may be much larger and more significant than the perturbation itself. The study of cascading failures in complex networks has been a strength of network theory (20), providing realistic models of cascading failures in electrical grids (21), traffic jams in transportation networks (22), as well as various models of the spread of epidemics (23). In each of these domains, cascade modeling could be used to identify the presence of phase transitions and values of critical parameters. In addition, the relationship between network topology and resilience can be explicitly examined. For example, scale-free networks (24, 25) may be resilient to small, frequent perturbations that may affect one or more nodes (26) but are extremely vulnerable to catastrophic failure should one or more hubs fail (27). In instances where two or more networks share components, as is frequently the case for modern infrastructure (for example, railway lines and power grids), cascading failures can propagate rapidly between coupled networks (28). If complex infrastructural networks in premodern urban environments share common topological and functional characteristics with modern complex networks, then it is possible that catastrophic infrastructural failure may also have occurred in the past. In some cases, systemic vulnerabilities may have played a role in the deterioration of social or economic conditions within past urban environments, or even triggered tipping points, leading to large-scale urban abandonment.

Angkor’s water distribution network

Angkor was the world’s most extensive city by the 13th century CE, covering an area of approximately 1000 km2 (Fig. 1) (29). Angkor developed from two urban centers around Ak Yum (7th and 8th centuries CE) and the Bakong (8th and 9th centuries CE). Epigraphic sources date the beginning of the Angkor period to 802 CE. Over the subsequent 600 years, the city established a number of infrastructure networks, the most extensive of which was the network of canals, embankments, moats, and reservoirs used to capture, store, and distribute surface water resources (2933). This infrastructural network served several functions, most critically flood regulation and irrigation for agriculture, and developed to become both spatially extensive [covering up to 1200 km2; (31)] and internally complex, with many thousands of individual components that were, to a greater or lesser degree, dependent on other components elsewhere in the system. By the end of the 11th century, this network had effectively reached its maximum lateral extent (30). The emplacement in the network of the last large reservoir (the Jayatataka) in the late 12th century CE represents a process of urban reorganization and infrastructural infill that mirrors the shift toward dense, highly structured city blocks within and surrounding the temple enclosures (33). The network was likely continually maintained and modified during its development, with some specific instances of modification known in the archeological remains of the network (31). Maintenance of the network appears to have ceased in concert with its breakdown and the coincident demise of the city. Angkor’s water management infrastructure was a true palimpsest, representing hundreds of years of extensification, intensification, and modification, which had developed to become a massive, complex, and interdependent infrastructural network.

Fig. 1 Location of the study site.

(A) Location of the medieval city of Angkor in modern Cambodia. Angkor was the primate city to a kingdom that incorporated most of mainland Southeast Asia at its peak in the 11th and 12th centuries CE. (B) Archeological map of Greater Angkor, after (29, 33, 46, 47). Limit of the study area indicated on (B) represents the watershed boundary.

Angkor foundered during the 15th century. The Cambodian Chronicles report that the king and his administration relocated to “Middle Period” cities in the southeast, close to the modern capital of Phnom Penh (34). While Angkor was never completely abandoned, its population fell precipitously, and many parts of the city were surrendered to the encroaching forest. Attributing this extraordinary social transformation primarily to the sack of Angkor by the Siamese in 1431 CE remains problematic (35), and the emerging consensus is that several factors may have played a role in the decision to abandon the city (33, 36, 37). Among these factors is a coincident period of weakened summer monsoon rainfall and increased intraseasonal variability revealed in tropical tree ring indices (37, 38) and numerous other proxy data sources through the tropical Asia Pacific (3941). Two multidecadal-scale droughts in the mid 14th and early 15th centuries CE bracket a period of unusually intense summer monsoon rainfall that resulted in significant damage to Angkor’s critical water management infrastructure (37). Historian Victor Lieberman (42, 43) argues that the demise of Angkor during this period was part of a larger territorial fragmentation across the whole of the region, in part related to climatic variability. However, even within this context, Lieberman recognizes that Angkor is unusual and that the collapse there was more wrenching and longer lasting than elsewhere on the mainland [(43), p. 27]. No other contemporary state-level polity in mainland Southeast Asia endured social upheaval of the same magnitude or duration.

For many years, it has been suggested that the city’s networked infrastructure played a role in the demise of Angkor (30) and, more recently, that the city’s massive infrastructure articulated in some way with climatic variability to produce that outcome (44). What we have lacked is a clear account of how that articulation operated from a mechanistic point of view and what characteristic(s) of the city’s infrastructure made it particularly vulnerable to catastrophic failure under the impact of internal and/or external perturbation.

Here, we investigate the topological vulnerability of Angkor’s water distribution network to erosion and sedimentation caused by (potentially extreme) variations in rainfall conditions. The model network interrogated here is based on the topology of Angkor’s historic water distribution network as it was in the 13th century CE. Modeling of the system in earlier time periods is not possible because the system is a palimpsest developed over 600 years or more, and earlier network configurations have been lost or obscured by subsequent modification. We can know the topology of the final network configuration, however, as there were no additional modifications to the system after its fragmentation. The model is composed of 1013 edges (primarily excavated canals, but also embankments, earthen dikes, roadways, and so on, that act to control and direct the flow of water) and 617 nodes (canal confluence or bifurcations, reservoirs, temple moats, and so on), which represent points at which edges meet or terminate. These numbers include the super source and super sink (see Materials and Methods). We note that topological and physical vulnerability may be distinct. That is, physical deformation of the structures composing the network does not necessarily change their function in terms of water distribution. Here, the term “topological vulnerability” refers to the network’s function in distributing water from the natural rivers to various use, storage, or discharge points within the network. We apply a simple numerical simulation model of erosion and sedimentation within the water distribution network (Fig. 2) to identify whether these processes could markedly alter the distribution of water through the network.

Fig. 2 Schematic of the numerical model of erosion and sedimentation, illustrated for a single bifurcation in the network.

Channels (edges) are depicted as dark blue lines and junctions (nodes) as pale blue circles. The flow on each edge corresponds to the thickness of the blue lines. Initially, flow is distributed evenly. In response to a heterogeneous increase in flow, one edge begins to erode, increasing its capacity and diverting water away from the neighboring edge. Sedimentation on this depleted edge, in turn, decreases its capacity and diverts more water along the eroded path, resulting in a damaged state in which the initial flow volume is unevenly distributed.


For a given flood magnitude (γ), three “phases” exist in topological damage (Q) as a function of erosion (α) and sedimentation (β) (see Materials and Methods and the Supplementary Materials for details on how these terms are defined and calculated; Fig. 3 and tables S1 and S2). The first phase, which we refer to as phase 1, is characterized by negligible damage to the water distribution network. In this region of parameter space, flooding does not produce erosion. As α increases (and erosion threshold decreases), we come to phase 2, in which erosion occurs. The transition between phases 1 and 2 is characterized by a local maximum in Q, which we attribute to heterogeneity in erosion patterns. Heterogeneity is required for flow distribution to be altered because of erosion (if all channels are eroded in proportion to their initial flows, no change in the overall flow distribution occurs, and Q = 0). If sedimentation were not included in our model, this transition would still exist. This can be seen as an independence of Q on β until after the transition between phases 1 and 2 (Fig. 3B).

Fig. 3 Scaling of topological damage (Q) as a function of α and β (κ = 1, γ = 1) demonstrates three distinct regions of parameter space.

The color plot in (A) represents the entire parameter space in α and β with the transition between regions (i) and (ii) dependent only on the erosion threshold, while the plots in (B) illustrate higher-resolution scans of “alpha” for fixed “beta” values, demonstrating a transition between regions (ii) and (iii) dependent on the thresholds for both erosion and sedimentation. Each datum point presented here corresponds to an average of 100 different flood configurations.

From phase 2, we can move into a feedback regime (phase 3) by increasing β (that is, increasing the propensity for sedimentation). In this phase, incision of the flow path and subsequent flow capture starves adjacent paths of flow, increasing the probability of sedimentation in those edges. This effect propagates to other dependent components of the network.

This unstable feedback regime is strongly dependent on flood magnitude (γ). Comparison of Figs. 3A and 4C demonstrates that increasing γ substantially will expand the region of parameter space occupied by the feedback regime (Fig. 4). For a wide range of α and β values, the damage mechanism robustly produces topological damage as long as the flood magnitude is high enough that the majority of parameter space is within the unstable phase. Very large floods will therefore quickly result in systemic instability in the water distribution network. The establishment of “pluvial” conditions over mainland Southeast Asia in the final decades of the 14th century (37), which is synchronous with radiometrically dated evidence of severe erosion within the network (see Materials and Methods), is likely to have forced the network suddenly into this unstable phase space and likely contributed to the fragmentation of Angkor’s infrastructure network.

Fig. 4 Effects of flood magnitude (γ) on Q.

(A and B) Plots of Q as a function of α for κ = 1, β = 0.001 (A), and β = 0.5 (B). (C) Color plot of Q as a function of α and β for γ = 5 and κ = 1.

While the positions of the transitions between phases 1, 2, and 3, as well as the practical bounds of Q, depend on the topology upon which the simulation is performed, the model is applicable to any directed, acyclic network (see the Supplementary Materials for results from a generic topology).


From a design standpoint, the spatial distribution of damage within this network is informative. While our primary result is valuable for assessing the validity of the hypothesis investigated, the water distribution network itself has unique properties that can only be understood by examining which parts are more or less susceptible to perturbation.

To localize changes to flow distribution within the network, we calculated the average of Q on each edge (Qi)n over n = 3000 repetitions. A map representing these results for a particular set of conditions (γ = 1, α = 0.9, and β = 0.9) is shown in Fig. 5. This map is representative of the relative distribution for other values of control parameters tested as well. That is, while the magnitude of the damage depends on α, β, and γ, its spatial distribution within the network does not. Damage is highest in upstream regions near the canals connected directly to source nodes and dissipates through downstream bifurcations and recombination points.

Fig. 5 A map representing the average damage distribution in the network.

This map was generated for control parameters γ = 1, α = 0.9, and β = 0.9 [region (iii)]. The color of each edge represents the average Q value for that edge 〈Qin, calculated from n = 3000 independent runs.

Some parts of the network appear robust to damage. Even when there is a large amount of water diverted from the main inputs, the downstream branches tend to redistribute the flow in such a way that it equilibrates in subsidiary regions such as the central urban area. On the other hand, the northern and western regions appear to accumulate more damage relative to the central urban area and southeastern regions.

These results indicate that damage is most severe near upstream bifurcations, especially those areas that serve as the primary distribution forks for main inputs. The most robust regions are located downstream of points at which damaged paths recombine. This is because damage can represent either decreased or increased flow. If a depleted and flooded edge recombine at some downstream node, their respective perturbations cancel.

It should be noted that the results presented in Fig. 5 are aggregate. That is, they cannot explain the particular archeological observations made on the real landscape, which were derived from one of an infinite set of possible scenarios. That said, if we are to take these results literally, they indicate that the infrastructural damage to waterways produced by monsoon-driven flooding during the late 14th and early 15th centuries would have had the most marked effects on the northern and western regions of the network, leaving the central and southeastern regions relatively untouched from a functional perspective. The most marked evidence of erosive damage to the water distribution network can be observed to the north and east of Angkor Thom, particularly focused on the line of the Siem Reap River. As the primary offtake from the source, this empirical evidence fits well with the model output. However, excavations in the northwest sector of Greater Angkor remain very limited, and it may be possible that further evidence of erosive damage is apparent in that area also.

Implications for stability of complex flow networks

We use a minimal three-parameter model of the water distribution network of Angkor to determine whether the network was topologically vulnerable to systemic failure when perturbed by low-frequency, high-magnitude weather events. We find that, above a critical flow volume, damage to the network passes a threshold and becomes highly unstable, resulting in the centralization of flow and a cascading functional breakdown of the network. Our model indicates that damage to the network is focused on major distribution nodes close to the flow source (that is, where the network captures flow from the natural river system). This agrees well with observations of the physical damage to the system, and particularly the deep incision on the Siem Reap canal, which appears to have captured the majority of the natural flow from the Puok River. In the model, we find that damage diminishes in downstream reaches because of redistribution of load and the consequent equilibration of damage through distributary reaches.

We conclude that Angkor’s urban infrastructure, particularly its very large and convoluted water distribution network, was prone to cascading failure under the impact of high-amplitude but low-frequency climatic perturbation. This critical behavior was likely to have been the product of centuries of iterative development of water management infrastructure, creating a system that was at once very large, tightly interdependent, and characterized by the reuse of older, redundant infrastructure. It follows that earlier, less extensive, or interdependent phases of the same network may demonstrate different characteristics in terms of topological damage in response to climate forcing, and may have been more resilient than the network as we model it, at the peak of its complexity in the 14th century CE.

We suggest that the functional breakdown of critical urban infrastructure played a role in the broader demise of Greater Angkor as a viable urban complex in the 14th and 15th centuries CE. It seems probable that Angkor’s water distribution network had been further primed for catastrophic failure by long periods of drought during the 13th and 14th centuries, to which the functions of the network, and presumably the agricultural practices of the Khmer, had become equilibrated. Marked changes to the function of important components of the network are apparent in the remains of the system [such as alterations to the Siem Reap River and the reconfiguration of the outlet to the East Baray reservoir at Krol Romeas; (31, 37, 45)]. The abrupt transition, in the middle of the 14th century CE, from prolonged drought to particularly wet years, meant that the network was disequilibrated with prevailing weather. The subsequent adjustment of the system, through erosion of the major offtakes and centralization of flow, exploited the systemic vulnerability shown by our model and resulted in probable breakdown of the network.

There are two primary outcomes of this finding. First, the basic pathology of Angkor’s collapse is analogous to the challenges faced by networked urban infrastructures in the modern world. This was not an exotic catastrophe with no modern analog, nor was it the result of unique exogenous forces. Second, this historic example of the cascading failure of critical infrastructure as a result of climate extremes reemphasizes the importance of building resilience into modern networks, resilience engendered through modularity and redundancy, to cope with high-amplitude/low-frequency climatic extremes, of which we may, at present, have no historical or instrumental precedent.


Map data

There is a long history of archeological remote sensing and mapping at Angkor, and comprehensive archeological maps representing the complete water distribution network have been produced and meticulously ground verified (29, 30, 33, 46, 47). Some components of the network (such as the large southwest-northeast distributor canal from the West Baray and several roads) continued outside of the mapped area, and these were truncated at the watershed boundary (Fig. 1B).

The map was divided into 400, 7840-m2 grid squares in ArcGIS 10, and all features functionally related to the distribution of water in each grid square were given a number, a flow direction, and a dimensionless weight value, according to the coding and classification rules described below (see fig. S1). Metadata for each component (such as its type and geographical location) were also recorded.

Coding and classification of map data

Access to and control of surface and ground water were critical for Angkor’s success as a low-density city, particularly in terms of flood control in the urban core and agricultural productivity in the vast low-density sprawl that surrounds it. Yet, the city’s water distribution network was not constructed with a single integrative function (as a planned, unified single system for water management and distribution) but was a palimpsest created over several hundred years. As such, parts of the system may have ceased to function, have become redundant, or have been repurposed while the city was active [the early centers of Mahendraparvata and Hariharalaya are obvious examples of such; (48, 49)]. The systemic complexity and interconnectedness of the network emerged from its iterative growth, as well as the overprinting, modification, and redundancy of existing infrastructure. While this is widely understood, the complete network is often treated as a single, functioning, engineered system with a set of shared purposes [see, in particular, (32), but also 31)].

That the network is a diachronic artifact means that the functional part of the network at any given point in time may have been much smaller than our model indicates. From a network theory point of view, this means that the network is not topologically static either during its growth or in its operation. Equally, the remains of the network that can be observed (and mapped) in the modern environment represent only the endpoint of the network’s growth, not its evolution. However, irrespective of the intentions of the people who built and maintained the network, or the use of its parts over time, the components of the network are still capable of influencing the flow of surface water because they remain physical topographical features in the landscape. The age or original purpose of an individual component (edge) within the network is then irrelevant if its existence continues to influence water flow. In most cases, individual components can be linked mechanistically to other components in terms of flow connectivity, irrespective of the age of the components or the specific intentions that may have motivated their construction. These simple principles—that extant topographic features (canals, embankments, dikes, reservoirs, moats, and so on) influence the flow of surface water and that captured water flows downslope—form the basis of our approach to coding Angkor’s diachronic water management network.

Here, we took the network as it was represented in the archeological maps as an accurate representation of the water distribution network, as it was (conservatively) well after the end of its growth, in the middle of the 14th century. We were specifically interested in this “time slice” because it coincided with the “pluvial” episode associated with significant disruption to and reorganization of the water distribution network (37). At that time, the Siem Reap canal incised up to 8 m below the land surface, in its middle reaches. This is known to have occurred in the 14th century [no earlier than 1268–1423 2σ calibrated years CE using the IntCal13 database (50), with a median probability of 1345 CE; 1286–1440 2σ calibrated years CE with the southern hemisphere SHCal13 dataset following (51), with an offset of 21 ± 6 years following (52) and (53), with a median probability of 1356 CE]. That event then represents a fundamental reorganization of the water distribution network at the end of the “Angkor 1” drought (37) and represents the final fragmentation and collapse of the network. The frequency of flood events that cause damage to the network must be lower than the rate of response in the receiving system (that is, adjustment of the network topology through erosion or sedimentation) for the model to be applicable. The superdecadal frequency of seasonal rainfall events likely to lead to significant flooding [as reflected by the high Palmer Drought Severity Index values in tropical tree ring indices (37)] indicates that this minimal requirement can be met. Our simulated network incorporates all surface archeological features that can be associated spatially and functionally with the water distribution network and is therefore the best approximation of the network as it was in the middle of the 14th century.


The cascade model used here is a generic, phenomenological model of erosion and sedimentation dynamics with a minimal number of control parameters. We apply this dynamic model to an acyclic, directed tree network representing Angkor’s water distribution network, to test the system’s sensitivity to cascading failure (see fig. S7 for a degree distribution of the network). Parameter validation is the most difficult aspect of modeling complex phenomena such as cascading failures in infrastructural networks, and this is particularly so for ancient networks. For this reason, we chose to design a generic, phenomenological model with well-constrained control parameters that integrate innumerable physical processes (such as rainfall, soil characteristics, human intervention capacity, local vegetation, channel design, etc.), which cannot be directly measured. In our network, initially, each edge i in the graph is characterized only by a weight wi, which is based on the estimated capacity of the corresponding canal in the archeological data (see tables S1 and S2 for a list of symbols and control parameters). These dimensionless weights are based primarily on the mapped width of the feature, which acts as a surrogate for the true hydraulic radius of the ancient and now largely infilled canals. We then used these weights to calculate the relative distribution of flow within the network in the preflood condition: Starting at the source, flow was partitioned at nodes (junctions) in proportion to the relative weights of the edges downstream from the node. This calculation provided an initial flow value for each edge Embedded Image (superscript denotes initial condition, subscript denotes any individual edge i).

We assume that erosion and sedimentation are negligible in this network of low-gradient alluvial channels under stable rainfall conditions and that flow velocity vo is initially constant throughout the network. Under these assumptions, we continue by assigning each edge an initial flow capacity (Embedded Image), which replaces wi as a surrogate for the cross-sectional channel area under equilibrium conditions (Embedded Image). For convenience, we set vo = 1, so that throughout the network, the initial fluid velocity on each edge Embedded Image.

After setting up these initial conditions, we introduced flooding as a gamma-distributed stochastic perturbation that (on average) increases flow rates throughout the network by a factor, γ, of their initial values. This flood scenario increases the volume of water moving through the network. Angkor’s canal and storage network functioned to capture and control surface waters and to minimize destructive flooding of the urban core by providing accommodation spaces for flood water (large reservoirs), fed by a network of excavated canals (30). As such, our model assumes that flood flows within Angkor’s water distribution network were primarily restricted to the artificially incised channels—often further constrained by artificial levees or dikes—rather than occurring as widespread overbank surface flow. We note, also, that damage to the network that can be observed empirically indicates that damage presents as erosion (vertical incision, with few instances of lateral channel migration or avulsion) or sedimentation (infilling of canals particularly in the southern distributary parts of the network) of the canals, rather than as evidence of broadcast flood deposits such as those associated with historical centers in northern Thailand (54, 55). Flood flow velocity, and thus stream power, is higher in straight channels than an equivalent flow volume in naturally meandering channels because the gradient of straight channels is higher (56). Relatively high stream power in straight channels and the inability of the canals to dissipate flood energy to floodplains result in rapid changes in the geometry of incised, channelized networks during peak flow. Our model anticipates this dynamic response via a threshold parameter α: The channel capacity increases if fluid velocity exceeds a threshold Embedded Image, where α∈[0, 1] is a free parameter describing the tendency for the channels to erode.

If canal capacity increases as a result of incision and/or widening during peak flow, then flow capture will create a feedback loop whereby the increased hydraulic radius of the incising canal increases the proportion of flow delivered from the junction immediately upstream of the eroding canal, thus encouraging further erosion. Where flow capture results in a centralization of flow along a single path, rather than as a distributed flow through a drainage network, there is a reciprocal starvation of flow and a decrease in fluid velocity in canals sharing upstream junctions with the incising flow path. This allows sedimentation to occur, thus reducing channel capacity in subsequent flow events. In our model, deposition of sediment due to decreased velocity was implemented via a second threshold parameter β: Channel capacity decreases if fluid velocity on an edge falls below a threshold defined as Embedded Image where β∈[0, 1] is a free parameter describing the tendency for sedimentation to occur.

The increment changes in channel capacity due to sedimentation or erosion depend on both the flow magnitude and the degree to which the threshold condition is exceeded. Because channel capacity adjusts to accommodate changing amounts of flow, the network will stabilize if changes in velocity are insufficient to produce further erosion or sedimentation.

The free parameters determining the extent to which perturbations are amplified in this way are the erosion and sedimentation thresholds (determined by α and β, respectively) and the flood magnitude factor γ. Our model also contains a parameter κ that determines the shape of the distribution from which flow perturbations are drawn. However, this parameter does not markedly alter the average behavior of the network, so we will ignore its influence here. For the results presented in the main text, we set κ = 1, which is equivalent to distributing the perturbations exponentially with the rate parameter Embedded Image.

To examine the model’s behavior as a function of these free parameters, we defined a value Q∈[0, 1] that represents the average factor by which flow on each channel in the damaged network changed relative to its initial value (Eqs. 3 and 4). That is, after the flood waters subside, and the initial flow levels are introduced at source points, a value of Q = 1 represents a condition in which a single path through the network draws all the flow (i.e., the network has been reduced to a single path from source to sink). On the other hand, if Q = 0, the distribution of water through the network is unchanged. This metric provides an unambiguous, well-constrained quantity representing the topological damage to the network. It is important to note that the value does not directly represent the physical damage to the network in the form of erosion and sedimentation. Instead, it quantifies the topological damage to the network’s ability to deliver water to use points.

Sources and sinks

We designated nodes with no input as sources and nodes with no output as sinks. We approximated the watershed feeding the irrigation system as a virtual supersource connected to all physical source nodes. We set the weights of edges connecting the supersource to real sources equal to the sum of all weights exiting the real source nodes to which they connect. Similarly, we connected all sink nodes to a single virtual supersink representing the aquifer, the Tonle Sap Lake, and all human use points (fig. S2). The weights of these edges can be set arbitrarily, since they have no neighbors (here, we set their weights equal to 1). Using this initial topology, we computed the relative equilibrium flow on each edge using our numerical model of flow distribution (fig. S1D).

Flood implementation

In our simulations of flooding, we defined the flood magnitude (γ) as the proportion to which the average rainfall exceeds the total flow into the network at equilibrium. We added this extra flow to the network in the following way: The flow from the supersource increases from its initial value (feq) by a factor of [1 + γ] (that is, ftot = feq + γfeq), and each edge in the network receives a unique perturbation (Xi) to its initial flow fo. In our implementation, X is a gamma-distributed random variable: Xi ~ Γ(κ, θi), where κ and θo = foγκ−1 are the “shape” and “scale” parameters defining the distribution, respectively (note that foγ is the average of the distribution from which the perturbation on each edge is selected).

Cascading failure simulations

In real systems, erosion and sedimentation rates depend on many factors such as the fluid velocity, the concentration of sediment suspended in the fluid, the vegetation in the channel walls, and the material into which the channel is formed. In addition, changes in channel capacity over time would likely be influenced by human activity such as repair works. Our model simplifies this very complex set of conditions to focus on the primary role of fluid velocity in channel erosion and sedimentation. It is worth noting here that our model considers the topology fixed to a specific set of edges and nodes and therefore cannot account for effects such as flood waters breaching the banks of channels or the erosion of new channels. The dynamics can change the weights (capacities) associated with each edge but cannot cause the formation of new edges in the current implementation.

In addition to the assumption of fixed topology, we assume only that fluid flow is conserved and that channel capacity (cross-sectional area) increases from erosion or decreases from sedimentation at a rate that is proportional to both the flow through the channel, and the degree to which the fluid velocity is above or below a predefined threshold for erosion or sedimentation, respectively. After specifying the control parameters and initializing the network to its equilibrium state, the flood flows were added to the supersource and to each of the edges. Then, the flow was distributed through the network, and new edge capacities were calculated from the increased flow velocities. The damage to the network was then computed by distributing the equilibrium flow through the damaged network and calculating Q. The simulation ended if Q converged to a specified degree. If not, the simulation incremented forward (tt′), and the process repeated with the new edge capacities. Here, we identified convergence asEmbedded Image(1)

During the cascade, erosion and sedimentation dynamics were not applied to edges between the supersource and real source nodes, as these represent disparate regions of the watershed. These edges maintain the same flood flows throughout the simulation.


Supplementary material for this article is available at

Supplementary Text

Fig. S1. Conversion of archeological maps of 14th century Angkor into a directed, acyclic network topology.

Fig. S2. A map representing the Angkor network with supersource (S) and supersink (T) included.

Fig. S3. A generic network showing equilibrium flow distribution at initialization.

Fig. S4. To induce damage, random flood flows are added to each edge.

Fig. S5. Final flow distribution resulting from the perturbation shown above, for α = 0.9, β = 0.9, and ΔQmin = 2 × 10−3.

Fig. S6. Results of the model on a generic topology.

Fig. S7. Probability distribution for indegree and outdegree of each node in the final network representation of Angkor’s water distribution network, with supersource and supersink omitted.

Table S1. List of symbols for control parameters accompanied by their purpose in the overall algorithm (role), brief description, and units.

Table S2. List of symbols for noncontrol parameters, accompanied by their purpose in the overall algorithm (role), brief description, and units.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge the support and collaboration of the APSARA Authority, Cambodia. We thank V. Yauw for her work on an earlier version of the model presented here. Funding: This research was supported by The University of Sydney’s DVC Research Strategic Research Excellence Initiative (SREI-2020) project, “CRISIS: Crisis Response in Interdependent Social-Infrastructure Systems” (IRMA 194163), and the Australian Research Council Discovery Project DP170102574. D.E.’s contribution was funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 639828). Author contributions: Conceptualization: D.P. and M.P. Methodology: D.P., C.Z., M.P., J.T.L., and N.F. Software: C.Z., N.F., and J.T.L. Formal analysis: C.Z., N.F., J.T.L., and M.P. Investigation: D.P., C.Z., D.L., D.E., and C.P. Writing, original draft preparation: D.P., C.Z., and M.P. Writing, review and editing: D.P., C.Z., R.F., J.T.L., D.E., C.P., and M.P. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article