Spatiotemporal processes and their implementation in Spatial System Dynamics models

Abstract

Computer-based modelling is indispensable to science. Be it in the field of social or environmental science, modelling is used to learn, to gain new insights into behaviours and to raise awareness. In addition to Cellular Automaton and Agent-Based Modelling, System Dynamics (SD) has evolved as an approach which opened up the field of mathematical modelling to non-mathematicians by establishing a simple Visual Modelling Language. Even though there have been recent developments in expanding SD to include space, Spatial System Dynamics (SSD) hasn't received much attention in geosciences. The typology of SSD introduced in this paper is intended to spur discussions on future possible uses. Models are categorised into three process groups: local processes, diffusion processes and processes that modify space itself. For the latter category spatial analysis may be incorporated into SSD applications. This has been illustrated by means of a model on the co-evolution of surface runoff and landform change.

Keywords::

• Spatial System Dynamics
• model typology
• co-evolution of process and space
• mathematical modelling

1. Introduction

Besides modelling paradigms such as cellular automata or agent-based modelling, SD evolved as a deterministic approach, which was introduced by Jay Forester in the mid-1950s (Radzicki & Taylor 1997). SD enables graphical model design by means of a Visual Modelling Language (VML) which facilitates modelling of complex dynamic systems.

Concepts and techniques developed in SD are frequently applied in the fields of spatial planning, ecology or hydrology (e.g. Forrester 1969; Filgueira et al.2014; Ahmad & Simonovic 2000). However, compared to CA or ABM, SD modelling has not received a great deal of attention in the geoscience community. This may be attributed to its initial use in modelling and simulation of non-spatial processes.

Efforts to make SD spatial resulted in a variety of home-made solutions based on the technical coupling between SD software and Geographic Information Systems (GIS). Moreover, sophisticated software products such as SME, pioneered by Maxwell and Costanza (1997), or SIMILE, introduced by Muetzelfeldt and Massheder (2003), accomplish the full integration of dynamic and spatial processes.

This article aims at revealing the potential of these tools and concepts for modelling spatial processes. Model applications are explored and categorised according to the spatiotemporal process type being modelled. This is intended to answer the following questions:

• What main types of spatiotemporal processes are modelled in SSD?

• How is space conceptualised for these types of spatiotemporal processes?

Apart from this general overview, we primarily focus on a specific process type which is not only influenced by its spatial environment, but also modifies space itself (cf. BenDor & Kaza 2012). The capabilities of SSD for modelling this type of process-space interaction are demonstrated by means of a hypothetical model on surface runoff, erosion and landform evolution.

In the following section SD and the incorporation of space in SD are briefly explained. This is followed by the introduction of a typology on SSD including a general description of process categories, approaches for space-time integration and practical application examples. The paper concludes by revisiting the two primary investigation questions.

2. System Dynamics

SD is an approach for modelling the behaviour of complex systems over time. SD models consist of differential equations which are formalised using a VML. The differential equation may be defined by making assumptions about the system being modelled. It can, for instance, be assumed that population growth is proportional to the size of the population (Blanchard et al.2002). This rate of change can be expressed by

(1)

where P is the population and k is a proportionality constant. The solution to this equation is an exponential function (primitive integral) which describes the system behaviour based on the simple model assumption of proportional population growth.

This mathematical model can be formalised in an SD model of stocks and flows (see Figure 1A). Flows (double line arrrow) represent the rate of change over time. Stocks (rectangle), also called state variables, describe the state of the system over time (population at a certain time) in the simulation. This means the stock provides the primitive integral as a result.

Figure 1 Stock and flow (a) and causal loop diagram (B) of population growth

The same example may also be represented by means of a conceptual diagram showing causalities in the system (see Figure 1B). Existing causal relations between system elements are indicated by an arrow with positive or negative polarity. A negative polarity indicates that the change of a variable in one direction causes the second variable to change in the opposite direction. A positive polarity, however, means that variables change in the same direction.

Causal loop diagrams help reveal feedbacks, which are a central feature of complex non-linear systems. The two types of feedback can be classified as positive and negative feedback. Positive feedback means that growth reinforces further growth (e.g. higher population leads to increased potential for further growth). In contrast, negative feedback acts to negate this growth (Ford 2010). For instance, a large population also has an increased number of deaths, which diminishes the population.

The basic idea of using a VML is to make mathematical modelling accessible to non-mathematicians and to facilitate collaborative and interdisciplinary work. Furthermore, SD models are aiming to improve knowledge of a systems's behaviour to reveal patterns, expose sources of undesirable system behaviour, and help avoid unforeseen consequences of future policy implementations. Practical implementations of real systems may consist of a large number of stocks, flows and feedback loops, which often leads to unforseeable model results.

Systems become even more complex if space is considered in SD models. In an SSD model stocks and flows are bound to locations. This makes cause and effect not only distant in time, but also distant in space (Sanders & Sanders 2004).

3. SSD typology

There are a variety of ways to organise models in a typology. We decided to base the typology in this article on spatiotemporal process groups, since its intention is to reflect on the capacity of SSD to model different types of processes. A literature review showed that the majority of SSD models are used to model either local processes (e.g. tree growth in a forest), diffusion processes (e.g. water flow or urban spread) or processes dynamically changing the underlying spatial structure (e.g. LULC change).

A similar organisation was proposed by Miller et al. (2005) for dynamic GIS models. Furthermore, Mazzoleni et al. (2006) also distinguished between SSD models on local processes and lateral processes, which correspond to the first two classes being used in the presented typology. A more conceptual view on spatial processes was proposed by Brown et al. (2005) with the differentiation between Eulerian and Lagrangian processes. The first view describes processes which influence properties at fixed locations (e.g. temperature). A Lagrangian perspective, however, tracks the changing location of particles through space (e.g. nutrient flow). This underlines the difference between diffusion processes which are just dependent on spatial properties and those processes which actually alter spatial properties. The typology classes 2 and 3 correspond to these categories.

In some cases the boundary between categories is not sharp and needs further explanation. Moreover, other distinctive features are described with respect to the three main classes. Further ancillary model features being discussed are implementation techniques such as spatial data models and software usage.

Local processes

In instances where spatial neighbourhood interaction is assumed to be irrelevant for the processes being modelled, spatial interaction is occasionally ignored in the model. This may, for instance, be an appropriate approach for modelling plant growth, if enough space and resources are available and thus spatial competition among plants is negligible.

A model of this process type was introduced by Scheffran et al. (2007) to better understand the dynamics of biofuel crop growth in Illinois. In their SSD model the landscape is decomposed into townships (agents), which are treated independently in the simulation. Whereas geographic position is important due to spatial differences in crop yields, neighbourhood dependencies can be ignored, since product competition takes place on a scale well beyond the township level. However, single agents are structured vertically into a crop yield model and a profit model. These submodels represent local processes which are interrelated and thus have to be coupled in the implementation. The modular structure of SD facilitates an implementation of this model type.

`In order to enable a spatial visualisation, the SD model needs to be linked to a spatial data model. In practice, formats such as raster or vector data provided with GIS are used. The association between SD and data model is established by replicating SD to assign stock variables to single cells of a raster grid (see Figure 2) or to point, line and polygon geometries of a vector data model. Since spatial interactions are being ignored, the data model does not require any information concerning relative positions (topological information).

Figure 2 Local evolution process type in SSD

In order to enable communication between SD and GIS data, these software modules may be coupled. For instance, Lowry and Taylor (2009) developed an SSD tool by coupling the SD software Powersim to Google Earth (GE). The simulation is entirely conducted in Powersim, whereas GE is used to visualise outputs. The link between those applications was established by relating stock values to colour attributes of the Keyhole Markup Language (KML) file in GE.

Diffusion processes

Diffusion models are typically intended to model lagrangian processes. This means particles or any other elements which change their spatial positions over time. In SSD this is expressed as flows between stocks. The assignment of stocks to spatial data is similar to the local process type (see Figure 3). However, in contrast to models on local processes, the relative position of cells (or any other geometric primitives of spatial data models) to each other is of particular significance. The SD model defines which cells are linked and where direct exchange can take place. The spatial data model determines their spatial relationship.

Figure 3 Diffusion process type in SSD

The diffusion process type is commonly characterised by rapid flows which require high temporal model resolutions. Modelling of rapid processes is especially relevant in the field of hazard risk assessment. For instance, Zhang (2008) simulated nitrobenzene concentrations in the Songhua River after a pollution accident. Similar to this implementation, hydrological discharge simulations are often linked to other processes such as nutrient cycling or organic decomposition (e.g. Voinov et al.2004; Reap et al.2003). This in turn requires a vertical integration of submodels to consider relationships between diffusion processes.

Besides physical diffusion processes, a number of models on the spread of animal species, infectious disease and human migration have been implemented in SSD. As in the previous example, systems being modelled are often composed of related diffusion processes. This also holds true for a model on fox rabies spread in Illinois implemented by Deal et al. (2004). In short, the number of foxes infected is modelled per raster cell as a function of fox densities and bite rates. The disease spreads as foxes emigrate as a result of limited resources (prey) in their host cell.

In contrast to animal spread, human movements are typically bound to traffic networks. Whereas the topology of a raster grid is characterised by regular neighbourhood dependencies and uniform distances between cells, a network model exhibits more complex neighbourhood relations. This is often bypassed by using euclidean (air-line) distances instead of network distances (e.g. Haller et al.2008) or networks are merely used to indicate relational space (e.g. Cruz & Olaya 2008).

Regardless of the spatial data models being used, diffusion processes in SSD place high demands on model developers. The simulation of diffusion processes not only requires an association of SD and spatial representation, but also the integration of flows between spatial locations. Two approaches for handling this established in SSD are:

• Stock array functions are provided with SD software packages (Sanders & Sanders 2004). The array is used like a two-dimensional matrix where each element is a separate stock linked by flows to neighbouring stocks (see Figure 3). One disadvantage of this approach is that advanced programming skills are required to couple the array to a spatial data model.

• Alternatively, specialised SSD software may be used to duplicate and link models developed as conventional SD models. This also facilitates the visualisation of outputs in a GIS-like map, since an interface to stand-alone GIS software is commonly foreseen or GIS data models may even be integrated in the SSD software (e.g. SME or SIMILE).

Structural change processes

Structural change processes modify the qualitative structures of space. This may be a shift in agricultural land use due to changes on the selling market or a temporarily closed road due to a car accident. Fawcett and Oliveira (2000) also mentioned this type of causality in the context of patient treatment after an earthquake disaster. Transport of injured people and evacuation may cause congestions. Moreover the earthquake may have damaged bridges and tunnels and thus modified the trafficable street network.

In many cases the relationship between processes causing a change in spatial structures and the structure as such is bidirectional. For instance, Birkin (2013) explained this mechanism based on the co-evolution of population and infrastructure networks. He concluded that this system is reasonably modelled by evaluating not only the influence of future population on the need for a new rail network, but also the effects of this new transportation network on population movements.

This type of feedback is often composed of a lagrangian diffusion process and an eulerian process of local change (see Figure 4). Sanders and Sanders (2004), for instance, modelled migration flows between regularly sized city zones and the city environment as a function of relative attractiveness. The movement of individuals, however, also affects local conditions due to considered effects of migration on the housing market. Thus, the lagrangian migration process has an impact on urban infrastructure and vice versa. This means that eulerian and lagrangian processes are combined and also interact in this system.

Figure 4 Structural change process type; land use change example

A model on invasive bark beetle spread introduced by BenDor and Metcalf (2006) exhibits a similar mix of considered process types. Bark beetles migrate dependent on population density (beetles per bark surface area) and distance between origin and destination raster cell. Tree die-off rates in a raster cell are directly affected by the number of ash borer larvae at this location. Beetle dispersion is in turn affected by carrying capacities per raster cell (bark surface area of host trees), which is reduced due to tree deaths. This closes the feedback loop between the lagrangian diffusion process (beetle migration) and the eulerian process (beetle-induced tree deaths).

In the models introduced so far spatial structures were expressed as quantity per raster cell (area of bark surface per cell) or zone (e.g. housing per zone) in the models. In other words, the models do not depict houses or trees as independent objects as they exist in the real world, but represent their properties in an aggregated way. A discrete representation of space is, however, reasonable when the relation between individual objects is of relevance. What makes the implementation of this type of model challenging is the fact that SD is an aggregated approach by definition. Schieritz and Milling (2003) used the metaphor of modelling the trees versus modelling the forest to put emphasis on the distinctive characteristics of ABM and SD modelling.

Nevertheless, the discrete object view may also be involved to a certain degree in SSD models as, for instance, demonstrated by Lauf et al. (2012) in a model on land-use dynamics. Their application consists of a supply and demand model implemented in SD, which is coupled to a raster map showing changes in land-use. The supply and demand model operates on a macro level and drives the growth of land-use objects. Even though space is discrete in this example, no feedback between spatial structures and supply-demand processes is considered.

A first step towards an integration of this feedback is presented by BenDor and Kaza (2012). The spread of an infectious disease is dependent on spatial connectivity in their model. Connectivity between locations is represented by means of a network. The peculiarity of this network is that links can be switched on and off to simulate a quarantine policy. Links are eliminated if the infected population at a node exceeds a certain ratio, and links are then restored if infection rate is below a certain threshold. In this way, the process of disease spread changes network structures, which in turn drastically alters the dispersion of the disease.

This system may be modelled more realistically by considering not only binary conditions (link is on or off), but also network distances. It can be assumed that the spread from one settlement to another depends on spatial proximity. In such an application GIS is needed not only for the purpose of visualisation, but also for calculating network distance. In the following section, the use of GIS analysis in SSD is demonstrated by means of a model on the co-evolution of surface runoff and hillslope geometry.

4. A SSD model on runoff-topography feedback

The reallocation of materials by flowing surface waters leads to temporal evolution of the topography. Material is eroded at one location and accumulated at another as a consequence of differences in spatial flow characteristics. The topography modified in this way affects magnitude and direction of subsequent runoff events and thus influences further material transport.

The method being used is equivalent to the diffusion model type presented unded ‘Diffusion processes’ in section 3. Water flow and erosion was modelled in two-dimensional space based on stock array to raster data association (see Figure 3). Water flow from stock to stock was modelled in x and y directions as a function of water stored in the stock and slope incline derived from a DEM (see Figure 5A). An additional factor was used to scale flow rates. Then water flow rates served as an input to the erosion model represented in a second stock array. The erosion model calculates material flow between stocks as a function of material present in a stock and water flow rates (see Figure 5b). As in the runoff module, an additional scaling parameter was used to adjust material flow.

Figure 5 Association of (A) water flow and (B) erosion simulation; 2D stock arrays are represented schematically

The flow rate from one stock [a] to a neighbouring stock [b] can be mathematically described for the water flow (equation 2) and erosion modules (equation 3) as

(2)

(3)

where F and V are flow rates and volumes of water (subscript w) and material (subscript m) stored in the stock. Flows are scaled based on the parameters p_m and p_w. Furthermore, slope S indicates the topographic incline towards the neighbouring stock and determines the flow direction. In order to take the evolving topography into consideration, the underlying DEM is updated with erosion data after each rainfall-runoff event.

The model was implemented by coupling Vensim to ArcGIS using a middleware programmed in Python. Functions available with the desktop software can be imported in the code to successively run the required SD and GIS calculations as follows. (1) Surface runoff and material flow is simulated for the time of one rainfall event. (2) The SD software is interrupted to update the DEM with changes in material dispersion in GIS. (3) Slope incline is calculated in GIS based on the updated DEM. (4) The SD software is restarted to simulate the next rainfall-runoff event. These steps are repeated in a loop for a predefined number of rainfall-runoff events.

System behaviour

A north-facing slope was chosen to demonstrate the spatial feedback involved in the rainfall-runoff system. The topography was represented in an 8 × 8-cell raster grid with a resolution of 10 m (see Figure 6A). Due to the high relief intensity of about 30 m, high runoff and erosion rates can be expected. The model was run for multiple rainfall events with duration of 10 hours and a rainfall intensity of 5 mm/h.

Figure 6 Evolution of resultant runoff vectors derived for runoff-erosion event numbers 1, 100 and 500

Resultant runoff vector fields indicate significant changes in topography over time (see Figure 6B,C,D). Material is eroded on upslope areas and accumulated at the footslope, which overall reduces slope inclines and flow rates. This is especially apparent in runoff event 500, which shows a sudden change of flow rates between upslope and footslope sections from relatively high to very low flow rates. This is due to the slope concavity, which causes a deceleration of flow velocities (Fullen & Catt 2004). Accumulation at the footslope and retrogressive erosion leads to the formation of a sharp boundary between a slightly inclined slope and almost flat areas. Even though this may be a consequence of ignoring riverbed erosion in the model, similar slope characteristics can also be observed for hillslopes under low riverbed erosion conditions (Zepp 2008).

The evolution of topography also affects flow directions. Convergent or parallel vectors as modelled for the first rainfall event gradually change to divergent flows. This holds especially true for the footslope section, where accumulation causes lateral runoff. Thus, rainfall-runoff events reshaped the terrain geometry from linear-linear to slightly concave-convex.

Furthermore, erosion is non-linear in this system. The steep slope at the beginning of the simulation causes rapid material transport and landform changes. However, erosion processes are continuously mitigated as slope inclines decrease. Even though the topography was levelled significantly as a result of 500 rainfall events, the evolution from steep slope to flat terrain would require a substantially higher number of runoff simulations. This means that spatial feedback exhibits a negative polarity (more runoff →  more erosion →  decreasing slope incline →  less runoff).

5. Conclusion

In this article we suggest that SSD models can be assigned to three main process groups: (1) local processes without neighbourhood interaction, (2) processes of diffusion and (3) processes which modify spatial structures. Table 1 shows the summary of models being reviewed structured according to this typology.

Table 1 Selection of SSD models

#	Type	Field of application	Implementation	Authors
(1)	Local	SSD methodological study	SD/GIS coupling	Lowry and Taylor (2009)
(2)	Local	Bioenergy farming	SSD software^a	Scheffran et al. (2007)
(3)	Diffusion	Flood damage simulation	SD/GIS coupling	Ahmad and Simonovic (2004)
(4)	Diffusion	Infectious disease spread	SSD software^a	Deal et al. (2004)
(5)	Diffusion	Casualty treatment after earthquake	–	Fawcett and Oliveira (2000)
(6)	Diffusion	Nitrogen accumulation	SD/spreadsheet coupling	Ford (2010)
(7)	Diffusion	Domestic migration	SD software^b	Haller et al. (2008)
(8)	Diffusion	Life Cycle Assessment	SSD software^a	Reap et al. (2003)
(9)	Diffusion	Ecosystem modelling	SSD software^a	Voinov et al. (2004)
(10)	Diffusion	Water pollution	SD/GIS coupling	Zhang (2008)
(11)	Structural change	Infectious disease spread	SD software^c	BenDor and Kaza (2012)
(12)	Structural change	Specie spread	SSD software^a	BenDor and Metcalf (2006)
(13)	Structural change	Land use simulation	SD/GIS coupling	Lauf et al. (2012)
(14)	Structural change	Urban migration	SD software^b	Sanders & Sanders 2004)
(15)	Structural change	Urban residential development	SD/GIS coupling	Xu (2011)

a SD and GIS modules are compiled by means of the Spatial Modelling Environment (SME).

b Euclidean distances are used as input variables in the SD software (no spatial module).

c The model uses relational space (network), instead of geographic space.

In general, a spatial view on processes is generated by duplicating SD components and assigning them to elements in a spatial data model as provided with GIS software. In models of process type 1 each SD component operates independently. Even though position is of relevance, there is no need to consider the position of elements relative to each other, since there are no neighbourhood interactions. Nevertheless, interactions may take place vertically between different local processes.

A vertical overlap of processes may also be given for type 2 (e.g. interrelation of surface runoff and erosion). The modularity of the SD approach meets the requirement for modelling and linking such processes and thus facilitates model-building. This vertical link is complemented with spatial links between single SD components for modelling process type 2. Furthermore, information on relative position is needed to consider effects of neighbourhood on processes.

In contrast to type 1 and 2, space is not static in process type 3. In this most complex process category spatiotemporal processes modify spatial structures, which in turn affect process behaviour. Spatial structures considered in a model on type 3 may be the length of a street, the area of a land-use patch or the geometry of the earth's surface.

In order to quantify these structures and to use them as feedback in the model, the spatial analysis functionality of GIS may be adopted into SSD applications. So far, the use of GIS in SSD applications has usually been restricted to visualisation and data management, since GIS's spatial analysis functionality is not necessarily needed for spatiotemporal processes typically modelled with SSD.

The use of spatial analysis in SSD is demonstrated by means of a hypothetical model on flow processes and landscape evolution. The successive change in slope incline is analysed throughout the simulation to provide feedback to runoff and erosion processes. This example reveals non-linearity imposed by negative feedback between surface runoff, erosion and landform evolution and also underlines the value of incorporating spatial analysis into SSD.

Acknowledgements

The authors appreciate the financial support provided by the Austrian Science Fund (FWF) through the Doctoral College GIScience (DK W 1237-N23) to the first author and the Natural Sciences and Engineering Research Council (NSERC) of Canada CGS doctoral scholarship to the third author.