Estimating Spatio-Temporal Risks from Volcanic Eruptions Using an Agent-Based Model

: Managingdisasterscausedbynaturalevents,especiallyvolcaniccrises,requiresarangeofapproaches, including risk modelling and analysis. Risk modelling is commonly conducted at the community/regional scale using GIS. However, people and objects move in response to a crisis, so static approaches cannot capture the dynamics of the risk properly, as they do not accommodate objects’ movements within time and space. The emergence of Agent-Based Modelling makes it possible to model the risk at an individual level as it evolves over space and time. We propose a new approach of Spatio-Temporal Dynamics Model of Risk (STDMR) by integrating multi-criteria evaluation (MCE) within a georeferenced agent-based model, using Mt. Merapi, Indonesia, as a case study. The model makes it possible to simulate the spatio-temporal dynamics of those at risk during a volcanic crisis. Importantly, individual vulnerability is heterogeneous and depends on the characteristics of the individuals concerned. The risk for the individuals is dynamic and changes along with the hazard and their location. The model is able to highlight a small number of high-risk spatio-temporal positions where, due to the behaviour of individuals who are evacuating the volcano and the dynamics of the hazard itself, the overall risk in those times and places is extremely high. These outcomes are extremely relevant for the stakeholders, and the work of coupling an ABM, MCE, and dynamic volcanic hazard is both novel and contextually relevant.


Introduction
. During natural disasters, the level of hazard varies over space and time, due to various factors. During volcanic eruptions, for example, the spatial extent of the impact is controlled by the contrasting nature of the volcanic sources, the type and magnitude of the explosive eruptions, and the local topography (Lirer et al. ). Likewise, with regards to floods, the specific hazard can vary depending on hydro-meteorological aspects and drainage basin topography (Merz et al. ; Yan et al. ). Such spatio-temporal variability of hazards means that the associated risks can be dynamic over time and space. This is especially true during volcanic crises, where the time and location are critical in defining the risk to human populations; 'being in the wrong place at the wrong time' can lead to fatalities. Under certain conditions, those residing in the vicinity of a volcano may need to evacuate hazardous areas quickly. This is o en the only way to reduce the risk from a volcanic impact because it is almost impossible to survive the hazardous material emissions that occur during an eruption, such as pyroclastic flows and toxic gases (d'Albe ). .
The combination of the mobile nature of people together with the spatio-temporal variability of hazards means that risks are dynamic. Better spatio-temporal modelling is required to estimate these dynamics rather than relying on a static risk map (Alcorn et al. ). Risk modelling is important in understanding the potential ABM is used here to account for the nature of the social processes in an emergency situation, which are complex (Dash & Gladwin ). Representing human behaviour in such a situation is extremely challenging due to the di iculties associated with modelling human behaviour. Specifically, the responses and behaviour of people during a disaster will vary according to their socioeconomic and demographic characteristics (Christia ; Donovan ; Dove ; Lavigne et al. ; Rianto ). This integration is conceptually able to model the dynamics of human risk to natural hazards in a spatio-temporal manner. .
The risk is considered here to be the probability of harmful consequences or expected losses that result from the interaction between hazards and vulnerable people or objects (Blaikie et al. ; UNISDR ). Risk is estimated to be a function of hazard and vulnerability using Equation (Sar et al. ; UNISDR ). Consequently, when the value of the hazard changes, the risk will change as well. For example, consider a population with a degree of vulnerability arising from living in the certain hazard zone (Figure ). During a volcanic emergency, the magnitude and intensity of the hazards vary with respect to the proximity to the summit, as well as the topographical conditions that determine the direction of the flow of volcanic material. As the population may be moving during the crisis, their hazard level will vary over time (t 1 to t 2 ). Simultaneously, the hazard will vary due to the changing intensity of the volcanic activity. Therefore, the degree of risk varies in both cases. The aforementioned concept of risk with regards to the mobile nature of people and the dynamics of hazard can be used to formulate the spatio-temporal risk model on an individual basis. . Individual risk (Newhall & Hoblitt ) is the probability that a particular individual, at known coordinates, will be killed or injured by the volcano within a specified period. In this research, we specify the hazard as a potentially damaging eruption that may cause loss of life, injury, and social disruption (UNISDR ). A map of the hazard can consist of several elements. Specifically, in the case of volcanic hazards, these elements include the types of hazardous material that are emitted during the eruption, such as the Pyroclastic Density Current (PDC), lava flow, and tephra fallout (Alcorn et al. ). In the study area, the historical records have been compiled to form a single hazard map (BNPB ). Meanwhile, we describe vulnerability as the characteristics of a person or group that influence their capacity to anticipate, cope with, resist, and recover from the consequences of a natural hazard (Blaikie et al. ). Vulnerability is multidimensional in nature and can be measured using a combination of many di erent variables (Lummen & Yamada ). Here, we quantify vulnerability using the Social Vulnerability Index (SoVI) (Cutter et al. ). The SoVI has been developed based on several attributes; namely, socioeconomic status (income, political power, prestige), gender, race and ethnicity, age, commercial and industrial development, employment loss, rural/urban, residential property, infrastructure and lifelines, renters, occupation, family structure, education, population growth, medical services, social dependence, and special needs populations (Cutter et al. ). This index has been widely used to measure social vulnerability regarding environmental hazards in various regions (Armas & Gavris ; Chakraborty et al. ; Garbutt et al. ; Letsie & Grab ; Schmidtlein et al. ; Tapsell et al. ; Yoon ). .
Individual risk assessment has di erent criteria from risk assessment for the community/region. A risk assessment on a regional scale, using a community or group of people as the assessment unit, uses criteria based on the characteristics of the community and the region where the population lives, although there is no literature describing risk assessment on an individual basis. Therefore, in this research, we define several characteristics of individuals (or 'agents') to generate the SoVI and assess the hazard degree of their specific location to provide the value of the hazard at a particular time. The concept used here to define individual risk consists of three parts: defining the socio-economic parameters for individual vulnerability; defining the hazard at the individual location; and measuring the risk. There are several socioeconomic parameters that are used in the vulnerability assessment through the application of an MCE. Meanwhile, the hazard is assessed based on the person's location within the given hazard zone. The risk is defined based on the measured vulnerability and the hazard level ( Figure ). . Thus, MCE is o en used to analyse the compromises which exist between multiple choice alternatives, while GIS provides the ability to analyse complex spatial problems from several layers of spatial data. MCE analysis starts with the construction of an evaluation matrix that contains elements that reflect the characteristics of the set of alternatives, based on a specific set of criteria. Each element can be weighted based on its contribution to the goal using several techniques; namely, ranking, rating, pairwise comparison (AHP), and trade-o (Malczewski ). Commonly, two techniques are employed to aggregate this element so that the final result can be achieved through a weighted linear combination (WLC) and Boolean overlay (Eastman ; Malczewski ).
. Finally, those concepts need to be implemented in an ABM to simulate the dynamics. Although potentially powerful, and successfully integrated into GIS (Carver ), the integration of MCE and ABM is rare. Bishop et al. ( ) discussed the potential use of MCE in exploring the various outcomes of decision-making processes by human agents with di erent preferences, while (Gao & Hailu ) applied MCE in an Agent-based Simulation as a Decision Support System regarding the selection of fisheries management strategies. However, we found no articles using such an approach for disaster risk modelling. ABM has emerged as a valuable alternative to the traditional aggregate mathematical modelling, as it can accommodate the complexity of a system through its ability to capture the interactions between agents on the same or di erent scales (Gilbert & Troitzsch ). While an ABM consists of discrete 'agents' that can interact within an environment (Gilbert ) and has the ability to incorporate the complex, multiple attributes of individuals, it lacks a method for evaluating those attributes into a single decision/value. MCE, on the other hand, appears promising with regard to solving this problem, and so has the potential to be integrated into the model.

Materials and Method
. The basic concepts of integrating multi-criteria evaluation (MCE) into Agent-based Modelling (ABM) for the Spatiotemporal Dynamics Model of Risk (STDMR) were theoretically discussed in the previous section. In this section, we provide an overview of the application of STDMR. It starts with an introduction to the study area, an outline of the process of collecting data from this area for use in the model, and the integration of STDMR to the ABM of evacuation. This ABM is developed based on the spatio-temporal volcanic evacuation model framework (Jumadi et al. ) that has been improved with the individual evacuation decision model and subjected to validation using real data (Jumadi et al. ). An overview on how the ABM works and integrates the STDMR is provided below.

Study area .
Merapi Volcano is located on Java Island, Indonesia, and poses a potential hazard to the surrounding communities. Recent work suggests that the potential for Merapi erupting is far higher now than has been found historically (Andreastuti et al. ; Camus et al. ). These risks were confirmed by the most recent event, that occurred in . This potential eruption style varies between either Sub-Plinian or Plinian. In disaster studies, volcanic explosivity indices (VEI) are o en used to describe the potential destructiveness of an eruption. The VEI ranges from (least destructive) to (most destructive) (Newhall & Self ). The VEI of Merapi's eruptions ranged from -but unexpectedly increased to in (Surono et al. ). As a consequence of the VEI eruption in , the area surrounding Merapi su ered the worst disaster of the last century.
. Eruption activities commonly produce diverse hazardous events, and Merapi has unique hazard types. Nue'es ardentes (Bardintze ) and lahars (Lavigne & Thouret ) are two particular hazards that are harmful to the communities besides the explosion (Thouret et al. ) and the slight impact of ash (Damby et al. ). Nue'es ardentes is produced by gravitational dome collapse, whereby the impact's extent is commonly controlled by its topographical character (Kelfoun et al. ). It can travel up to . km from only a few individual events (Abdurachman et al. ), whereas lahars events are usually initiated by high rainfall intensity (Lavigne & Thouret ). Lahars are the overbank flows of pyroclastic material coupled with rainwater and are considered the most dangerous part of the material flow system in Merapi (Charbonnier & Gertisser ). The direction of this flow is strongly influenced by the initial flow direction and the topography, that a ect the subsequent spatial extent of the hazard (Itoh et al. ). .
The eruption in strongly a ected the geomorphological structure (Saepuloh et al. ) and geological character (Gertisser et al. ) of Merapi, which has implications for the spatial extent of the hazard map (BNPB , ). Also, the eruption changed the potential direction of the pyroclastic or lahar flows. It can be predicted that the southern flank of Merapi will probably be more heavily impacted by the next eruption than elsewhere ( Figure ). Based on this, we use the Sleman Regency, an administrative region located on the southern flank of Merapi, as a case study. Sleman ( Figure ) is geographically located between • ' " and • ' " east longitude, and • ' " and • ' " South latitude. Sleman covers , hectares/ . km ) or about % of Yogyakarta Province. Administratively, this region consists of districts, villages, and , hamlets.

Data collection
.
The data used here are collected from primary and secondary sources. The secondary data consist of admin-istrative boundaries, volcanic hazard zones, shelter locations, land use, census microdata, and road networks. The primary data are drawn from an extensive questionnaire survey undertaken in . The survey data were used to complete the variables of the census microdata (Center ) in order to create the population of human agents as well as develop the evacuation decision-making process (Jumadi et al. ). .
The questionnaire was developed to gather information regarding the mechanisms used in decision-making and the interaction of people during eruption crises in the Mt. Merapi region. A literature review was conducted to explore the variables that influence the decision-making and interaction. Five primary variables were assessed via the questionnaire survey; namely, socio-demographic characteristics, perception of volcanic hazard, decision-making behaviour, interaction during a crisis, and willingness to accept and act on an alert. The question list was developed based on these variables. The demographic characteristics are used in this research to characterize each agent as well as identify its social vulnerability. Meanwhile, the decision-making process takes place when people start to evacuate. It explores the variability of the population in terms of making decisions during a crisis. The main indicator of this behaviour is the start time, related to the onset of the enhanced activity of the volcano. The questionnaire explores the factors that might motivate or demotivate people when making decisions about whether or not to evacuate. The interaction is also taken into account on the questionnaire. In this case, data on the probability that people will forward their information about alerts and the associated impacts on their decisions are needed.
. In order to collect these variables, stratified random sampling was applied. Household member samples, represented as building units, were selected randomly for each building block (dusun). This area segmentation is based on the fact that each dusun has one village chief who mobilises people (Rukun Tangga) and, commonly, in the rural areas of Indonesia, has homogenous social characteristics. Twelve villages were selected within a radius of km. Several ring bu ers with distance ranges of km were created to define the sampling areas, with three villages selected from each range. Furthermore, ten participants from each village were selected randomly, resulting in participants in total.
. The results of the survey were statistically analysed to develop the evacuation decision model (Jumadi et al. , see). The data from the survey were tabulated and analysed using SPSS. The result of the regression analysis is used to develop the driving forces governing the decision to evacuate or stay that is implemented as the decision-making rule of the agent in the ABM.

Agent-based STDMR of volcanic evacuation
General framework .
The ABM of the volcanic evacuation simulation was developed based on the relationship model between the volcano and the surrounding population. The basic model and its complete documentation is provided in a previous publication (Jumadi et al. ). This model marries the physical environment and social interactions to generate the value of risk (Cova ; Pons et al. ; Sengupta & Bennett ), as presented in Figure . The physical variables are generated from the characteristics of the volcano and its hazard zone, including the VEI, activity length, activity level, and the spatial extent of the hazard (Figure ). The VEI can also be used to estimate the spatial extent of the impact, and is generated based on the eruption record of the volcano. The probability distribution of the VEI is used to define the VEI in the simulation. Eruption records can also be used to estimate the length of the crisis (activity length). Considering that volcanic activity fluctuates within a period of crisis, the activity level from the rest condition (i.e. period of inactivity) to the climax of the event can be estimated. This activity level is also related to the hazard. Figure : General Framework of the ABM. The le -hand box (Physical Variables) illustrates how the VEI and activity level are used to estimate the spatial extent of the volcanic hazard. The VEI and the length of the activity are the physical characteristics of the volcano as recorded in the literature. These feed into the socioeconomic variables (right-hand box) that are the attributes of the human-agents and are used to assess the overall risk. The hazard (which varies spatially and temporally) is used to estimate the exposure to the population and, subsequently, the overall risk. The activity length is used to estimate the evacuation time (the period during which the risk changes dynamically as a result of people's movement) and, subsequently, the spatiotemporal risk dynamic which quantifies the risk at every time period (t).

Synthetic population generation .
To create human-agents, census microdata (Center ) and a separate survey of households are used. We use conditional probability (Monte Carlo Simulation) to generate the synthetic population of agents (Heppenstall et al. ; Moeckel et al. ), a complete description of which is provided elsewhere (Jumadi et al. ). In this model, human-agents are generated for each sub-district of Sleman in individual units, grouped as households. It would be advantageous to create an agent to represent every person, but the AnyLogic PLE limits the total number of agents to , . Therefore, approximately , ( % of the total population of study area) representative agents were created. The characteristics of the people are used in order to calculate the SoVI variables; they, together with the other variables, influence the decisions of the agents (for further details on this, see Jumadi et al. ). Each agent is initialised with characteristics that are randomly drawn from probability distributions generated from the census microdata and survey data.

Model Description
. The model is documented using Overview, the Design Concept, Detail (ODD protocol) (Grimm et al. ; Polhill ), as presented in the following sub-sections. Further details are provided at: https://tinyurl.com/ stdmr.

Overview
Purpose .
This model is designed to demonstrate how the risk varies over a spatio-temporal range of volcanic activity.
Here, the risk dynamically changes along with the changes of hazard due to the di erent volcanic activity levels.
Entities, state variables, scales, and environment .
The model consists of three agents, namely, the volcano, stakeholder, and people (population), that interact within the geographical environment (Jumadi et al. ). The volcano acts as an agent, which initiates the hazardous situation and influences the environment as a potential threat to the surrounding population. The other agents in the interactions are the stakeholder and the population. The stakeholder, in this case the authorities (government), plays a role in observing and analysing the activities of the volcano and in issuing warnings  to the population. In the ABM simulation, each agent displays specific behaviour and mechanisms when interacting with others as well as with the environment. The environment is represented through spatial data with dynamic hazard properties. An overview of the attributes of the model is provided in Table (extended from  Jumadi et al.  ).
Process overview and scheduling .
The human agents in this model are the focus of the observation. They should act to protect themselves from danger in a crisis. Therefore, the human-agents utilise a decision mechanism that allows them to take the decision to evacuate in dangerous conditions. This evacution decision is based on the Normal -Investigating -Evacuating state model that is provided elsewhere (Jumadi et al. ). Conceptually, this evacuation decision is influenced by several factors, including: ( ) risk communication and warning; ( ) perception of risk; ( ) community and social network influence; and ( ) disaster likelihood, environmental cues, and natural signal (Ahsan et al. ; Dash & Gladwin ; Lim et al. ). The mechanism of the individual decision during evacuation, based on the literature review, is provided in Figure . The transition between the di erent states is based on threshold-based rules (Kennedy ; Robinson et al. ) that evaluate the strength of the force to evacuate based on various factors. .
The risk to individuals is evaluated based on the hazard and vulnerability variables. The hazard level is measured by the environmental properties in the agent's location, whereas the vulnerability of individuals is based on the SoVI, which is itself evaluated based on its socio-demographic character. The following subsection elaborates on this risk model in detail. The risk to the individual might change a er the decision to move is made since their location changes. When people move during an evacuation process, the hazard of their environment will change, resulting in the risk changing also. Therefore, the value of the risk is dynamic over time as an individual moves.

Design concepts .
The following concepts will be used in this model: . Emergence: the population risk patterns emerge from the behaviour of individuals in making the decision whether to evacuate or not. The individual risk is analysed using the SoVI attributes (Chakraborty et al. ; Cutter et al. ) and the individual's environmental hazard susceptibility (BNPB ).
. Sensing: Every individual in the population understands its position which they use to analyse the hazard.
. Interaction: there are various interactions in the model: (a) interaction between the VEI with the spatial extent of the hazard; (b) interaction between the hazard and the people; (c) interaction between the volcanic activity and the alert level of the stakeholder (the alert level will increase with increasing volcanic activity).
. Stochasticity: the characteristics of the volcano, such as its VEI, type of hazard and spatial extent of hazard zone are summed up based on the literature.
. Observation: the spatial pattern of risk of people will be recorded at the end of each experiment. These data will be further used in spatial analysis using GIS.

Details
. Based on the concept of the individual risk model (Section ), risk comprises two main components: hazard and vulnerability. Here, we provide a calculation procedure based on MCE to implement this into the ABM of evacuation. Consequently, Java functions are designed based on this concept and implemented into the previous model (see Appendix A) (Jumadi et al. ).  We implement a dual hazard model to set up the environment of the agent-based evacuation model ( Figure  ). The first hazard model (a) is the actual spatial extent of the hazard for Merapi, based on several historical records of eruptions, including the one that occurred in (actual hazard) (BNPB ). The distribution of hazards in this map is strongly based on the physical record present in volcanic material deposits, while the second hazard model (b) is the one used by the government to alert the population at risk during the eruption (perceived hazard) (Mei et al. ). This hazard model is a rough estimation based on the distance from the volcano, as well as being closely related to the administrative boundaries. This makes it easier for the authorities when delivering the evacuation command or for practical purposes. For example, it will be easier for people to remember that "people at a distance up to km are in danger" (hazard model b) rather than "people in a medium hazard zone are in danger" (hazard model a). The first hazard model will be used to define the individual risk, while the second hazard model will be used for the decision-making regarding evacuation. This is based on the assumption that using the second hazard model will result in fewer errors compared to using the first hazard model (Jumadi et al. ), while the second hazard model does not directly mean the actual hazard, and so is inappropriate for assessing the true risk. Therefore, we implement a dual hazard in order to obtain a better outcome for the evacuation decision while retaining the precision of the risk. .
The hazard is classified into three zones (Figure ). The hazard level of each zone will dynamically change over the course of the simulation as the volcanic activity changes. The rules governing hazard levels are based on the function of several variables concerning the characteristics of the volcano. The VEI and Volcanic Activity Level (VAL) are provided in Table (Jumadi et al. ). The changing hazard levels within those zones are illustrated in Figure . Finally, the hazard level in the agent's location (based on its coordinates) is classified and scored   ). The attributes of the agent (person) are scored and the score for each attribute is defined using a pairwise comparison weighting based on local knowledge and a literature review (Saaty ). Each theme of the attribute is ranked based on the vulnerability level (Cutter et al. ). The result of the pairwise comparison analysis is shown in Table . In the absence of detailed data indicating otherwise, the classification used by Cutter et al. ( ) is interpreted as a simple binary value of "vulnerable" or "less vulnerable". This gives rise to values in the Pairwise Comparison Index of . and . , respectively. While this may simplify the SoVI, as a demonstration of the utility of this approach this is deemed su icient for the purposes of this paper. Finally, we aggregate the social attributes using Equation to calculate the SoVI (Chakraborty et al. ).
The calculation of the individual risk is based on the risk concept proposed by UNISDR ( ) and Sar et al. ( ), as previously explained in Section . We express the individual risk as a certain quantification that can be classified. The formula for providing the value is presented in Equation with a weighting rule based on Table . This equation generates the risk of individuals as a risk index (R i ) using a weighted linear combination (W LC) (Malczewski ). Once the index has been obtained, it is classified into three classes, based on Table  . R

Implementation, experimentation and analysis .
The model is implemented using AnyLogic, a multimethod (Agent-based, System Dynamics, and Discrete Event) simulation modelling tool developed by the AnyLogic Company (Borschev ). An overview of the agents' statechart to express the behavioural rules of the agents is provided using Unified Modelling Language (UML) in Figure , while detailed documentation of the ABM application is provided at the OpenABM repository (https: //tinyurl.com/stdmr). These statecharts show the implementation of the model in AnyLogic (Jumadi et al. , see) for the conceptual design and OpenABM repository (https://tinyurl.com/stdmr) for the code). A statechart is typically a state transition diagram used to define event-and time-driven behaviour in AnyLogic (Grigoriev ). There are three main statecharts for human-agents (Figure a), including observing the hazard, evacuation decisions, and alerting the community. The 'observing hazard' statechart enables the humanagents to sense the hazard in their location and classify its level, based on Table . The 'evacuation decision' statechart is used by the agents to decide whether or not they need to evacuate (see Section . and Jumadi et al. for details on the evacuation decision model). When the human-agent feels in danger, decisions are made by the individual regarding whether to evacuate or to stay where they are. Meanwhile, the volcano and stakeholder agent have only one statechart each. The volcano is utilised with the statechart of volcanic activity generator, thereby giving the stakeholder agent the ability to make decisions based on the volcanic activity. The stakeholders will alert people when the volcano displays a high level of activity. .
Using the developed model, experiments can be used to show how the spatio-temporal dynamics of risk vary with the magnitude of the hazard (VEI) and the crisis length. Here, we explore the eruptive activity via a scenario using VEI and a Crisis Length of days to mimic the eruption that occurred in (Mei et al. ; Surono et al. ). The results, presented as the coordinates of people and their associated risk level, are then analysed using kernel density analysis to identify the risk hotspots (Lin et al. ; Thakali et al. ). To produce the final risk hotspot, we run the experiment times to provide su icient samples for statistical analysis based on the central limit theorem (Ghasemi & Zahediasl ; Haneberg ) and spatially average the results.

Verification and validation
. In implementing the model structure discussed above, we need to ensure that the model works in regard to the overall concept (verification) as well as ensuring a close fit to the real world (validation). Verification was conducted through visualisation (Hawe et al. ), by means of graphical output monitoring and inspecting the statecharts. For validation, we used a retrodiction approach among various other validation techniques (Hawe et al. ). This approach focuses on measuring the replicative validity of the model -i.e. the ability of the simulation to generate output that matches the real data (Troitzsch ). We compare the spatial pattern of risk in the predicted maximum level of hazard intensity with the real spatial pattern of casualties that occurred in . For this purpose, we used a qualitative approach to provide the comparison of both real data and the simulation outputs (Crooks & Hailegiorgis ) in which the plausibility of the output was judged by the authors.

Results
. Section provided a technical overview of the model; we now describe the simulation experiments and spatiotemporal analysis of the results, and discuss the outcome with reference to related works, highlighting the potential implementation with regard to supporting emergency management.
The purpose of the experiment is to highlight the veracity of the approach (coupling an ABM with a physical hazard model and an MCE to determine individual vulnerability and, consequently, individual risk) and to show that the overall outcomes are valuable for practical emergency management. The outcome of the experiment can be saved for further analysis as well as directly reviewed during the simulation (e.g. Figure ). Figure shows the result, illustrating the spatial distribution of the people at risk as well as the dynamic of the risk level. Subfigures a-d illustrate the changing level of risk as the emergency develops. Initially, most individuals have low (or negligible) risk (see Figure a) but, as the hazard spreads, the risk becomes far greater (see Figure b and c). In Figure d, few people at risk remain due to the mass evacuations when most people move away from the hazard zone during a high level of volcanic activity. The remaining people are considered reluctant to move as a consequence of the variability of individuals' decisions (see Jumadi et al. ). . The saved output from the simulations is used to provide the spatio-temporal densities of those at risk by demonstrating the dynamic responses of the agents to the changing hazard levels. The Kernel density maps provide a better approximation of the spatial distribution of those at risk compared to the distribution of points (Figure ), because the agent population is distributed randomly using the Monte Carlo approach (Section . ). Therefore, using the point distribution of those at risk directly to understand the risk can be misleading. We used GIS (see Sections . -. ) to explore the dynamic risk over time by analysing individuals' locations and attributes (risk level attribute). Figure presents the results of the varying density of the individuals at risk at various time points in the simulation. From these results, we can see that the level of risk to humans can change dynamically depending on the state of the volcano agent. The risk values not only depend on the level of hazard but also the number of people at risk. This model can show the impact of the evacuation procesess on the outcome of disaster risk reduction responses. This result confirms the importance of providing risk assessments through a dynamic rather than static model. Figure : The STDRM Analysis. The densities of people at risk (a) are applied to the people at risk's distribution graph (b), the simulated temporal volcanic activities graph (c) and the temporal curve of the percentage of the evacuees' accumulation (d).

The risk hotspots
.
In this study, we use the term 'risk hotspots' to indicate the geographical locations regarding people at risk who are reluctant to evacuate during the simulated crisis. Hotspots are defined as relatively high-risk locations (Thakali et al. ). To create a hotspot, we analyse the density of the individuals (using kernel density analysis) who are at risk as the volcanic activity increases. A risk hotspot is, therefore, a place with a substantial concentration of people who are at risk and reluctant to leave at a time when the activity of the volcano is high. We captured the distribution of individuals who remain in their location until the end of the volcano's highactivity period (see Figure ). The Risk Hotspots are provided in Figure . From the figure, it is clear that the risk hotspots are mainly located in three areas. The first is in Cangkringan District, the second around Ngaglik District, and the third around Tempel District. These are in the high and medium hazard zones, where individuals are at risk of direct volcanic events, such as toxic gases, nue'es ardentes, and PDC, especially where these hazards are concentrated lower down the volcano slopes through topographic channelling e ects long river valleys. Figure : The risk hotspots.

Validation of results
.
The most striking result of the simulation is that we can now highlight the risk hotspots as an emerging result of the evacuation decisions of individuals during a crisis, that has a similar spatial pattern as the real casualties' distribution data. This can improve the decisions of disaster managers in focusing resources to mobilise and facilitate evacuation processes in hotspot areas. These hotspots can also be used to identify the areas that may attract potentially high casualties when disasters occur. .
A qualitative comparison between the risk hotspots and the distribution of casualties in the eruption reveals close similarities (Figure ). This map was derived from casualty data provided by the local government of Sleman (Pemkab Sleman ). The distribution of casualties map also shows a relatively high percentage of casualties in Cangkringan District, but there are discrepancies in Pakem and Turi Districts.

Discussion and Conclusions
. In this paper, a new approach of Spatio-temporal Dynamics Model of Risk (STDMR) is proposed and a case study using a pre-developed agent-based evacuation model of Mt. Merapi is provided (Jumadi et al. , ). This approach first creates an individual-level population (synthetic population) of agents who live in the area surrounding a volcano. Each agent has a unique vulnerability and, since vulnerability comprises several factors (Cutter et al. ), MCE is used to create a single social vulnerability index for each individual. This is coupled with a dynamic hazard model to capture the dynamics of risk. The model is able to highlight a small number of high-risk spatio-temporal positions where, due to the behaviour of individuals evacuating the volcano and the dynamics of the hazard itself, the overall risk at those times and in those places is extremely high. Nevertheless, the outcomes are interesting and extremely relevant for the stakeholders, and the work of coupling an ABM, MCE, and dynamic volcanic hazard is both novel and relevant for decision-makers and disaster management experts alike. .
The STDMR integrates the MCE-based individual risk model into an ABM simulation, and can show the impact of the evacuation procesess on the disaster risk reduction outcome. This can be used to measure the e ectiveness of evacuation plans in reducing risk. For instance, this can to be used to compare how well a staged evacuation strategy works compared to a simultaneous evacuation strategy (Jumadi et al. ). The staged strategy has been proven to be a better strategy compared to the simultaneous strategy in reducing clearance times by minimising the potential tra ic congestion (Chen ; Chen & Zhan ), whereas the e ectiveness of this strategy, considering various groups of risk, has been demonstrated by Jumadi et al. ( ) using the STDMR approach. .
Moreover, this approach can improve the existing static GIS-based risk analyses that are commonly conducted in hazardous areas/regions (Alcorn et al. ; Martins et al. ). It achieves this in two ways. Firstly, enabling the population at risk to move during the crisis the model creates a considerably more realistic spatial distribution of the population. Secondly, accounting for the individual risk to people as well as the dynamic volcanic activity makes the resulting pattern of risk far more realistic. This model can provide an individual-level of risk that can produce a more detailed spatial pattern of risk compared to regional-level risk analysis. .
The integration of MCE-ABM for STDMR has been successfully presented in this paper to illustrate the dynamic changes in volcanic risk in a spatio-temporal manner. The ability of the model to demonstrate the e ect of the evacuation processes on the risk reduction outcome can be used to measure the e ectiveness of various evacuation plans to reduce the risk. Moreover, from the simulation, we present the risk hotspots that reflect the concentration of those at risk at particular sites as the outcome of the evacuation decision of individuals. This simulation can potentially be used to increase the decision-making processes regarding evacuation. Knowing the hotspots can help the decision-makers to allocate more resources towards managing and mobilising these areas.
. However, there some limitations that can be addressed for future research. These include the weighting method, general applicability, and further improvement of the model rules. For the weighting method, we suggest adding more discrimination to the Social Vulnerability Index (SoVI) beyond a binary classification where detailed data are available. In addition, this model needs further analysis of its robustness by testing how the approach might be applied in the context of either other volcanoes or di erent hazard types, such as flooding, etc. where there is a strong spatio-temporal component. Moreover, this model can be improved in terms of the decision-making rules of the agents and the agent interactions. For example, the destination choice rule should be improved, as this has not yet been calibrated or validated. It is important to compare the distribution of evacuees with the real data as this reflects the validity of the destination choice rule of the agent. In , the populations within the danger zone in Merapi evacuated to temporary shelters (evacuation centres) distributed outside the danger zone. These shelters were commonly public facilities such as stadiums, schools, mosques/churches, etc. Moreover, the probability of unsuccessful contact following the agent interaction might a ect an outcome that is ignored in this model. This model assumed that all agents always successfully make contact with their connections and always follow the commands given. Also, it is possible for people to ignore the evacuation order altogether.