Mount Merapi eruption: Simulating dynamic evacuation and volunteer coordination using agent-based modeling approach

Humanitarian logistics operation is regarded as unregulated system, a sudden surge of demand, high uncertainty with respect to location, timing and intensity, of which then leads to either resource scarcity or oversupply. Challenges in the operation are corresponding to types of disasters, phases of disaster relief and types of humanitarian organization (Kovacs and Spens, 2009). For instance, disease surveillance is of high importance after the flood disaster, whereas evacuation operation after disaster is the most crucial phase to avoid great losses and fatalities in volcanic disaster. When it comes to volcanic disaster, unpredictable behavior of volcanic eruption takes along another challenge. In the 2010 Mount Merapi eruption, the unpredictable pyroclastic behavior which had passed through nearby shelter caused the closure of the shelter and the movement of Internal Displaced Persons (IDPs) away to a new safe shelter, eventually resulted in the evacuation dynamics. Furthermore, different evacuation behavior of the IDPs complicated the evacuation operation as some IDPs were reluctant to evacuate due to low hazard perception and strong trust to cultural leader who refused to evacuate (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018). Similar evidence was also observed in volcanic disaster of Mount Karthala (Morin and Lavigne, 2009). Evacuation decision making and cultural aspect had made the evacuation operation hard to envisage, and thus making even more difficult to plan for the aid distribution as the number and the location of the IDPs were volatile.

The emergency response during the evacuation, i.e. last-mile relief distribution, also required extra handling due to mobilization and a large supply of resources in such a short time. Instead, the goal of the last-mile relief distribution is the short delivery time and the optimal fulfillment of demand. Subsequently, the role of volunteers to channel aid is of importance. The movement of volunteers to provide humanitarian logistics assistance is not without obstacles. The involvement of numerous actors, such as non-governmental organizations (NGOs), which generally worked independently and may compete for the resources, confronted efficient aid distribution operation. Therefore, coordination is required to increase efficiency and responsiveness of the logistics operation (Jahre et al., 2009). Although coordination has been well-understood and successfully implemented in a commercial context, unfortunately, only a few success stories can be found in humanitarian context. Collaboration seems to be one of the fundamental weaknesses in humanitarian operations (Kaynak and Tuğer, 2014; Rey, 2001). The coordination-related problem in humanitarian context is associated with inefficient (waste valuable resources) such as uneven distribution – a surplus of personnel/supplies in one area, and a dearth in another, and ineffectiveness (mistarget beneficiaries) such as disparities between capacities and a number of needs – inappropriate personnel and aids. Without coordination through information sharing, invalid and incomplete information make aid operations less optimal – inefficient information sharing might cause relief supplies do not match the needs of survivors or relief supplies being sent are blocked due to lack of last-mile distribution arrangement. Kaynak and Tuğer (2014) have investigated that coordination appeared to be the main problem in volunteer movement. The coordination is therefore required to manage different players in aiding people in need to reduce redundant effort and provide better service (Moore et al., 2003), thus becoming an integral aspect of emergency response efficiency.

The present study addresses two research objectives; first, understanding the individuals’ evacuation decision making underlying the evacuation dynamics; and second, exploring appropriate coordination mechanism and other potential strategies in favor of efficient and effective last-mile relief distribution. The aim of the study is therefore to develop an empirically founded agent-based simulation model of the evacuation dynamics which is to be used as an experimental tool. An agent-based model was used as a modeling platform to take into account individuals’ evacuation decision making, dynamics aspect and interactions (i.e. between pyroclastic movement and IDPs, between IDPs and volunteers, among IDPs and among volunteers). The 2010 eruption of Mount Merapi in Sleman district is taken as a studied case because the unexpected pyroclastic movement had made it a suitable case to represent the evacuation dynamics, which provides appropriate setting to evaluate the effectiveness of coordination and other operational strategies. The 2010 eruption of Mount Merapi, the most active volcano in the world, was a historic natural phenomenon for Indonesia because the eruption was the largest in the last hundred years. It was recorded that there were 61,154 IDPs, 341 of casualties – the highest recorded casualties and 5 trillion Rupiah of the losses suffered (Mei et al., 2013). Small eruptions have still occurred in every 4–6 years.

The characteristics of humanitarian logistics are very much different from those of commercial logistics. Uncertainty, dynamics and adaptability are attributable to humanitarian logistics (Rey, 2001). Agent-based modeling (ABM) approach, which is commonly used to model complex system, was therefore selected because the ABM has the capability to capture heterogeneity of individuals with respect to attributes and behavior, to transform the behavior of every individual into a simpler logic, to model interactions (between individuals, between individuals and environment) and to facilitate a more representative spatial environment of the simulation setting by embedding Geographical Information System (GIS) data in the agent-based simulation model (Wilensky and Rand, 2015).

Coordination in humanitarian logistics has been gained substantial consideration in the literature due to its critical role to improve efficiency and effectiveness of response effort. The fact that a single organization is unlikely capable of handling relief operations has made coordination a necessity. The large number and diverse actors involved in relief operations have raised challenges in coordination in order to support relief mission while fairly distributing risk and benefit (Balcik et al., 2010).

The literature focusing on coordination in humanitarian logistics has varied from qualitative studies (e.g. theoretical framework and empirical studies) and quantitative studies (e.g. optimization and simulation). The qualitative study of coordination literature attempt to explore coordination mechanism from organizational perspective such as network centrality Moore et al. (2003), warehousing collaboration (Balcik and Beamon (2008), cluster approach (Jahre and Jansen, 2010), centralized vs decentralized coordination (Dolinskaya et al., 2011), the concept of chain coordinators (Akhtar et al., 2012), transitional coordination and coordination structure (Clarke, 2013), The International Search and Rescue Group and Urban Search and Rescue (Tatham and Spens, 2016), to name a few. Holguin-Veras et al. (2012) points out that the coordination is addressed as a socio-technical process whereby a social network of actors orchestrates a set of technical activities. It implies that understanding the functioning of the entire system requires proper consideration of its elements. Therefore, some literature (for instance Heaslip, 2012; Kabra et al., 2015) conducted empirical studies to identify barriers on coordination which highlights that cultural and people barriers exist in addition to managerial, technological and management barriers. While all of the aforementioned literature are useful in understanding the different mechanism of coordination, the literature do not attempt to model the actual system and to quantitatively evaluate the effectiveness of the collaboration mechanism. Only few literature have focused on modeling humanitarian operation quantitatively. Due to the characteristics of humanitarian operations which are stochastic, dynamics and adaptive, the simulation approach is preferable than analytical models (Holguin-Veras et al., 2012). The ABM, as a simulation approach, has gained more attention to be used in humanitarian logistics. For instance, Comfort et al. (2004) built an abstract agent-based model to simulate the role of information to support coordination, Das and Hanaoka (2014) modeled resource allocation at upstream operations, Mochizuki et al. (2015) addressed sourcing options in horizontal cooperation, Suarez-Moreno et al. (2016) modeled coordination between donors and humanitarian staff to understand the effect of priority donations, the number of warehouses and stock-out probability toward material convergence, and Aros and Gibbons (2018) simulated communication media options among inter-organizations. In addition to modeling the upstream operations, the agent-based model was also used to model downstream operations, namely, evacuation; for instance, earthquake/tsunami evacuation (Mas et al., 2012), building evacuation in case of fire (Tan et al., 2015), Haiti earthquake (Crooks and Wise, 2013).

Henceforth, there exist gaps in the aforesaid literatures. First, many existing agent-based models were however set in a hypothetical environment, only very few were built to represent a particular system, i.e. Mas et al. (2012). Second, most of the simulation efforts have been focused on either upstream operations (i.e. sourcing) or downstream operations (i.e. evacuation), so that an integrated model which addresses both evacuation process and demand fulfillment operation, particularly during last-mile relief distribution is still lacking. On the other hand, the last-mile relief distribution is uniquely challenging in relief operations (Balcik et al., 2010). Third, quantitative evaluation of the effectiveness of the coordination mechanism is not yet available in the existing literature, despite that many qualitative literature studies have proposed various potential coordination mechanisms.

To fill the gaps, the present study, therefore, contributes in three ways. First, the present study considers human and cultural factors which were acquired from empirical studies when modeling the evacuation. GIS data are employed to facilitate real environment setting in the simulation model. Hence, from the methodological point of view, the present study demonstrates the combined approach of empirical studies and simulation within the context of humanitarian logistics. As the simulation model is developed to mimic the real system, it is expected that accurate prediction can be attained. Second, the developed model integrates the evacuation process which determines aid demands (downstream) and coordination mechanism which influences aid supplies (upstream) during last-mile relief distribution in humanitarian context. Third, the present study complements existing literature by providing quantitative evaluation of various strategies proposed by previous qualitative literature.

The paper is structured as the following: the Section 1 has highlighted the background, the research gaps and the contributions of the present study with respect to methodological and empirical contributions. Section 2 presents literature review on evacuation decision making and coordination. Section 3 describes the methodological approach of the empirically founded agent-based model, followed by the agent-based model including verification, calibration and validation tests, sensitivity analysis and scenario development in Section 4. Section 5 presents results and discussion, which is then followed by managerial and policy implications in Section 6, before presenting the conclusion in Section 7.

This section reviews literature on evacuation decision making and coordination mechanism. The evacuation decision making emerges the evacuation dynamics, and eventually the demand for aids (downstream). Coordination mechanism describes the processes of demand fulfillment by relief actors such as NGOs and government (upstream). The section delivers sufficient foundations to develop a conceptual model of the agent-based simulation model and to provide a definition and scope of coordination within a humanitarian context.

Following Vorst (2010), it is necessary to develop a simulation model which addresses human behavior in evacuation model because it provides a realistic prediction and, thus better problem solving. It was also evidenced that in the 2010 Mount Merapi eruption, the evacuation was not easy because many people refused to evacuate due to low perception of hazard, belief, social engagement with cultural leader and a strong sense of belonging toward their valuables and livestock (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018). Although the local government has provided shelters and evacuation route, it was recorded that there were 341 of casualties, which was the highest recorded casualties. Two factors are identified as the major causes of the inefficient existing contingency and evacuation plan, i.e. the behavior of people during Merapi disaster (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018) and the uncertainty of eruption behavior of Mount Merapi (Mei and Lavigne, 2013). Consequently, it is, therefore, a necessity to design an evacuation plan which takes into account individuals’ evacuation behavior and addresses the dynamic nature of the eruption.

A literature review has been conducted to explore underlying factors of evacuation behavior as well as types of evacuation behavior. It was found out that individual people evacuate at a different time. According to Perry (1979), several factors determining people decision to evacuate, on a general level, are adaptive planning, threat perception, hazard level, family, kin relationships, community involvement, age and culture. Further, Perry (1979) mentioned that the higher the individual’s perception of the threat, the higher the probability for evacuation. Furthermore, the more often a person experiencing a disaster, he/she will be more confident about the threat that occurred. As more and more warnings are perceived, more threats will be felt. Warning letters from trusted sources increase the degree of confidence that the threat is true. Perceptions of threats depend on the extent of the danger. The higher the level of disaster, the higher the people desire for evacuation. Furthermore, the more dangerous warning messages will be given to the location of the disaster impacts, the perceived level of personal hazards will increase. Previous studies specifically on the Mount Merapi eruption (e.g. Lavigne et al., 2008; Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018) have indicated that disaster experience, disaster training, trust on the cultural/local leader, vulnerability, ownership of livestock/valuables are all influencing the evacuation decision making. The findings are also supported by Rahman (2017) who statistically analyzed the significant relationship between hazard perception and time to evacuate, which were acquired through a survey toward population residing surrounding Mount Merapi. Once a person decides to evacuate, the next decision is then how to evacuate, referred to as evacuation behavior.

When it comes to the evacuation behavior, a number of researches have studied evacuation behavior in the face of natural disaster, for instance, Lavigne et al. (2008) for volcanic hazards and Pan (2006) for earthquake disaster. Sopha, Asih, Ilmia and Yuniarto (2018) and Sopha, Asih, Nurdiansyah and Maulida (2018) have explored that there are four groups of evacuation behaviors; non-adaptive (incapable to self-evacuate, dependent of others), adaptive (capable to self-evacuation, so-called adaptive behavior because they can change their behavior depending on the situation, for instance, when searching for the route), altruistic (capable to self-evacuation and willing to help others) and leader-following. Adaptive and altruistic behaviors use hazard perception as a trigger to start evacuating, whereas non-adaptive behavior depends on the decision of neighboring agents performing adaptive behavior, and leader-following behavior depends on the evacuation decision of the cultural leader. During the eruption in 2006, the cultural leader faithfully waited for Mount Merapi eruption until it was completed and survived from the danger of hot clouds/pyroclastics. Courage figure of the cultural leader in 2006 had become a reference for some residents not to evacuate in the 2010 eruption. The cultural leader was believed as a guard of Mount Merapi.

Following the aforementioned empirical evidence, the present study models two stages of evacuation decision making, namely, “when to evacuate” and “how to evacuate” decisions. The “when to evacuate” decision is influenced by hazard perception which is heterogeneous among individuals. The present study adopts that “how to evacuate” decision consists of three types of evacuation behavior, i.e. adaptive behavior (evacuating to the nearest aid shelter), non-adaptive behavior (following the behavior of nearby individuals with adaptive behavior) and cultural leader-following behavior (following the behavior of the cultural leader). Individuals with cultural leader-following will evacuate according to the evacuation decision of the cultural leader – if the cultural leader decides to stay, the people will stay regardless of obvious signs of hazard. Individuals with adaptive behavior evacuate to the nearest shelter through the shortest path.

Coordination is very common in commercial supply chains, however, coordination in humanitarian context is a challenging task due to its complexity when it comes to uncertainty, rapid change in humanitarian environment (difficult to get accurate and real-time information), as well as, large number and diversity of involved organizations (Sarcevic et al., 2012). Literature, i.e. Rey (2001), Moore et al. (2003), Jahre et al. (2009), Balcik et al. (2010) and Kaynak and Tuğer (2014), have demonstrated that coordination plays an important role in facilitating efficient and effective logistics in the humanitarian context. Coordination is a must in humanitarian operations because a single organization is incapable of fulfilling the needs of affected people without working together with other organizations (Akhtar et al., 2012). Moreover, different and multiple responding entities, being it national, international, governmental, non-governmental, require better inter-agency coordination mechanism (Tatham and Spens, 2016). Due to fundamental differences between commercial setting and humanitarian setting, coordination mechanism working for commercial setting may not be appropriate in a humanitarian setting.

Humanitarian logistics literature has defined coordination in the humanitarian context. Coordination is defined as the relationships, interaction/interdependency among different actors in the humanitarian environment (Balcik et al., 2010). The environment involves international organizations, host governments, aid organizations, the military, private companies and others that may have different interests, capacities and logistical expertise. There is usually no single actor who has enough resources to respond effectively to major disasters. Further, Balcik et al. (2010) explained the interaction may refer to information sharing, resource sharing, joint decision making, joint projects and regional division of tasks, with various strength of relationships among players involved. Another definition by Ozlem et al. (2014) mentioned that coordination is the alignment of humanitarian actors in providing aids in the most possible effective and efficient way in which spatial equity distribution and resource efficiency are the two important performance indicators in humanitarian operations. The alignment could be in the form of resource sharing of tangible and intangible assets (processes, people, management skills, experiences and benefit sharing) (Tatham and Spens, 2016). Among the types of coordination, i.e. information sharing, resource sharing, financial sharing, information sharing seems to be the most important and time-sensitive type of coordination, because information sharing has multiple dimensions with respect to timeliness (collected and distributed), validity and relevancy to the field. The present study, therefore, applies information sharing as a coordination type in the simulation model.

According to the type of engagement, coordination can be categorized to strategic level (long term), tactical level (medium term) and operational level (short term) (Balcik et al., 2010). Furthermore, coordination can be classified into vertical coordination – the extent to which an organization coordinates with those performing different operations/roles in serving relatively similar end customers, and horizontal coordination – the extent to which an organization coordinates with other organization performing similar operations at the same level (Balcik et al., 2010). In disaster response situations, coordination should be carried out by all stakeholders such as government, private sectors, NGOs and communities.

Coordination mechanism, also known as coordination structure, refers to a set of methods to manage interdependence between the organization (Xu and Beamon, 2006). Dolinskaya et al. (2011) categorized humanitarian coordination mechanisms in humanitarian operations into two mechanisms, i.e. centralization and decentralization.

Centralization is the process of coordination with a single agent that has the authority to direct coordination activities. The centralized system, also known as top-down style, is characterized by the presence of a central player. The central player functioned as a director of the relief efforts. The player also gathers information and makes the decision for collaborative participants. Coordination is always done with a command approach in which the central agent controls all logistical resources and makes decisions for the agencies involved. In centralized coordination systems, often the central government plays the role of central actor and is responsible for logistical coordination decisions. The advantages are easier responsibility placement, centralized data, control, reduced duplication, improved security, rapid system development, information flow, reduced cost and fast movement of resources and information. The associated disadvantage is that the time required to implement the decision becomes longer. This is because the decision of the center must go through several stages to get to the field so it takes more time and makes the movement of aid slower. An example of a centralized system is relief operation in 2000 Mozambique flood in which the World Food Programme (WFP) and UN High Commissioner of Refugees centrally arranged the usage of available vehicles and the schedule of supply deliveries (Sami and Van Wassenhove, 2003). In 1991, the United Nation adopted the centralized system by establishing an organization focusing on facilitating coordination, namely, the Office for the Coordination of Humanitarian Affairs (OCHA, 2014), which conducts need assessment, plans operations, shares information and divides tasks among organizations.

Decentralization involved multiple organizations in which logistics coordination is regulated by agreement. A decentralized system is characterized by which each relief organization makes its own decisions. Most of the disaster response is inherently decentralized. There is no single organization having full authority over the person or organization in humanitarian operation. Each humanitarian aid organization makes its own decisions to share information, disseminate and be accountable to other volunteers. The decentralized system is usually implemented in the case of diverse organizations. The advantage of such a system is faster decision making as each region is given the opportunity for decision making. Disadvantages are related to increased administrative costs, uneven distribution, less uniformity and consistency, the absence of control or definite procedures and less optimal decision. An example of a decentralized system is medical work in disaster response during the Haiti earthquake 2010. When every player has the same objective and different abilities, the decentralized system may cause potential competition of players for scare resource (Kent, 2004).

There is currently no agreement on which system, either centralized or decentralized, is the best approach. Some studies, e.g. Kehler (2004) documented the success and shortcomings of the UN centralized system. Wex et al. (2011) stated that deficiencies of past emergency response are due to the lack of centralized coordination to control over decentralized actions, whereas Stephenson (2006) argued that the decentralized system is preferred to facilitate coordination, particularly when the magnitude of the disaster is overwhelming. Further, instead of top-down control, relief actors should be reconceived as social networks. The decentralized system is also supported by Sarcevic et al. (2012). Kehler (2004) summarized that successful and effective coordination depends on capacity and contextual conditions.

Coordination occurs at different logistics processes, i.e. procurement (supplier–buyer alliances, collaborative procurement), warehousing (standardization of methods, warehousing through a third party) and transportation (shipper collaboration, the use of 4PL). Some simulation studies in relief operations have addressed collaborative procurement (e.g. Das and Hanaoka, 2014), some other quantitative studies have evaluated warehousing collaboration (e.g. Mochizuki et al., 2015; Suarez-Moreno et al., 2016). Given that a simulation study on transportation collaboration in relief operation is, according to the best knowledge of the authors, still lacking, the present study focuses on collaboration during last-mile relief distribution. Moreover, coordination in the last-mile relief distribution was not implemented in the 2010 Mount Merapi eruption so that it is worthwhile to evaluate different coordination mechanism in order to increase the efficiency and effectiveness of future relief operations.

ABM approach has been widely applied in a number of logistics/supply chain domains such as agricultural supply chain (e.g. Utomo et al., 2018) and city logistics (e.g. Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018). However, the application on humanitarian logistics is still very few. ABM is particularly useful in humanitarian logistics because the humanitarian logistics is characterized by high uncertainty, heterogeneous actors/agents, adaptive agents, interactions (agent-agent and agent-environment), geographic and location-dependence, and dynamic behavior of the system. Those characteristics are best to be modeled by ABM than other modeling techniques (Wilensky and Rand, 2015).

ABM was, therefore, used as a modeling platform to implement the heterogeneous individuals’ evacuation behavior and the interactions among individuals and between individuals and environment, which generates/emerges the evacuation dynamics. The simulation model addresses the evacuation decision making by modeling the movement of IDPs from their actual home location to the nearest shelters. Constraints such as transportation capacity, limited mileage at one-go, limited shelter capacity, limited road capacity are considered. Once the number of the evacuated IDPs exceeds the road capacity, bottleneck occurs and consequently slows down the IDPs in reaching the shelters. Volunteers and IDPs meet at the shelters where demands of the IDPs are fulfilled by the supplies from the volunteers. Limited shelter capacity causes the dynamic movement of the IDPs from one shelter to another shelter. The volunteer movement is started from the aid centers to the shelters. Due to limited transportation capacity and limited mileage at one-go, the volunteers should go back forward from the aid centers – the shelters – the aid centers. Hence, the simulation model well-represents lead time in logistical activities. The volunteers move based on the information. The mechanism of information sharing is developed through scenarios. Experiments are conducted to explore coordination mechanism and other potential strategies favoring efficient and effectiveness in last-mile relief operations.

Different from abstract agent-based models which do not represent any specific system in the real world, the empirical agent-based model should be constructed in such as a way to adequately represent the real system. The procedure by Boero and Squazzoni (2005) is, therefore, adopted. The development of the empirically founded agent-based model consists of three stages, i.e. model specification, model calibration and model validation. Model specification was conducted through field study and secondary records analysis to qualitatively select empirically grounded agents’ attributes and behaviors during the 2010 Mt. Merapi disaster (e.g. Rahman, 2017; Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018; Sopha and Asih, 2018). Model calibration aiming at selecting the values of model parameters was conducted by direct calibration (Windrum et al., 2007) through tailored empirical study which was conducted by interviews in 2014 to explore underlying factors and types of evacuation behavior (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018), and by a survey in 2017 (Rahman, 2017) to collect quantitative data based on the 2014 study. Model validation, further discussed in Section 4.2.1, was carried out by comparing the simulation results against government historical records and another independent data by Mei et al. (2013). Validation tests were conducted for both overall model and sub-model (i.e. pyroclastic movement), as well as, at both macro-level (i.e. the number of IDPs over time and the number of casualties) and micro-level (i.e. the number of IDPs by the shelters).

This section presents the development of a simulation model following ABM approach. The purpose of the model are twofold, first, to understand the mechanism of individual behaviors underlying the evacuation dynamics, and second, to evaluate the effectiveness of coordination mechanism and other strategies through what-if scenarios. The model has focused on last-mile relief distribution in which the volunteers deliver the aids to the IDPs as final recipient of the chain at the shelters.

This sub-section describes the implementation of the evacuation decision making and coordination mechanism discussed in Section 2. The theories were translated into an agent-based simulation model, referred to as a model specification. Following Wilensky and Rand (2015), an agent is defined as an autonomous entity which has heterogeneous attributes and behavioral rules, whereas the environment is the landscape on which agents interact. Agents can interact with other agents or with the environment. The interactions influence decision making and agents can update their internal state. Due to the interactions, a system exhibits a property at the macro-level that is not defined at the individual level, so-called emergent properties (Wilensky and Rand, 2015).

The agent-based model has three different types of agents, i.e. IDPs referring to individuals of the affected population, volunteers referring to aid organizations which deliver the aid supplies for the IDPs, shelters where the demands is fulfilled by the supplies and pyroclastics which determines when the individual IDPs evacuate. The individual agents (the IDPs), the volunteers and the pyroclastics are movable agents, while the shelters are static agents.

Based on Sopha, Asih, Ilmia and Yuniarto (2018) and Sopha, Asih, Nurdiansyah and Maulida (2018), individual agents are grouped into children, elders, disabilities (vulnerable people) and adults. These groups have different needs so that the group (consisting of children, elders, disabilities and adults) is one of attribute variables of individual agents. Other attribute variables are location on which the individual agents originally resides, state which reflects to the situation of the agents (stay, evacuate, or settled), demand which reflects the need for food, water, hygiene kits, medical equipment, medications (demand in this present study is approximated using WHO calorie requirement (World Health Organization, 2003)), and hazard threshold which is defined as the perceived hazard level which triggers when-to-evacuate decision. Previous studies on Mount Merapi disaster (e.g. Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018; Rahman, 2017) have indicated that people have different hazard threshold. Based on the 2017 survey, 20 percent of individual agents will evacuate when they perceive that the disaster is dangerous and 80 percent of individual agents will evacuate when they perceive that the disaster is very dangerous.

Volunteer agents perform last-mile relief distribution, delivering the supplies from the aid centers located at both Universitas Gadjah Mada (UGM) and Maguwoharjo stadium, to the shelters. The attribute variables of the volunteer agents are aid capacity which is defined as the number of on-hand aid supplies, and energy which indicates the available energy for delivery. Each volunteer uses pick-up truck with a certain mileage limit due to limited fuel-tank capacity. As the individual agents and volunteer agents meet in the shelters where the demands of the individual agents are fulfilled by the volunteer agents, hence location, unmet demand, and state (i.e. closed or opened) are then specified as the attributes of aid shelters. Table I describes the agents’ attributes in detail.

Using GIS data, the environment of the model is the territory of Cangkringan and Pakem villages, Sleman district, which was affected by the volcano eruption. The GIS includes the actual location of shelters, road networks, aid centers, as shown in Figure 1. With respect to spatial scale, one patch in the simulation model represents 100 meters in the actual system. With respect to the temporal scale, one time-step in the simulation model is corresponding to 2 hours in real system.

The evacuation behavior of the IDPs basically consists of two decision-making processes; “when to evacuate” and “how to evacuate,” as presented in Figure 2.

Following Sopha, Asih, Ilmia and Yuniarto (2018) and Sopha, Asih, Nurdiansyah and Maulida (2018), the evacuation decision is initiated by either a disaster warning or natural phenomenon. “When to evacuate” decision in the simulation model is implemented using pyroclastic movement as a trigger to evacuate. The pyroclastic movement was modeled using seismic intensity (see Table III) to generate affected area in the simulation model (represented as a red patch in the simulation model). The movement can thus be observed by individual IDPs as an incoming threat which in turn influences the individuals’ decision to evacuate. The individual IDPs decide to evacuate when the distance with pyroclastic is less than 15 patches in the simulation model (corresponding to 1,500 meter in the real system). The validation test of pyroclastic movement (sub-model validation) is detailed in Section 4.1.2.

Likewise, the affected area in the simulation model determines the hazard level (disaster intensity). However, the disaster intensity can be perceived differently by the IDPs (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018; Rahman, 2017). The higher the perceived hazard, the more probable the IDPs will evacuate. The simulation model, therefore, uses hazard threshold which is heterogeneous among the IDPs, to represent the perceived hazard level (see Table I). When the hazard level exceeds the individual hazard threshold, the IDPs decide to evacuate. Each individual IDP evacuates at a different time, depending on its individual hazard threshold. The IDPs which are still far from the pyroclastics may decide to evacuate when the hazard level exceeds the hazard threshold. Probabilistic techniques based on by Woo (2008) and Doyle et al. (2014) and empirical distribution of hazard threshold by Rahman (2017) were used to implement when to evacuate decision in the simulation model.

When it comes to “how to evacuate” decision, the decision making of an individual depends on the group in which he/she belongs. Leader-following behavior is simply following the cultural leader. In the present study, the cultural leader decided to stay regardless of the increased intensity of the disaster. The behavior was performed due to a strong attachment to the cultural leader. During the 2010 Mount Merapi eruption, the leader-followers accounted for 0.17 percent of the population, corresponding to 19 people. As the 2010 eruption was more intense and unpredictable than that in 2006, fatalities of the cultural leader and his followers occurred.

Vulnerable groups (children, elders and disabilities) who are unable to self-evacuate follow nearby adults with adaptive behavior. The adaptive behavior refers to self-evacuate behavior, performing high rational thinking when it comes to when to evacuate decision and high behavioral control when it comes to the ability to identify the nearest shelter and move toward the destined shelter (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018). With respect to “when to evacuate” decision, the adaptive behavior employs hazard threshold which was compared against the hazard/disaster level.

In the 2006 and 2010 eruptions, due to the narrow road, the population was displaced on foot. Moreover, road density has also influenced the dynamics of the evacuation process. The simulation model implemented road capacity variable to model the carrying capacity of roads, meaning that the IDPs or volunteers cannot use the road if the road has already occupied and reached its carrying capacity.

Volunteers come from various places, but generally, they pick the supplies from the aid centers at either UGM or Maguwoharjo Stadium. To accommodate the movement of volunteers, 96 pick-up trucks are available for relief distribution to the eight shelters. Each volunteer’s pick-up truck can accommodate relief items to a certain capacity. Each pick-up has a maximum capacity of aid supplies that it can carry in one-go. When volunteers arrive at the shelter, they deliver the aid supplies according to the demand in the shelter. If the capacity is exhausted, the volunteers return to the aid centers for reloading. Volunteers should also return to the aid centers if they run out of fuel. The pick-up truck has certain mileage limits, i.e. maximum mileage until it has to refuel. Similarly, if the fuel is exhausted, volunteers should look for a fueling station. The process of reloading and refueling requires a certain amount of time of one day in the simulation model.

Shelters are the places where the demands generated by the IDPs and the supplies carried by the volunteers meet. Each shelter updates demand for newcomer IDPs and settled IDPs following the required daily calorie requirement by WHO (see Table I). The total unmet demand (D_t) is updated over time to calculate the summation of the unmet demand of new arrived IDPs and the daily unmet demand of the settled IDPs following Equation (1) as below:

Once the amount of total unmet demand and the received supplies are calculated, the balance is evaluated. If the total unmet demand (D_t) is greater than the received supplies (C_t), the shelter updates the unmet demand (D_t+1) and the volunteer updates the remaining supplies/capacity (C_t+1) as the following:

Otherwise, if the unmet demand (D_t) is lesser than the received supplies (C_t), the shelter updates the unmet demand (D_t+1) and the volunteer updates the remaining supplies/capacity (C_t+1), respectively, as the following:

Interactions among agents (i.e. between the IDPs and the shelters, between the volunteers and the shelters) and interaction between the agents and environment (i.e. the IDPs and pyroclastics) occur in the simulation model. The change at one agent affects the decision of other agents, and thus generates unpredictable behavior of the system (emergent properties). For instance, on November 3, 2010, the pyroclastic flow with the temperature of 600°C can move up to 60 km/h and destroy whatever it passes. The uncertain and dynamic pyroclastic movement thus governs uncertainty in the evacuation decision (i.e. when and where to evacuate), as well as, which shelters are still functioned. Due to the pyroclastic movement and increased intensity of the eruption, the SD Jiwan shelter was closed and the IDPs were moved to a safer nearby shelter, i.e. Wukirsari shelter. The number of shelters is reduced to seven shelters and there were increased IDPs in Wukirsari shelter, and so are the demands. Due to the limited capacity of Wukirsari shelter, many IDPs should move to another safe shelter. Furthermore, the closed shelter due to pyroclastics, therefore, creates more evacuation flow, road blockings, which in turn hinders efficient and effective evacuation.

Figure 2 shows the evacuation behavior of agents and interaction among agents to be implemented in the simulation model. The emergent properties to be measured are the unmet demand over time, which indicate the effectiveness of last-mile relief distribution.

Parameterization deals with the selection of the values of model parameters. In addition to two empirical studies (i.e. Rahman, 2017; Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018), other previous studies such as Lavigne et al. (2008), Mei et al. (2013), Budi-Santoso (2014) and secondary records of the number of volunteers from Badan Penanggulangan Bencana Daerah BPBD (Regional Disaster Management Agency managed by the government) were used for parameterization. Table II indicates the parameterization techniques used for the simulation model.

Based on historical data of about 960 volunteers were participated to help the affected population. In total, 96 volunteers are thus generated in the simulation model. Other supporting data were also collected such as WHO standard on daily calorie requirement during disaster which specifies a range of 1,500 kcal up to 3,000 kcal (World Health Organization, 2003). Maximum capacity carried by volunteers was estimated based on the vehicle capacity, which corresponds to 358,000,200 kcal in the real system or 358,000 kcal in the simulation model. The maximum energy available for one-go distribution was estimated from the capacity of the fuel tank which corresponds to 656.6 km in the real system and 6,566 patches in the simulation model (Note: 1 patch represents 100 meters in the real system). Reloading time – time required to return to aid centers and reload the capacity and fuel – was set to be one day.

During the disaster period, the seismic intensity increases and the increase varies over time. The seismic intensity during Mount Merapi eruption was measured daily. Based on seismic data over time (Budi-Santoso, 2014), the pyroclastic movement was estimated and implemented as the affected area in the simulation model as shown in Table III.

During the 2010 eruption, pyroclastics was presence from Gendol river (northeast of Cangkringan sub-district) to the southeast covering Opak river and Kuning river, as indicated in Figure 3 in which pyroclastics move from northeast to southeast. The non-significant difference (t(18)=12.500, p=0.822) between seismic intensity data and generated affected area in the simulation model indicates that the generated pyroclastic flow in the simulation model is able to represent the observed pyroclastic movement in the real system.

Summarizing, all parameter values in the model were acquired from empirical data so that none of the parameter values was based on hypothetical data. It can, therefore, be argued that the model has been directly calibrated (Windrum et al., 2007).

The agent-based model specified in Section 4.1 was implemented in Netlogo (Wilensky, 1999). The spatial data were prepared using ArcGIS software and then exported to Netlogo. The first eruption occurred on October 26, 2010. The government issued a warning and gave an evacuation order on that day. However, many individuals had not yet been evacuated. Only 20 percent of the population had decided to evacuate on that date (Rahman, 2017). The biggest eruption which exploded volcanic ash as high as 17 km and released pyroclastic up to 16 km from the center of the explosion occurred on November 5, 2010. The evacuation order was again issued. Based on BPBD records, 11,259 people including 3,720 vulnerable groups who reside in 23 hazard-prone sub-villages were affected. Eight shelters based on government evacuation plan were put in operation (Table IV).

The simulation was initialized by setting the environment (road networks and pyroclastics), the shelters, the volunteer agents and setting the individual agents on their actual geographical location and assigned their attributes specified in Table I. As one individual agent in the simulation study represents 10 individuals in the real system, the simulation model, therefore, generates 1,130 individual agents and 96 volunteer agents. Table V presents the distribution of individuals based on residential locations (see Figure 1(a) for spatial distribution during initialization).

As one time-step in the simulation is equivalent to two hours in the real system, one day is corresponding to 12 time-steps (ticks). The simulation which models post-disaster operations is run for 20 days (October 26, 2010, until November 14, 2010) which is corresponding to 240 ticks. In each time-step, each agent follows the decision making specified in Figure 2. By the end of each time-step, the number of IDPs at each aid shelter and the total unmet demand were recorded.

Due to the presence of randomness in the simulation model such as when generating hazard threshold, the simulation model produces different results in each simulation run. Therefore, replications of 50 simulation run with identical initial settings were conducted to obtain representative simulation results. Simulation results are hence presented using the mean value of these replications.

Verification which ensures that the model is implemented as the intended design was conducted through structured debugging walk through following Macal and North (2007). Moreover, extreme condition test was also executed (for example by setting the individual agents having no energy, and the result should indicate that none of the individuals reach the shelters). Validation examining to what extent a simulation model is representing the real system was also conducted. In ABM, validation comprises of conceptual model validation, data validation and output validation (Macal and North, 2007). To test the extent to which the model is able to resemble the actual system, three validation tests are therefore employed. With respect to conceptual model validation, the conceptual model was empirically founded and validated toward the experts (Sopha, Asih, Ilmia and Yuniarto, 2018; Sopha, Asih, Nurdiansyah and Maulida, 2018) and a sampled population. Second, the model was parameterized using the empirical study (Rahman, 2017) which was specifically designed for model parameterization, and no hypothetical data were used in the simulation model (see Section 4.1.5 for details). Third, output validation was conducted by contrasting the number of IDPs resulted from the simulation to that based on the historical record of BPBD – the Regional Disaster Management Agency – (as shown in Figure 4) and another independent data of Mei et al. (2013) as shown in Figure 5.

Statistical analysis was conducted to quantitatively test whether or not the distribution of simulated data on the number of IDPs at each shelter and the number of casualties is comparable to that of empirical data. It is found that, based on χ² test, there is no statistically significant difference (χ²=11.631, df=8, p=0.168) between the accumulated IDPs arrived at each shelter and the estimated casualties generated by the simulation and those based on BPBD records (after being down-scaled by 10 as one simulated agent represents to 10 people in the real system), indicating that the simulation model is able to reproduce a stylish fact of the number of IDPs arrived at each shelter and the number of casualties.

In order to increase model confidence in representing the real system, the total generated number of IDPs is contrasted with another independent data provided by Mei et al. (2013) as shown in Figure 5. Similar result of statistically non-significant difference between the normalized IDPs pattern of Mei et al. (2013) and the normalized simulated data is also evidenced (t(38)=−0.043, p=0.966). It implied that the simulation model exhibits reasonable representation for the actual system, and therefore can be used to explore potential strategies toward efficient and effective last-mile relief distribution.

It is worth mentioning that the conceptual model of the evacuation process in the simulation model is very much influenced by the “when to evacuate” decision and “how to evacuate” decision. Therefore, it is necessary to analyze the simulation results in the presence of changes on key variables. Sensitivity analysis was conducted to evaluate the robustness of the simulation model with respect to “when to evacuate” and “how to evacuate” decisions. Furthermore, this section provides insights on managerial standpoints toward effective evacuation process.

With respect to “when to evacuate” decision, two simulation settings, namely, immediate evacuation and late evacuation are analyzed. The base model employs 20 percent of the population evacuate when they perceived that the disaster is dangerous, whereas the rest of the population evacuate when their perception toward the disaster is very dangerous. Vorst (2010) has indicated that although the signs of disaster were recognizable, a significant proportion of the population refused to evacuate. It is therefore interesting to assess if 100 percent of population evacuate when the disaster perceived to dangerous (referred as immediate evacuation), and on the contrary, 100 percent of population evacuate when the disaster is perceived to be very dangerous (referred as late evacuation).

Figure 6 shows that the earlier the people decide to evacuate, the more people will be at the shelters by November 14, 2010. Moreover, immediate evacuation is corresponding to a relatively steady increase of the IDPs. On the contrary, late evacuation generates a dramatic increase in IDPs at a certain time period. Statistical test indicates that the evacuation pattern of immediate evacuation and that of late evacuation is statistically significantly different (t(38)=−2.388, p<0.05). It encapsulates that when to evacuate is a significant factor to influence the evacuation. It is worth noting that the base model has similar evacuation pattern with that of late evacuation (t(38)=−0.694, p=0.694).

“How to evacuate” decision can be classified as adaptive behavior (self-evacuate behavior), non-adaptive behavior and leader-following. It appears that children, disabilities and elders are unlikely to perform adaptive behavior. However, based on the empirical distribution (Rahman, 2017), not all adults perform adaptive behavior, some perform leader-following. The sensitivity analysis is thus conducted to assess the effect of the “how to evacuate” distribution of the adults on the evacuation pattern. Two simulation settings, i.e. all adults performing adaptive behavior and 90 percent of adults performing adaptive behavior, were conducted. Figure 7 indicates that similar evacuation patterns among the three simulation settings are identical, indicated by statistical tests of t(38)=−0.055, p=0.956 for all adults performing adaptive behavior, and of t(38)=0.276, p=0.784 for 90 percent of adults performing adaptive behavior. It is also interesting to note that adaptive behavior results in the highest number of evacuee reaching the shelters. The more people performing adaptive behavior, the more people are able to evacuate at the shelters.

Moreover, Figure 7 indicates that the effect of evacuation behavior can be noticeably observed in the long-run (after November 5, 2010). It implies that the evacuation behavior affects the effectiveness of evacuation process – the IDPs reaching the shelters. It is then interesting to evaluate the time required for reaching the shelter. The next sensitivity analysis is hence comparing the three simulation settings using evacuation time as a performance indicator.

Figure 8 shows the average evacuation time for the three simulation settings. When using average evacuation time, a significant difference is observed between the base model and the all adults performing adaptive behavior model (t(38)=3.513, p<0.001), and non-significant difference between the base model and the model of 90 percent adults performing adaptive behavior (t(38)=−1.745, p=0.089).

Although the base model has already consisted of 99.83 percent adults performing adaptive behavior, the evacuation time is however not significantly different than that of 90 percent of adults performing adaptive behavior. Figure 8 further demonstrates that it requires 100 percent of adults performing adaptive behavior, corresponding to 69 percent of the total population, to significantly reduce average evacuation time. It implies that educating people to be able to self-evacuate is essential to reduce the evacuation time.

Summarizing, both “when to evacuate” and “how to evacuate” decisions are influencing factors in the course of evacuation. It is interesting to note that the earlier period of evacuation is highly affected by “when to evacuate” decision, whereas the later period is influenced by “how to evacuate” decision. Given the fact that “when to evacuate” decision is depending on hazard threshold and “how to evacuate” decision is depending on human-related factor (such as age, physical ability) and cultural-related factor (such as trust to cultural leader), the results imply that both human and cultural factors should be considered when modeling the evacuation and/or designing the evacuation plan.

Coordination in the simulation model refers to the interaction among aid organizations in the form of information sharing. Actors involved in post-disaster operations were government, military and aid organizations. Based on Indonesian Regulation No. 21 of 2008 on disaster management, government acts as the main actor in the humanitarian supply chain. The National Disaster Management Agency (BNPB) and the Regional Disaster Management Agency (BPBD) serve as the government agency that handles a disaster in the emergency response phase. The military was also an actor involved to provide major help such as communal kitchens and road repair. The military was highly coordinated within the internal sphere so that all actions were undertaken on the basis of instructions from the leadership center. Meanwhile, relief agencies/aid organizations such as NGOs, humanitarian organizations and some other profit companies were actors which help the government to reduce losses.

Unlike the government agency and military which are centrally coordinated within themselves, aid organizations work independently. Given that the number of aid organizations is quite significant during the 2010 Mount Merapi post-disaster, it is argued that coordination among aid organizations, and between a government agency and aid organizations, is required in favor of efficient relief operations. However, in reality, many aid organizations had worked independently. Coordination among aid organizations was absence. The scenario is therefore developed in order to evaluate the effectiveness of coordination among the aid organizations. Two scenarios are developed as follows:

• Scenario 1: no-coordination among the aid organizations. The aid organizations work independently. If each volunteer agent represents an aid organization, each volunteer makes its own decision in determining the destined shelter and delivering the aids. The scenario denotes, to some extent, decentralized coordination.

• Scenario 2: coordination is implemented by allowing information sharing among the volunteers on the shelter to be visited so that the aids can be delivered evenly to the shelters (horizontal coordination in downstream logistical activities). The scenario applies centralized coordination in which it is assumed that the information is centralized and the government is responsible for coordinating the volunteers. In addition to volunteer coordination, the government is also responsible for allocation of the aids between the aid centers (Maguwoharjo and UGM) so that the supplies can be allocated evenly to the two aid centers (coordination in upstream logistics activities). The coordination scenario is thus further detailed into two sub-scenarios. Consequently, Scenario 2a simulates the coordination among volunteers within an aid center in delivering the aids to the shelters. Scenario 2b simulates the coordination among volunteers within an aid center and the coordination between the two aid centers to allocate the supplies.

  Both Scenarios 1 and 2 address coordination mechanism, i.e. decentralized coordination and centralized coordination, respectively. On the other hand, some literatures have pointed out other operational strategies relevant to last-mile delivery, such as clustering strategy, using small-sized truck and pre-positioning strategy. Clustering strategy has been implemented in the commercial supply chain in order to improve response time, e.g. Schmitt et al. (2015). Comfort et al. (2004) have demonstrated that the spatial size influenced the duration of disaster response activities. The larger the size of the disaster area, the more time needed to meet the demand. For that reason, the regional division increased the efficiency of response operations (Comfort et al., 2004; Balcik et al., 2010). Second, relief operation is characterized by responsiveness particularly during the early period of after-math disaster. The responsive supply chain is normally using a small vehicle rather than a large vehicle in the last-mile delivery. The scenario is implemented by reducing the capacity of the vehicle by half, but increasing the number of vehicles as twice as much, so that the total available capacity of volunteers is similar to that of the base model. Third, preposition strategy has been discussed in many literatures such as Balcik and Beamon (2008) and found to be useful to anticipate sudden and surge demand during an after-math disaster. The simulation model assigns as much as 10,000 kcal as inventory at each aid shelter. Hereafter, the best scenario among the two coordination mechanism is combined with the operational strategies in Scenario 3 as the following.

• Scenario 3: combined coordination and operational strategies. Scenario 3a combines the strategic coordination with clustering of the service area. Scenario 3b applies the joint strategic coordination with the small-sized truck. Scenario 3c assesses the combination of coordination with pre-positioning strategy.

The simulation is run for 50 days (until December 15, 2010). The average of unmet demand based on 50 replications is then reported.

Figure 9 shows the simulation results of the non-coordination and coordination scenarios. During the earlier period of evacuation, the unmet demand presents and can be fulfilled by the volunteers. When the number of IDPs starts to increase dramatically, the available supplies are not sufficient to meet the demands. It is observable that scenario without coordination gives the worst performance when it comes to demand fulfillment. The increasing unmet demand pattern of non-coordination indicates under-supply at the shelters, an indication of ineffective aid delivery. This is due to supply redundancy at some shelters. Without information sharing, each volunteer delivers supplies without considering supplies delivered by other volunteers. It causes repeated effort in one shelter and gaps in services in other shelters. The result is supported by Aros and Gibbons (2018) who demonstrated that communication is corresponding to more effective response at the disaster site. Limited information sharing and communication was evidenced to be one of the important barriers against efficient and effective relief operations (Bharosa et al., 2010). In short, the coordination scenario performs much better than the non-coordination scenario. Furthermore, Akhtar et al. (2012) evidenced that the coordination among aid organizations is a necessity when organizations individually cannot effectively respond to disaster because coordination produces better competencies and resources, and thus more effective and efficient to respond.

While aforementioned previous studies have only highlighted the necessity of coordination based on empirical analysis, Figure 9 provides quantitative and dynamic assessment of the effect of coordination on the demand fulfillment. The figure specifically provides insights on the mechanism of coordination with respect to who should coordinate and where the coordination should be conducted. Through Scenarios 2a and 2b, it is demonstrated that coordination among volunteers at downstream operations, which has been much discussed in the previous studies, is not sufficient to handle uneven distribution and long-term demand fulfillment. Although Scenario 2a (coordination among the volunteers) has faster response as the disaster strikes than Scenario 2b, however, it fails to meet the demands in the long-run. It appears that the pattern of unmet demand is repeated (steady state) after November 25, 2010. It indicates that Scenario 2a fails to adapt with the dynamic change of the evacuation. The pattern seems unlikely to change unless there is another intervention on the system. It indicates that uneven supplies at the aid centers exist so that it will further constrain the delivery operations at the downstream. Conversely, Scenario 2b (coordination among volunteers and between aid centers) has one high peak on November 16, 2010, but demonstrates better performance in the long-term. It appears that the pattern of a high peak followed by decreasing unmet demands indicates an attempt of the supply chain to adjust and adapt its configuration in allocating and delivering aids to keep up with the dynamics of the IDPs at the shelters. It is therefore worth to note that coordination conducted at both upstream and downstream operations is required for sustaining demand fulfillment.

Figure 10 shows the simulation results of the combined strategy between coordination and other operational strategies of clustering strategy, small-sized truck strategy and pre-positioning strategy, respectively. The figure reveals that the combined strategies, in general, perform better performance than the coordination strategy alone.

Previous studies have proposed various combined strategies, such as Comfort et al. (2004) with clustering strategy, Kunz et al. (2014) with pre-positioning strategy and increasing disaster management capabilities. However, none of literatures has conducted comparison analysis of the effectiveness among the combined strategies. The present study therefore contributes to evaluate the effectiveness of the combined strategies, particularly for the context of volcano eruption of Mount Merapi.

Among the evaluated strategy combination, it appears that the combined strategy of coordination and clustering strategy performs best when it comes to long-term demand fulfillment. The clustering strategy, in which customers are grouped into a number of groups, has usually been implemented in the commercial supply chain to reduce response time and transportation cost. It is found that an increasing number of region operating as self-organizing units is corresponding to more efficient disaster response. It can be argued that clustering strategy reduces delivery lead time, and consequently providing better response. Furthermore, clustering strategy allows to embrace local knowledge which can further facilitate efficient communication.

Transportation plays a critical role with respect to lead time and meeting customer needs (Balcik et al., 2010) so that the choice of transportation mode is essential. Nevertheless, literatures exploring the transportation mode in humanitarian operations are very few. The most recent study, i.e. Maghfiroh and Hanaoka (2018) has evidenced that based on the optimization of last-mile delivery routing under uncertainty, the number of small-sized truck used in distributing aids is higher than the number of larger-sized trucks. The simulation results indicate that the usage of the small-sized truck seems to be promising, which is in line with Maghfiroh and Hanaoka (2018). Using two small-sized trucks is preferable than one large-sized truck because the small-sized trucks provide flexible accessibility (when roads are narrow), and flexible routing particularly when information gaps during last-mile delivery exist. Along with the increased demands, the required additional number of the small-sized trucks in relief distribution was significantly increasing by the increase of demand, while it was not the case for the large-sized trucks. However, it is interesting to notice that the combination of coordination and clustering strategy performs slightly better than that of coordination and using small-sized trucks.

The combined strategy with pre-positioning strategy (Scenario 3c) is able to meet the demands during early phase of evacuation, but fails to meet the subsequent demands. The pattern of pre-positioning strategy of the present study is actually similar to the pattern produced by Kunz et al. (2014) who examined the effectiveness of pre-positioning strategy using system dynamics approach. The present study demonstrates further that the pattern of the combined strategy with pre-positioning strategy (Scenario 3c) appears to be similar to that of the coordination strategy (Scenario 2b) unless that pre-positioning strategy can fulfill all the demands during the early phase of evacuation. It implies that the pre-positioning strategy only contributes to the 100 percent service level during the early days of relief operations (short-term impact). On the other hand, the pre-positioning strategy requires high investment due to the uncertainty when it comes to the location and the timing of the next affected areas. It is unlikely that the pre-positioning strategy would be able to meet the demands for long-fused disaster. Hence, the pre-positioning strategy should be combined with other strategies such as coordination which has been demonstrated in the present study for sustaining demand fulfillment.

Derived from the simulation results, insights are gained and managerial and policy implications are suggested as the following.

First, given that low hazard perception and leader-following behavior are the two factors responsible for late evacuation, the human and cultural factors are of importance in the evacuation process, and therefore should be considered in designing evacuation plan. The fact that only 20 percent of the population evacuated early indicates that effective early warning system is required to facilitate early evacuation, enabling the affected population to make decision to evacuate early. The earlier the evacuation, the more IDPs arrive at the shelters. Moreover, the early evacuation produces less dynamics of the evacuation pattern, and subsequently more predictable demands. Consequently, the extreme burst of increased demand, which results in variability in the supply chain, can be avoided. As the cultural aspect (i.e. trust to cultural leaders is strong particularly in Javanese community), the early warning system can be developed by not only technology through alert system, but also community empowerment through cultural/religion leaders. Likewise, the early evacuation should also be accompanied by self-evacuation behavior, independently capable of identifying the nearest shelters and the shortest routes toward the shelters. Specifically, at least 69 percent of the population (which is corresponding to 100 percent of adults in the simulation model) should perform self-evacuation behavior in order to significantly reduce the evacuation time. Furthermore, a local government may facilitate education/training on the importance of early evacuation and self-evacuation behavior, as well as, the standard operating procedure of hazard (including location of shelters and the evacuation routes) to cultural/religion/community leaders who then influence the rest of the population. It is believed that investing on social capital to handle hazard would facilitate effective community-based hazard management, which ultimately inhibits human casualties.

Second, coordination among volunteers at downstream operation is not sufficient to deal with uneven distribution and long-term demand fulfillment. This is due to that the capability of downstream operation to meet the demands is however constrained by the available supplies at the aid centers. It is therefore coordination at both upstream and downstream operations should be conducted simultaneously. Different actors such as donors and volunteers (NGOs) are involved in the coordination at the upstream and downstream operations, respectively. It implies that the government who has the authority could help to form the coordination among the different actors in humanitarian operations and take part as a coordinator in managing and coordination NGOs. Given the fact that NGOs may compete among themselves for a scarce resource, the presence of authority facilitates collaboration among NGOs. From the perspective of NGOs, the coordination would be beneficial as the NGOs can plan aid delivery better and ensure that aids have been delivered without wasted, which may eventually generate more resources from donors.

Third, the combined coordination and operational strategies are promising to overcome the challenge in the last-mile delivery such as uneven distribution. Coordination provides transparency, facilitating broad perspective rather than local perspective, so that better decision for resource allocation can be made, whereas the operational strategies such as clustering strategy, using small-sized truck, pre-positioning strategy, reduce lead time so that responsiveness can be improved. The combined strategies allow to achieve both efficiency and responsiveness at the same time. It appears that the combined strategy of coordination and clustering strategy seems to be the most effective and efficient strategy to be implemented for the context of the present study, especially Javanese community which is characterized by strong culture-bound community perfectly fits to support community-based hazard management. The combined strategy is also beneficial in covering large disaster areas. On the other hand, the combined strategy with pre-positioning strategy seems to be the least preferable option due to high uncertainty of the affected areas and high investment cost.

The present study has conveyed methodological and empirical contributions in the existing literatures in humanitarian context. With respect to the methodological aspect, the present study has developed an empirically founded agent-based model to represent the evacuation dynamics and last-mile relief delivery during post-disaster of the 2010 Mount Merapi eruption, which is still lacking in the literatures. The developed agent-based model has been verified, calibrated and passed a set of validation tests. The patterns resulted from the simulation have been able to reproduce the stylized facts of historical data with statistically non-significant difference, indicating the practical validity of the model representing the actual system. In other words, the developed model has met acceptable confidence level to be used as an experimental tool. With respect to empirical contributions, the present study has also provided an understanding on the mechanism of individual evacuation decision making at the micro-level underlying the evacuation dynamics at the macro-level.

Moreover, the findings highlight that both human and cultural factors (i.e. low hazard perception, evacuation behavior and trust to cultural leader) are of significance in the evacuation process, and therefore should be pondered. Early warning system and community-based hazard management are required to facilitate early evacuation and self-evacuation behavior for effective and efficient evacuation process. On the other side, coordination in relief operations is a must and should be carried out at both sides of upstream and downstream of last-mile delivery operations. To enable the coordination at both sides, the local government plays a role as a coordinating body. The combined strategy of coordination with other operational strategies, particularly the combination of coordination and cluster strategy, seems to be promising toward more efficient and effective relief operations.

Although the model allows for quantitative evaluation of the effectiveness of relief strategies, some limitations of the model should be highlighted. The mechanism of reloading to the aid centers has not been modeled in detail because the simulation used time-based reloading (i.e. one day) to represent the reloading process. In reality, volunteers may take more time to replenish the aids due to traffic congestions, or unavailable supply. The avenues for future research could enhance the developed empirical agent-based model to include traffic congestion during volunteers returning to the aid centers. The present study disregards the urgency of the demand points, so that the urgency of demand points is also worth to explore as it may lead to increases of social benefit. Last but not least, the present study focuses on last-mile relief distribution. The model could hence be enhanced by incorporating sourcing operations (such as collaborative procurement) for a better understanding of the overall supply chain system.