1 Introduction

Tsunamis are one of the deadliest natural hazards in the world and have resulted in more than 200,000 fatalities and economic losses greater than 14 trillion U.S. dollars in the last 20 years (Rossetto et al. 2007; Charvet et al. 2014; Fukutani et al. 2015). One of the most important actions to mitigate tsunami risk is to develop effective tsunami evacuation plans (Taubenböck et al. 2009; Mas et al. 2012; Charvet et al. 2014; Chiou et al. 2015; Kamat 2015). In general, current tsunami evacuation plans are developed based on (1) simplified evacuation routes and evacuees' movements/behavior and (2) detailed evacuation modeling representing evacuation routes and evacuees' movements in a realistic manner.

The simplified approach typically considers evacuees moving directly to evacuation areas designated a priori. In this approach, the evacuation route is defined as a straight line connecting the evacuee to the evacuation point without considering avaliable transportation networks and obstacles (e.g. slope of the evacuation route; Taubenböck et al. 2009; Mas et al. 2012). Subsequently, the simplified approach calculates the evacuation time based on the shortest horizontal route to reach evacuation areas (e.g. higher ground) and hence, may be considered as a static evacuation approach (Muhammad et al. 2017; De Risi et al. 2018). On the contrary, the agent-based evacuation modelling can simulate dynamically the escape of individual evacuees to safe areas by accounting available transportation networks and obstacles. The dynamic approach also considers the realistic evacuation routes used by the agent to the evacuation point and time-dependency of the accessible routes prior to tsunami arrival (Zarboutis and Marmaras 2007; Chen and Zhan 2008). Such a dynamic approach has been implemented to develop realistic evacuation plans for various hazard mitigations (e.g. flood, urban fire, and dynamic urban growth; Zhang et al. 2015; Augustijn et al. 2016; Nishino 2019).

An agent-based model can be used for developing realistic tsunami evacuation scenarios. However, in such evacuation studies, a rigorous hazard analysis considering the uncertainty of future tsunami events is often neglected (Wang et al. 2016; Muhammad et al. 2017; Makinoshima et al. 2018) as it may require substantial computational efforts. Typically, open areas with an elevation higher than 10 m are designated as tsunami evacuation areas without assessing possible inundation scenarios thoroughly (Lämmel et al. 2010; Di Mauro et al. 2013). Moreover, a deterministic worst-case scenario is generally adopted, leading to potentially inaccurate or inadequate tsunami evacuation plans (Mostafizi et al. 2017; Takabatake et al. 2020; Kitamura et al. 2020; Faucher et al. 2020; Kubisch et al. 2020).

In this study, the evacuation is simulated using the agent-based modeling package MATSim (Balmer et al. 2009; Lämmel et al. 2009). MATSim requires the identification of links, nodes, and agents within a transportation network (left panel of Fig. 1a). Final destinations in evacuation simulations are represented by evacuation nodes of the network (red dots in the left panel of Fig. 1a) and are further connected to virtual evacuation links (dashed line in the left panel of Fig. 1a) and nodes (black dots in the left panel of Fig. 1a) which have infinite evacuation capacity (Lämmel et al. 2009, 2010).

Fig. 1
figure 1

a Agent-based tsunami evacuation modeling in MATSim and simulation scenarios, b study areas, c illustration of line capacity determination for motorbike (left panel) and pedestrian (right panel) in agent-based simulation, d distribution of agents in Padang, e evacuation modeling parameter set up

Full size image

Recently, probabilistic earthquake source models have been successfully used to capture the uncertainty of tsunami hazard in Mexico, Indonesia, and Japan (Goda and Song 2016; Grezio et al. 2017; Griffin et al. 2017; Goda et al. 2017; Goda 2019) and have been further implemented to quantify the tsunami risks (Muhammad et al. 2017; Goda 2019). However, no tsunami evacuation plan using such probabilistic models have been implemented. Subsequently, in this paper, a new way of propagating tsunami hazard uncertainty into evacuation scenarios is proposed. The methodology approaches the problem in a holistic manner considering (a) a stochastic source model, (b) numerical tsunami inundation simulations, and (c) an agent-based simulation for modeling the movements of evacuees. Implementing such probabilistic modeling is useful for assessing current tsunami evacuation approaches and for developing a new holistic framework for tsunami evacuation. Specifically, each simulated tsunami scenario allows distinguishing safe and unsafe areas clearly and identify potential evacuation points; for each stochastic inundation map, the evacuation nodes are updated according to the varying safe zones (red dots in Fig. 2b, e, h). This approach, for the first time, allows propagating the tsunami hazard uncertainty into the evacuation models (Taubenböck et al. 2009; Mas et al. 2012; Muhammad et al. 2017; De Risi et al. 2018).

Fig. 2
figure 2

Probabilistic tsunami hazards: a, d, g earthquake slip model, b, e, h inundation model, c, f, i tsunami arrival time for magnitudes Mw 9.0, Mw 8.75, and Mw 8.5, respectively

Full size image

As a case study, a region of Teluk Bayar in Padang (Fig. 1b), Indonesia is selected. It is exposed to a significant tsunami risk from the Mentawai segment of the Sunda subduction zone (Nalbant et al. 2005; Borrero et al. 2006) that may affect 6622 people who live there (BPS 2018). Several essential infrastructures, e.g. oil company, harbor, and main fish market, are located in this area with a typical road width of 6 m divided into two lanes, i.e. 3 m (see Fig. 1c; Hoppe and Mahardiko 2010; Alhadi 2014). Local authorities of Padang suggested two evacuation destinations for tsunami evacuation in this area: Air Manis Street and Bypass (red dots with white-dashed circle in Fig. 1d; Alhadi 2014); these evacuation points have been used during evacuation drills since 2010 (Hoppe and Mahardiko 2010; Alhadi 2014; Yosritzal et al. 2018). However, no tsunami evacuation shelters (TES) have been constructed and assigned in this area compared to other areas in Padang (e.g. more than ten TES have been built in the west of Padang region). Hence, it is essential to evaluate current existing tsunami evacuation plans in Teluk Bayur rigorously.

In this study, three cases are considered (right panel of Fig. 1a): (a) the simplified evacuation using a static model (simplified case), (b) the agent-based tsunami evacuation model considering the two official evacuation points (existing case), and (c) the agent-based tsunami evacuation model using stochastic tsunami inundation maps (optimized case). The optimized case is more rigorous and is used as a benchmark to evaluate current tsunami evacuation plans modeled through the simplified and existing cases. The starting time of evacuation is set to 9.00 a.m. (daytime; the busiest time). Moreover, two evacuation modes (i.e. the way evacuees moves/are transported) are considered independently: (a) pedestrian and (b) motorbike. These two modes are the two modes suggested by local authorities for short and long evacuation paths, respectively (Hoppe and Mahardiko 2010; Yosritzal et al. 2018), and are assessed independently as it is easier to draw practical conclusions.

To carry out such integrated stochastic tsunami evacuation assessments, first, the stochastic earthquake source models are generated. Second, the Monte Carlo tsunami simulation is carried out to develop stochastic inundation maps. Then, the agent-based tsunami simulation is implemented. Finally, the results in terms of tsunami evacuation time and the number of affected people are obtained. Eventually, comments and conclusions are drawn.

2 Materials and methods

In this paper, a new way of propagating tsunami hazard uncertainty into evacuation scenarios is proposed. The methodology approaches the problem in a holistic manner by considering (a) a stochastic source model, (b) the numerical tsunami inundation simulation, and (c) the agent-based simulation for modelling the movements of evacuees. The proposed methodology is compared with the simplified approach and the existing deterministic scenario that is considered by the local government.

2.1 Stochastic earthquake source models and tsunami simulation

The earthquake source models are generated stochastically by considering regional seismological characteristics (Muhammad et al. 2016) and by calculating earthquake source parameters including geometry (fault length and width), slip statistics (mean slip, maximum slip, and Box–Cox parameter), and spatial heterogeneity parameters (correlation length along strike direction, correlation length along dip direction, and Hurst number). Those parameters are calculated using the effective dimension analysis method (Mai and Beroza 2000), Box–Cox analysis (Goda et al. 2014), and spectral analysis (Mai and Beroza 2002). In this study, a total of 900 earthquake slip models are generated by considering 300 scenarios for moment magnitudes (Mw) of 8.5, 8.75, and 9.0 (Fig. 2a, d, g).

The initial deformation of the seafloor due to earthquake fault rupture is calculated for a given source model by considering both horizontal and vertical displacements of the seafloor (Okada 1985; Tanioka and Satake 1996). Tsunami propagation and inundation are then carried out by solving nonlinear shallow water equations with run-up (Goto et al. 1997). For the tsunami propagation, the combination of a high-resolution bathymetry (i.e. < 10 m) dataset and GEBCO2014 (http://www.gebco.net) is used. Moreover, the DEM adopted to carry out the tsunami simulation in this study was developed from GPS measurements and high-resolution satellite imagery with a resolution of < 10 m. The vertical elevation error for the local DEM and bathymetry datasets are 0.2 m and 0.5 m, respectively (Taubenböck et al. 2009; Schlurmann et al. 2010). Such a local DEM dataset allows a more accurate tsunami inundation model with respect to the global DEM datasets (SRTM1 and GDEM2; Muhammad et al. 2017; Muhammad and Goda 2018).

However, instead of using a 10-m grid cell, a 50-m grid is considered to run the tsunami simulation. This choice was made on the basis of the following: (i) a large number of tsunami simulation (i.e. 900 models) need to be run to obtain the stochastic inundation maps. A 10-m grid cell results in a computational time that is too demanding (i.e. 2 h vs 2 days for each tsunami simulation model for the 50-m vs. the 10-m grids); (ii) the 50-m is found to be sufficient to produce results as accurate as the 10-m grid specifically in defining the inundation areas for a relatively flat topography region like in Teluk Bayur, Padang; and (iii) the 50-m grid produces very close tsunami arrival time at the same point compared to the 10-m resolution (< 1 min) (Griffin et al. 2015; Muhammad and Goda 2018). Therefore, the computational grid used to run the tsunami simulation consists of four nested grids with increasing resolution: 1350, 450, 150, and 50 m. The duration and time-step of the tsunami simulation are 2 h and 0.5 s, respectively. Such a grid resolution is only used for tsunami simulation that is carried out separately with respect to the agent-based evacuation modelling and hence, the evacuation time of the agents is not affected by the change of tsunami simulation grid.

For each tsunami slip model, a tsunami simulation is carried out to generate a tsunami inundation map (Fig. 2b, e, h) and tsunami arrival times at agent’s locations in the network (Fig. 2c, f, i; see also Figs. 5, 6 and 7 of the “Appendix” for detail). The tsunami arrival time is defined as the time for which the inundation depth at the location of interest exceeds 0.3 m; it has been demonstrated that the evacuation might become difficult for this level of water depth (Mas et al. 2012). A total of 1800 tsunami evacuation time maps, i.e. three magnitude scenarios × 300 earthquake slip models × two evacuation modes, is generated. These evacuation maps are then used for the assessments of the evacuation plans.

2.2 Agent-based evacuation modelling

The agent-based evacuation modelling is carried out using the MATSim package. The surveys carried out after the 2009 Padang earthquake confirmed that about 70–80% of the people in Padang evacuated on foot and by motorbikes; the remaining population used bicycles and cars (Hoppe and Mahardiko 2010; Yosritzal et al. 2018). The evacuation during the 2009 event clearly demonstrated that the use of motorbikes and cars caused heavy congestion on the roads. Therefore, local authorities have suggested using either foot for short-route evacuation (e.g. to tsunami shelters) or motorbikes for long-route evacuation (Hoppe and Mahardiko 2010; Yosritzal et al. 2018). Based on these suggestions, in this study, two types of evacuation modes are considered: pedestrian or motorbike. The two modes are studied independently as it is easier to draw conclusions considering the two cases separately.

In MATSim, the speed can be defined both at link and agent levels (Balmer et al. 2009; Lämmel et al. 2009). In this study, the speed is implemented at link level as the two modes are considered separately, and the convergence to stable results is more straightforward to achieve [see simple examples in De Luca et al. (2021)]. The link in MATSim represents the road used by agents to evacuate. If a road is two-way, it is defined as two links, i.e. one per direction (Balmer et al. 2009; Lämmel et al. 2009). For the pedestrian evacuation, each evacuee is an agent. Conversely, for the motorcycle transportation mode, two to three evacuees are defined as an agent, especially if the evacuees belong to the same family/household or are located in middle/elementary schools. For instance, if one building consists of 5 people, then two motorcycles that are available in that building carries two and three people on each motorbike.

The agents are distributed based on the population data of each district including the number of average household size (i.e. 5 people) and building footprints in Teluk Bayur, Padang (Fig. 1d; BPS 2018). Population data are extracted from the Statistics Indonesia containing detail population information up to a district level (BPS 2018). Integrating these data with the building footprints from Open Street Map (https://www.openstreetmap.org), the agent's distribution is determined by considering a daytime period (i.e. 9:00 a.m.). The final building occupancy for residential buildings varies between 1 and 3 people, as most agents are going out in the daytime. The number of people at school, office, and market is generated based on the number of students (including the teacher/staff), workers, and people at the market, respectively (BPS 2018; Yosritzal et al. 2018). The remaining population is assigned to non-residential buildings.

2.3 Evacuation scenarios

In terms of evacuation points, the evacuation nodes for the three considered cases are modelled in a different manner depending on the cases. For the simplified case, all points in the areas outside the inundated zone are considered as safe (i.e. points beyond the red curve in the right panel of Fig. 1a). The safe area (i.e. the red curve in Fig. 1a) changes for each stochastic inundation map. The evacuation time (\(t\)) for this simplified case is calculated as a function of the distance (\(S\)) over the agent's speed (\(v\)), i.e. \(t=S/v\). The distance S is defined as the closest distance (along a strainght line) between the agent's location and the closest point in the non-inundated area, while the agent's speed is set to 1.66 m/s or 5.66 m/s for pedestrians and motorbikes, respectively.

For the second case considered in this study, two evacuation points (i.e. Air Manis and Bypass streets) are used for the existing evacuation scenario (red dots with white-dashed circle in Fig. 1d). For this case, only one evacuation plan is defined; therefore, only one hazard scenario (i.e., worst case) is considered neglecting the hazard uncertainty and acknowledging the current practice of using deterministic worst-case inundation maps. The chosen inundation model for the two evacuation points case corresponds to a slip model producing the largest inundation in terms of inundation level and extension inland. The adopted inundation model, also used for the assessment of the existing evacuation case, is consistent with the tsunami inundation map from the Germany Indonesia Tsunami Early Warning System (GITEWS) project that is used to develop the current local tsunami evacuation plan (Schlurmann et al. 2010).

Finally, the optimized case adopts various evacuation points located in the non-inundated zone that are unpaved and connected with the accessible links (see the right panel of Fig. 1a, e). The unpaved networks to those various evacuation points are not urbanized and prepared as proper evacuation area and are available from the beginning of evacuation simulation (i.e. 9:00 a.m.). In reality, the evacuee may continue to move to a safer area (i.e. non inundated area) if the tsunami inundates those evacuation points. However, in this agent-based evacuation simulation, those defined evacuation points are used as the final destination of each evacuee. Moreover, the extra evacuation points in the optimized case do not correspond to officially indicated safe zones, but they can be considered as safe on the basis of comprehensive tsunami inundation simulation results (Muhammad et al. 2017; Muhammad and Goda 2018). The assessment of the possibility of using those points for evacuation may provide valuable recommendations for the local authorities as to how more effective tsunami evacuation plans shall be developed in the future.

Except for the two existing evacuation points (i.e. Air Manis street and Bypass), the capacity of the extra evacuation points defined for the optimized case has been considered limited, i.e. a limited number of agents can be  accomodated in those areas. The threshold capacity is defined based on the capacity of the lane (i.e. the dimensions of the road) and of the space needed by the agents. Figure 1c illustrates the lane capacity for both pedestrian (right panel) and motorbike (left panel) corresponding to the real network (i.e. solid line in Fig. 1a). For the two considered evacuation modes, an area of \(\sim 0.5\, {\mathrm{m}}^{2}\) and \(\sim 3\, {\mathrm{m}}^{2}\) is needed for a single pedestrian and motorbike, respectively (Lämmel et al. 2009). Considering a length and width of a motorbike of ~ 3 m and ~ 1 m, respectively (Hoppe and Mahardiko 2010; Yosritzal et al. 2018), and considering the width of the road equal to 3 m, the lane capacity can be defined as 3, and the threshold for the motorbike capacity at the evacuation links can be set equal to 30 (i.e. represented by the dashed line in Fig. 1a). This number corresponds to a heavy congestion (i.e. a queue of 10 lines of 3 motorbikes), a situation that may lead to a rerouting of the evacuees to move to other evacuation points. The space for the pedestrian in each lane is at least 6 times smaller than the space of the motorbikes (i.e. \(\sim 0.5\, {\mathrm{m}}^{2}\) vs. \(\sim 3\, {\mathrm{m}}^{2}\); therefore, the capacity threshold for the pedestrians may be set to 200 (i.e. about 6 times larger than the motorbike).

2.4 Network

About the transportation network, the length of the road (link) varies from ~ 50 to ~ 500 m with a width of 6 m for two lanes (Hoppe and Mahardiko 2010; Yosritzal et al. 2018). Therefore, the link (one lane) is modeled with a width of 3 m. As illustrated in Fig. 1c, distances from one agent to another agent for the pedestrian and motorbike modes in the real network (solid line in Fig. 1a) are 0.6 m and 1 m, respectively. Subsequently, the maximum number of agent that can have access to the link at the same time for the pedestrian and the motorbike modes are 5 and 3 agents, respectively. For evacuation modeling, the following constraints are considered for the network. First, the maximum flow speed is restricted depending on the mode of the agents, i.e. 1.66 m/s and 5.66 m/s for pedestrian and motorbike modes, respectively (Hoppe and Mahardiko 2010; Muhammad et al. 2017; Yosritzal et al. 2018). Second, the maximum flow capacity (\(FC\)) is limited and is calculated using the following equation:

$$FC=w\times {C}_{max}$$
(1)

where \(w\) is the width of the link and \({C}_{max}\) is the maximum flow capacity per unit width based on the mode used for evacuation simulation. For instance, the maximum flow capacity for the pedestrian mode is 1.3 \(\mathrm{persons}/\mathrm{m}/\mathrm{s}\). Such a value is determined from Weidmann's fundamental diagram (Weidmann 1993) that defines a relationship between density and velocity (Lämmel et al. 2009). For the motorbike mode, the maximum flow velocity is set to 4 \(\mathrm{persons}/\mathrm{m}/\mathrm{s},\) taking into account that the flow speed of motorbikes is about 3 times faster than pedestrians. Such a limitation restricts the agents to leave the link when the maximum flow capacity is reached. Finally, the link capacity (\(\mathrm{persons}/{\mathrm{m}}^{2}\)), i.e. storage capacity (\(SC\)) is defined to limit the number of agents present on the link and is calculated based on the areas (\(A\)) and maximum density per unit areas (\({D}_{max}\)) (see Eq. 2; Lämmel et al. 2009). Consequently, if the storage capacity is reached, the agent is not allowed to proceed, but it is forced to spill back, leading to congestion (Lämmel et al. 2009). \({D}_{max}\) for pedestrian and motorbike are 5.4 \(\mathrm{persons}/{\mathrm{m}}^{2}\) and 0.4 \(\mathrm{motorbike}/{\mathrm{m}}^{2}\) (Lämmel et al. 2009).

$$SC=A\times {D}_{max}$$
(2)

Moreover, the speed of the pedestrian model in the simulation depends on the slope of the corresponding links and thus is updated accordingly depending on the position of the agent and the terrain topography. Tobler's hiking function is adopted to change the speed of the link (Tobler 1993). A total of 250 iterations are performed in MATSim simulation that is sufficient to reach an equilibrium condition in generating optimal plan for each agent movement (Lämmel et al. 2009, 2010; Wood et al. 2020). Wood et al. (2020) have emphasised that the actual evacuation is likely to be reflected in the iteration range #5 to #100. This aspect can be relevant to evaluate how preparing local communities to evacuation can improve the evacuation times (e.g., comparing optimal vs. sub-optimal evacuation scenarios). Wood et al. (2020) have considered 100 iterations to exclude results of the initial iterations with unlikely estimate of evacuation clearance time and to account on the network knowledge of the evacuees and their capability of seeking out additional information before and during the travel. In this study, the number of iterations assumed relies on a good preparation of the population.

Two movement strategies representing a motion plan of each agent (Balmer et al. 2009; Lämmel et al. 2009) are implemented and weighted: ReRoute (10%) generating a new plan in each iteration with a new evacuation route based on the information of experienced travel times from the last run and BestScore (90%) producing a new plan in each iteration based on the best score. Moreover, the scores representing distance and time are also adopted to define effective evacuation routes [see Lämmel et al. (2010) for the details].

2.5 Evacuation time

The tsunami evacuation time is defined as a sum of the initial reaction time (IRT) and the actual evacuation time of an agent. The IRT represents the starting time of evacuee's response, whilst the actual evacuation time is the needed time for the agent to evacuate to a safe zone. In general, the IRT for Indonesia takes 10 to 15 min (Taubenböck et al. 2009; Hoppe and Mahardiko 2010; Yosritzal et al. 2018). In this study, the quickest response time (i.e. 10 min) is adopted. A 10-min initial reaction is reasonable as the average initial reaction for people in Indonesia and more specifically, in Padang, is about 3–5 min (Hoppe and Mahardiko 2010; Yosritzal et al. 2018). The remaining 5 to 7 min account for the time needed to people to get access to the network link. This is relatively quick considering that the location of the agents to the links is close (typically less than 50 m).

2.6 Number of affected people

To understand the effect of the variable tsunami evacuation plans on the population, the number of affected people is also calculated. To evaluate the number of affected people, the tsunami evacuation time and the tsunami arrival time at each agent's location in the network during the evacuation are compared.

The agent \(j\) at any node \(i\) of the evacuation path \({K}_{j}\) is considered to be trapped by the tsunami if the tsunami arrival time at node \(i\) (\({T}_{i}\)) is less than the agent's arrival time at node \(i\) (\({t}_{ij}\)) (Lämmel et al. 2009). Therefore, the agent \(j\) is considered as affected during the evacuation if the minimum of the differences between \({T}_{i}\) and (\({t}_{ij}\)) for all the nodes of the evacuation path \({K}_{j}\) of the agent \(j\) is less than zero:

$$\underset{\forall \in {K}_{j}}{\mathrm{min}}\left[{T}_{i}-{t}_{ij}\right]<0$$
(3)

By analyzing many tsunami hazard scenarios and associated evacuation plans, the estimation of the number of affected people becomes a stable metric that can be used as a robust decision variable. Such risk assessment results are useful for developing effective tsunami evacuation maps in coastal regions.

3 Results and discussion

The comparisons of the tsunami evacuation times for the three different cases and for the two transportation modes are firstly presented. The uncertain tsunami evacuation time using probabilistic inundation models is also discussed. Finally, the number of affected people for the three evacuation cases and different earthquake magnitude scenarios are examined to understand how different tsunami evacuation approaches affect the coastal communities in Teluk Bayar in terms of potential casualties.

3.1 Tsunami evacuation time

Figure 3 shows tsunami evacuation times for the three evacuation cases and for the two types of evacuation modes considered in this study: pedestrian (Fig. 3; left panel) and motorbike (Fig. 3; right panel) that are simulated separately. For the existing evacuation plan, there is no need for performing any tsunami hazard simulations. For the simplified and the optimized cases, the evacuation nodes change with the tsunami hazard; only the 50th percentile of the 300 evacuation times obtained changing the evacuation nodes according to the 300 inundation scenarios for Mw 9 is presented as a representative result. Complete results of tsunami evacuation maps for the other magnitude cases with different percentiles are presented in Figs. 8, 9, 10 and 11 in the “Appendix”.

The maximum tsunami evacuation time for the pedestrian mode in the simplified case is about 20 min, whilst it is less than 13 min for motorbike. On the other hand, the evacuation time for the existing and the optimized cases are longer (> 30 min) with respect to the simplified case (see Fig. 3c–f vs. Fig. 3a, b). As far as evacuation mode is concerned, since the motorbikes' speed (5.6 m/s) is faster than pedestrians' speed (1.06 m/s), for the majority of the cases the motorbike agents are quicker to evacuate than the pedestrian agents, for both optimized and the existing cases (see Fig. 3e vs. Fig. 3f). However, some of the motorbike agents have longer evacuation times than pedestrians (> 20 min; see Fig. 3d, f), especially those who are located close to the extra evacuation points. This is due to the capacity limitation of the evacuation destination to 30 motorbikes (i.e. \(\sim\) 75 people) for the motorbike mode. In the same locations, the pedestrian agents show relatively quicker evacuation time (< 30 min) since the capacity of the evacuation destination is larger (i.e. 200 people).

Fig. 3
figure 3

Tsunami evacuation time considering the 50th percentile of the Mw 9.0 simulation for pedestrian (left panel) and motorbike (right panel): a, b simplified case, c, d existing case, e, f optimized case. Box plots of the 300 cases using the optimized and the existing cases with pedestrian mode (g) and motorbike mode (h)

Full size image

Figure 3g, h show the distributions of evacuation times for the pedestrian-mode and motorbike-mode, respectively. The grey markers are the results of the 300 agent-based evacuation simulations for the Mw 9.0 optimized case (i.e. written as “Simulation—Optimized Case” in Fig. 3g). The number of people with a specific range of evacuation time (e.g. 10–15 min) are counted for each of the 300 evacuation simulation results and plotted in the figure. Moreover, the box plot presents the variation of the number of people for each category of evacuation time, demonstrating in a clear manner the effect of the uncertainty of future tsunami events. Such a representation helps capturing the uncertainty of the tsunami hazard from the evacuation perspectives and may provide with a better understanding of the possible evacuation time when the significant tsunamigenic event occurred in the study area. Results show that for the evacuation time shorter than 15 min, the pedestrian case shows higher variability ranging from 2500 to 6000 people (Fig. 3g). The variability decreases with the increasing of the evacuation time (except for the evacuation times exceeding 30 min).

The comparison between the optimized and the existing case shows that, for short evacuation time (e.g. 10–15 min), the existing case leads to an underestimation of the number of pedestrians who are successful in evacuation (see red circles in Fig. 3g, h). On the contrary, the number of people having longer evacuation times (> 20 min) is overestimated for the existing case, in particular for the pedestrian mode.

3.2 Affected people

The current tsunami evacuation planning approaches are evaluated in terms of the number of affected people by referring to the optimized case. The number of affected people is calculated by considering different magnitude scenarios and evacuation mode types. The number of affected people offers a significant additional piece of information. In fact, although the evacuation times of the existing and the optimized cases are relatively similar (Fig. 3c–f), especially for the motorbike case, the influence of extra evacuation points using the optimized case is proven significant when the number of affected people by the tsunami is taken into account. Figure 4a–d present the number of affected people for the existing and optimized case, respectively. The optimized case in Fig. 4c, d is developed from the same inundation model used to create the existing case. The figures clearly show that the number of affected people significantly decreases when extra evacuation points are considered in the simulation. Specifically, for the pedestrian mode, the number of affected people reduces from 1642 to 253, passing from the existing to the optimized case. Such a reduction is also found for the motorbike mode (i.e. the number of affected people for the existing and optimized cases are 1321 and 251 people, respectively).

Fig. 4
figure 4

Plot of affected people by the tsunami for pedestrian (left) and motorbike (right) modes for the existing (a, b) and optimized (c, d) cases. Probability of the number of affected people using two types of modes: e pedestrian and f motorbike with different magnitude cases: Mw 9.0 (red), Mw 8.75 (black), and Mw 8.5 (blue) along with the existing deterministic case (magenta line)

Full size image

Repeating the calculation of number of affected people for all magnitude case, it is possible to calculate the probability that a given number of people is affected (i.e. probability of tsunami arrival time is smaller than the escape time—\(P({T}_{tsunami}\le {T}_{escaping})\)). Such a probability is derived using the frequentistic approach from the simulation procedure for all the considered evacuation cases and modes (see Fig. 4e, f). For the simplified case, less than 100 people is affected (i.e. < 10% of the other cases) by the tsunami and hence, can not be clearly distinguished from the other two cases and therefore can not be captured through the plot of the probability that a given number of people is affected (i.e. Figure 4e, f). Subsequently, the results from the simplified case are excluded from the calculation of the probability of having a given number of affected people.

The figures clearly show that the probability of having a greater number of the population affected escalates with the increase of the earthquake magnitude. For instance, for the case of pedestrian evacuation, the probability of having 50 people affected by the tsunami for Mw 8.5 is only ~ 1%, and it increases up to ~ 4% for Mw 8.75 and up to ~ 20% for Mw 9.0. The probability of having a given number of affected people for the pedestrian case is smaller than the motorbike case (e.g. ~ 20% for having more than 50 affected people for the pedestrian case vs. ~ 50% for the motorbike case). This is mainly due to the capacity limitation at the extra evacuation points and the traffic congestion occurring at the network level. Moreover, for the simplified case for all three magnitude scenarios, less than 100 people are affected by the tsunami compared to more than 1000 affected people for the other cases. This highlights how unrealistic the simplified case can be, and hence, it is not recommended to be adopted for developing a tsunami evacuation plan. The current local evacuation routes defined by the local stakeholders (i.e. existing case) produce a larger number of affected people of about 23% of the total population for both pedestrian and motorbike modes.

3.3 Discussion

Developing an effective tsunami evacuation plan in coastal environments is a complex process. It is essential to take into consideration the uncertainty of the hazard of future tsunamis and realistic evacuation processes. In this study, the uncertainty of future tsunami hazard is taken into account by adopting a probabilistic inundation model (De Risi and Goda 2016; Muhammad et al. 2017; Muhammad and Goda 2018). Three magnitude values (i.e. Mw 8.5, Mw 8.75, and Mw 9.0) are used to generate 900 probabilistic inundation maps; this large number of simulations allow quantifying and visualizing the uncertainty of future tsunami events. In previous studies, only the worst-case scenario was used for the tsunami evacuation design and planning (Taubenböck et al. 2009; Muhammad et al. 2017; De Risi et al. 2018).

It has been observed that increasing values of magnitude lead to an increased number of affected people (Fig. 4e, f). The results also highlight that even for the Mw 8.5, the coastal community can be severely affected by the tsunami. For the case of Padang, a significant earthquake with ~ Mw 9.0 is anticipated based on the paleo-geodetic and paleo-tsunami studies (Nalbant et al. 2005; Borrero et al. 2006). Results highlighted that risk in terms of the number of affected people due to the Mw 9.0 events (Fig. 4e, f) may reach up to more than 1000 people for the optimized case. Therefore, it is imperative to develop comprehensive mitigation strategies to reduce the number of fatalities in Padang.

The effectiveness of the current tsunami evacuation plan approaches can be evaluated quantitatively from the results in terms of tsunami evacuation time (Fig. 3) and the number of affected people (Fig. 4). In the simplified case, no agent-based simulation is implemented, and hence, the interaction among the agents is not taken into account. Such an approach results in less number of people (< 100 people) affected by the tsunami in comparison to the agent-based evacuation modelling (i.e. the optimized case; see Fig. 4e, f). The two routes that have been assigned as the evacuation points by the local government (i.e. Air Manis and Bypass streets) cause longer evacuation times (> 30 min) and, consequently, about 1500 people may be affected by the tsunami (~ 25% of the total population). Therefore, additional evacuation points need to be made available to the coastal community provided that they are safe on the basis of stochastic hazard simulation. Moreover, by adding extra evacuation points, the evacuation on foot may be sufficient for the evacuee to arrive to the safe areas. Adopting the pedestrian mode can reduce the vehicle demand (i.e. less number of people using motorbike or car during evacuation) that is essential to reduce the possibility of traffic jams and to speed up the evacuation time (Henry et al. 2017; Wood et al. 2016, 2018, 2020).

The results of the optimized case highlight realistic evacuation scenarios through spatial variations of tsunami evacuation time (Fig. 3e, f). For the optimized case, the evacuation times are no longer uniform from the coast to the higher elevation because of the locations, the availability of evacuation links, and the interaction among agents (unlike the simplified case; Fig. 3a). As presented in dot-dash-magenta rectangle areas in Fig. 3e, the red circles showing longer evacuation times (> 30 min) are due to the congestion during the evacuation process. Such realistic evacuation phenomena cannot be captured through the simplified case. In addition, the 50th percentile of the optimized case for the pedestrian mode shows that two types of areas are highly prone to tsunami: (1) the zone closer to the coast (e.g. oil company and harbor; see dot-dash-magenta rectangle areas in Fig. 3e) and (2) high-density population areas (e.g. market; see dash-magenta areas in Fig. 3f).

Finally, as also observed by Hoppe and Mahardiko (2010) and Yosritzal et al. (2018), shorter tsunami arrival times in the areas close to the coast (< 20 min) and heavy congestion during the evacuation in the high-density population areas pose a significant risk to the coastal communities and therefore it is imperative to conduct a proper assessment.

4 Conclusions and recommendations

This paper proposed a novel computational framework to develop effective and reliable tsunami evacuation plans by integrating an agent-based modelling and stochastic tsunami inundation maps. A total of 900 stochastic inundation maps generated from three magnitude scenarios were used to define extra evacuation points. Such evacuation points are updated according to time-varying inundated zones that are produced from the stochastic tsunami simulation. The proposed model was further adopted to evaluate existing tsunami evacuation models in terms of the number of affected people by the tsunami.

This study demonstrated that the agent-based tsunami evacuation modelling together with the probabilistic inundation scenarios (optimized case) is an effective tool to develop tsunami evacuation plans as:

  • It implements a realistic tsunami evacuation plan by focusing on each agent's movement and interaction with other agents; none of these realistic features are considered in the simplified case.

  • It allows the identification of the most suitable and safer evacuation points. From the urban planning perspectives, in particular for the deterministic tsunamigenic event that has a return period of more than 200 years (Nalbant et al. 2005; Borrero et al. 2006; Muhammad et al. 2016), a large number of evacuation points is not always available. The proposed framework, considering the uncertainty of tsunami risk (e.g. evacuation time and the number of affected people) may provide stakeholders with a better understanding to define effective evacuation points. Therefore, local communities will access to multiple optimized evacuation options in comparison to the only two evacuation routes suggested by the local government.

  • It allows stakeholders to prepare effective tsunami mitigation plans providing suggestions as to how effectively local evacuation plans can be developed and implemented (e.g. foot or motorbike).

  • The agent-based modelling approach through the selection of the most suitable number of iterations can allow the simulation of different level of preparedness of the population.

The results from this work can inform action plans for the stakeholders and the coastal community to reduce human losses due to future tsunami events. Still, conducting regular evacuation drills is essential in the coastal areas of Padang. The significant number of affected people by the tsunami found in this work is a clear warning for the stakeholders and local communities. Various evacuation routes outside the inundated zone and assessed as safe through stochastic tsunami simulations can be defined and this information can to be disseminated to the coastal communities. Rather than having a limited number of evacuation routes, it is highly suggested to increase the number of evacuation routes for Teluk Bayur and Padang areas. The local government and coastal community can make the areas more accessible and hence, can increase the number of evacuation routes for the coastal community. Tsunami evacuation shelters can also be built in the Teluk Bayur areas specifically at the harbour areas where the horizontal evacuation system may not be sufficient to save people.

Providing rich information of tsunami risk uncertainty resulting from the proposed framework may guide the related stakeholders in coastal community to develop reliable mitigation strategies and effective communication plans to the coastal community. A worst risk scenario is primarily used in Indonesia to inform measures for the coastal communities; such an approach may lead to an over-fear situation. This issue can be solved by an efficient communication of the risk, providing a range of future tsunami risk scenarios, to coastal communities. A further development for policy purposes of this model can be the comparison of agent-based models with different numbers of iterations that can also allow including different level of preparedness of the population leading to a quantification of the benefits in terms of evacuation time. Finally, the proposed framework developed in this study may be implemented in the other tsunami prone areas (e.g. Japan and Mexico) to quantify the uncertainty of the tsunami evacuation-risk metrices (e.g. evacuation time and number of affected people).