GIS-based study of tsunami risk in the Special Region of Yogyakarta (Central Java, Indonesia)

The island of Java is located above a plate boundary within the Indonesian archipelago as part of the volcanic arc resulting from the subduction of the Indo–Australian plate beneath the Eurasian plate. Due to its emplacement, it is exposed to several geohazards, including active volcanoes and earthquakes, and its secondary effects (e.g., landslides, tsunamis). Tsunamis have repeatedly hit Java; for example, the 17th of July 2006 caused by an Mw 7.8 earthquake off the coast of western Java, or the 22nd of December 2018 triggered by the collapse of the flank of Anak Krakatau during its eruption. In most cases, tsunamis have a destructive coastal impact and can cause a significant loss of lives. Therefore, it is of high importance to determine areas situated in tsunami risk zones and estimate plans for risk reduction and prevention. We present an ArcGIS-based method to calculate tsunami risk zones for the Special Region of Yogyakarta (Central Java), including the calculation of hazard zones, vulnerability zones, and the estimation of highly exposed areas. The final risk map reflects the current risk zones in the occurrence of a tsunami and should be recalculated and re-evaluated if the environmental and socio-economic conditions change. The results can be used as a base for tsunami evacuation plans and site planning.


Introduction
Indonesia is an archipelago located in Southeast Asia, in the Pacific ring of fire, the most seismically and volcanically active region in the world. Java Island is located a hundred kilometers north of Java Subduction Zone (figure 1). This area accumulates part of the shortening between the Eurasian and Indo-Australian Plates, which converge at the rate of 5-7 cm/yr ( [1]; [2]). This plate-tectonic configuration builds up a high density of active volcanoes along the island and the occurrence of large numbers of earthquakes and its secondary processes (e.g., landslides, liquefaction, tsunamis).
The Special Region of Yogyakarta is located in the central south of the island of Java. With 1207 inhab./km² it is the fourth highest densely populated province of Indonesia and one of the fastestgrowing provinces regarding social and economic aspects [3]. In the past years, a high urban expansion rate was caused by a rising migration rate and economic growth. Reasons for this include the excellent reputation of Yogyakarta as "The City of Education" and "City of Culture", and the rising number of tourists; Yogyakarta is the second largest tourist destination besides Bali ( [4]; [5]). As a consequence, the hazard and exposition of society to natural hazard rises.
The northern region of the province is located on the Merapi volcano's feet, while the southern region is delimited by the coast occupied by settlements, industries, tourism areas, agricultural fields, and the new Yogyakarta International Airport. The province is frequently hit by volcanic eruptions from the active Merapi volcano and earthquakes related to the subduction regime and shallow faults in the vicinity of the area (e.g., Opak fault). Regarding tsunamis, two events affected the coast of Yogyakarta in the International Conference on Geological Engineering and Geosciences IOP Conf. Series: Earth and Environmental Science 851 (2021) 012007 IOP Publishing doi: 10.1088/1755-1315/851/1/012007 2 historical period. On the 2nd of June of 1994, an Mw 7.9 earthquake occurred in Banyuwangi (East Java) and caused a tsunami with run-up heights over 10 m leaving more than 250 fatalities (e.g. [6]). On the 17th of July 2006, an Mw 7.8 earthquake off the coast of western Java caused a tsunami that affected over 300 km of coastline and caused more than 600 fatalities, with average run-up heights of 5 to 7 m [7]. Subduction zone earthquakes are by far the most typical source of tsunamis. However, there are other sources capable of generating a tsunami, including volcanic eruptions and (submarine) landslides. On the southern coast of Java, no tsunami caused by a submarine landslide has ever been historically documented [8]. Regarding volcanic eruptions, the Krakatau (located to the west of Java, between Java and Sumatra Islands) triggered several tsunamis due to its eruption. In 1883 the run-up height was higher than 20 m in Java, causing more than 35000 fatalities (e.g. [9]). More recently, on the 22nd of December 2018, the flank of Anak Krakatau collapsed during eruption and generated a devastating tsunami on the eastern coast of Java, leaving a death toll of 437 (e.g. [10]).
The previously mentioned local tsunamigenic events and large-scale events such as the 2004 Indian Ocean Tsunami (one of the most disastrous events in recent Indonesian history); together with the recent studies of [11] proposing the likelihood of megathrust events along the Java trench, suggested by seismic gaps (from seismic and geodetic data) in the subduction trench, and the study of [12] that found tsunami deposits evidence at many sites along the southern coastline of Java suggesting a return period of approximately 500 years for this megathrust earthquakes, demonstrate the importance of tsunami risk assessment in this area.
Several worldwide efforts were put into the improvement of tsunami risk assessment in Southern Asia after the 2004 Indian Ocean Tsunami. Rakowsky et al. [13] computed a complete database with tsunami scenarios covering the Sunda Trench as an aid for tsunami inundation studies and the Indonesian tsunami early warning system. Tsunami hazard maps are available for various municipalities of the southern coast of Java including Bantul and Parangtritis in the Special Region of Yogyakarta ([14]; [15]). However, tsunami risk maps for the Special Region of Yogyakarta are not freely available. We present a static approach to model risk maps based on hazard, vulnerability, and exposure. The results can be used as a base for tsunami evacuation plans and site planning.

Seismotectonic setting
The Island of Java is located on the Sunda volcanic Arc, due to the subduction of the Indo-Australian Plate beneath the Eurasian plate. The seismicity of Java is characterized mainly by interplate earthquakes that occur along the subduction slab interface in the Wadati-Benioff zone. The Wadati-Benioff zone slides gently, almost horizontal, beneath a melange wedge for the first 150 km from the Java trench, where the dipping angle increases rather sharply to 45° [16]. The slab prolongates to depths of about 650 km dipping approximately 65° northward [17]. This change in the steepness might cause a northward pushing and stress accumulation in the overriding plate [16]. Earthquakes of intermediatedepth (40-250 km) are generated beneath the Java Island, and deep earthquakes (400-650 km) at the back-arc region located north of the island and beneath the Java Sea. Regarding the historical seismicity, the earthquake magnitudes are likely to be in the range of Mw 7.0-7.5 [18]. These earthquakes generally have high tsunamigenic potential. In some cases, these events have demonstrated slow moment-release and have been defined as 'tsunami' earthquakes, where rupture is large in the weak crustal layers very close to the seafloor [18].

Methods
The term "risk" is described as the probability of harmful consequences or expected losses (deaths, injuries, property, livelihoods, disrupted economic activity, or damaged environment) arising from the interaction between a hazard, vulnerable conditions, and the exposed elements (people, infrastructure, among others [19]. Risk can be expressed in an equation (1) [20] as the sum of 1) hazard, an occurring physical event which can cause damages; 2) vulnerability, the capacity of conditions to resist hazards; and 3) exposure, different features located in hazard-prone areas [19].
We use a multi-criteria GIS-based method (figure 2) using the Geographic Information System Software ArcGIS 10.6.1 to determine the tsunami risk in the Special Region of Yogyakarta. We first compute a map for each component of the risk equation (eq 1), considering different criteria and combining them into a final risk map. We followed two approaches: 1) the hazard and vulnerability map were combined to produce an inherent risk map to determine the areas that would be hit by a tsunami; and 2) we include the exposure factor to define the areas that would be more affected in case of a tsunami to assess disaster management. The input data includes Landsat 8 satellite imagery [21], digital elevation model DEMNAS [22] and SRTM [23], land use data ( [24]; [25]; [26]), and population density [27]. All data are freely available online.

Hazard
To estimate tsunami hazard zones, we determine how far tsunami waves with different wave heights at coastline will inundate inland, the so-called inundation limit. For calculating the tsunami inundation in ArcGIS, the analysis toolbox "TsunamiInundation_10.1-10.2" from [28] was used. For each raster cell of the digital elevation model of the study area, the included tool in the toolbox, "Tsunami Inundation," calculates the loss of wave height per meter with a usage of equation (2), described by [29]. Next, the maximum inundation width is calculated with equation (3). The Manning roughness coefficient (n) indicates the ground resistance to flood flows in channels and flood plains [30]. This coefficient has a high influence on the tsunami inundation formula (2) [31]. It has been proved that for a lower roughness coefficient the tsunami would travel four times further inland than for higher coefficients [32]. The reason for this is that higher roughness involves higher friction, which directly influences the velocity of the water flowing inland. For calculating the roughness for the tsunami inundation, we use the Manning coefficient formula developed by [32] and modified from [30], which includes several roughnesses to represent the floodplain characteristics: Where n is the Manning coefficient, nb is a base value for the flood plain´s natural bare soil surface (geology), n1 is a correction factor for surface irregularities (morphology), n2 is a value for variation in shape and size of the floodplain cross-section (0 for floodplains), n3 is a value for obstruction, n4 is a value for vegetation and flow conditions (land use) and, m is a correction factor for sinuosity of the floodplain which takes a value of 1 for floodplains [30].To estimate the parameter n1 we determine the surface roughness with the Terrain Ruggedness Index (TRI). For parameter nb and n4 we consider the geological and land use characteristics of the area and we assign a Manning value proposed by [30] (e.g. coarse gravel -n= 0.03). It is not practicable to determine the parameter n3 for the study area since it is nearly impossible to distinguish small scale obstructions (debris deposits, exposed roots, logs, isolated boulders, among others) from satellite images. Therefore, we set the value to zero to avoid an overestimation. Then, the single parameters are summed up to a final value n, whereby the parameter nb is taken into account only with 25% to avoid an overestimation of the geology since the land use overprints the geology.
To compute different hazard tsunami scenarios, the tsunami inundation is calculated using formula (2) for different wave heights on the shoreline. The wave heights are divided into five classes and sorted according to their recurrence period into a scale of hazard classes. By comparing the wave heights of historical tsunamis [33], wave heights of 0 -5 meters are more probable than 10 -15 meters. Therefore, the class 0 -5 meters wave height is attributed to the hazard class "Very High," and so forth.

Vulnerability
According to the definition of vulnerability we use an environmental multi-criteria approach based on topographic elevation, slope, and roughness derived from land use following the workflow analysis proposed by [34] to estimate the vulnerability of the study area in the event of a tsunami. We produce a map for each criterion in ArcGIS, and then the criteria are weighted regarding their relevance against each to construct a final vulnerability map. The topographic slope is directly linked to the tsunami inland flow. For gently sloping areas, the vulnerability is higher than for steep sloping areas. A gentle slope represents less resistance for the water to flow; therefore, the tsunami wave's flow velocity and force are not reduced as much [35]. For the slope vulnerability map, we use an older version of SRTM data (3 arc-second) from [23] with a resolution of 90 meters. Newer DEM versions with higher resolutions like DEMNAS [22] are not suitable since vegetation and buildings/structures change the general slope on a very small scale, which leads to a wrong calculation of the topographic slope. We use the "Slope" tool in ArcGIS to compute this parameter. For the vulnerability classification, the classification from [36] were used (also implemented by [34]), where the slope is divided into five groups from low slope gradients associated with a very high vulnerability to high slope gradients associated with a very low vulnerability. Regarding elevation, high topographic regions are less vulnerable to the tsunami inland flow than low topographic areas. The tsunami wave needs a more considerable wave height and energy to reach high topographic areas. The elevation height is based on the DEMNAS DEM. For the elevation vulnerability classification, we use, as proposed in [34], the classification of [37], which classifies the elevation height in five classes, from low elevations associated with a high vulnerability to high elevations associated with a low vulnerability. As an additional criterion to the model of [34], we include the surface roughness and excluded the criteria coastal shape, proximity and direction, since the method is not taking into account a single tsunami point source and the coastal shape is almost straight everywhere. Due to the fact that the roughness influences the vulnerability (areas of low roughness are related to higher vulnerability), the parameter is additionally included here. The roughness coefficient is calculated with the Manning coefficient (as described in section 3.1) and represents its own criterion. The roughness vulnerability classes are set, based on [31], from very high for 0 -0.01 to very low > 0.1.
The three calculated parameters (elevation, slope, and roughness) are merged using a weighted overlaying analysis to generate a final environmental vulnerability map. To estimate the percentage weighting factor for each criterion, we use a method called "Analytical Hierarchy Process" (AHP) introduced by Thomas Saaty in the 1970s [38]. With this method, the criteria are compared, due to their importance, with each other pairwise [39]. Then, the three parameters are merged in ArcGIS using the tool "Weighted Overlay".

Exposure
People or infrastructures situated in a tsunami hazard zone are directly exposed to a potential tsunami and, therefore, at risk. We use various criteria to evaluate the socioeconomic exposure to a tsunami in the study area, which includes the categories: buildings and facilities, tourism areas, areas of critical land use, critical infrastructure and, critical transport. All criteria are integrated with ArcGIS into one map. Then each criterion is given an exposure class ("Very High" to "Very Low") for the day-time and one for the night-time. Since the exposure spatially varies from day-time to night-time (e.g. more public places like the beach site are visited during-day time than during night-time), we generated two maps to depict this fact. This evaluation is based on statistical data from [24], [25], [26]; and field observations from the fieldwork carried out in June and July of 2019.

Results
The resulting tsunami risk map for the first approach integrates hazard and vulnerability. In general, most of the Special Region of Yogyakarta is categorized as a "Low" and "Very Low" while 1.7% of the total area as "High" and "Very High". High-risk areas are located in the southern coastline of the districts Bantul and Kulon Progo, extending up to five kilometers inland following the main rivers. Contrarily, almost the entire coastline of Gunung Kidul is categorized as "Low" tsunami risk area due to its steep cliffs of a high of 40 meters, in exception of small bays ( figure 3).
The tsunami risk map (figure 4 and 5) obtained with the second approach, which incorporates the exposure factor, shows a similar pattern as the inherent risk map. Approximately 134,502 people (3.5%) [27] are located in a "Very High" to "Moderate" risk class area during day-time and approximately 122,974 people (3.2%) [27] during night-time. The areas categorized as "Very High" and "High" are adjacent to the coastline and belong mainly to the socioeconomic categories: "Beach", "Airport" and "Tourism Area". A comparison between the results for day-time and night-time shows that the percentage amounts in the risk class "Very High" decreases from dayto night-time. It is essential to point out that 25.1% of the medical facilities and 17.3% of the schools of Kulon Progo and Bantul districts are located in a categorized "High" and "Moderate" tsunami risk area [27].

Discussion
For the estimation of tsunami risk in the Special Region of Yogyakarta, we made two approaches. The first approach is based on the given natural conditions (slope, elevation, roughness, and proximity from the coastline) where risk areas are calculated considering their inherent terrain condition (what we name here as vulnerability) and location (hazard -inundation). The second approach includes the factor socioeconomic exposure, where risk areas are calculated considering the terrain conditions (vulnerability), the location (hazardinundation), and where and how many critical infrastructures and people are situated in these areas during day-and night-time. As mentioned above, it is important to differ for day-time and night-time because public places like the beach site are more exposed during the day than during the night whereby at night-time the settlements are more exposed. Both approaches consider different perspectives. If only the second one is taken as a final risk map, it could lead to a misunderstanding.
Risk is complex and depends on a number of parameters. It is different understood depending on the perspective. For example, those construction companies interested in knowing where to build up new constructions in a risk-free tsunami area focus on the inherent tsunami risk map, which is only dependent on the site's natural conditions. Alternatively, local authorities and agencies tasked with managing and The same applies to vulnerability. There are several types of vulnerability: physical, economic, social, and environmental. In this study, we only consider the environmental vulnerability and include some of the socioeconomic exposure's factors. We understand our data as a preliminary and fast method to overview the risk area in the Special Region of Yogyakarta. Accurate studies in vulnerability should include other elements such as the construction standards of buildings, socioeconomic status, level of understanding and hazard perception, the warning systems, and preparedness (ability to move away from the flood zone to a safe place), among others. The method we use is strongly dependent on the available input source data. We use open-source data that may not be up-to-date anymore. Therefore, the method is prone to error for example, in calculating the roughness parameter and the exposure factor, which are mainly based on land use data. The results will reflect the potential general tsunami risk zones for the current environmental and socioeconomic situation since we live in a changing world (and more so today).
Many studies focused on numerical modelling on tsunami propagation in the Java arc [13] which is an important parameter for a risk calculation. Since we did not take into account specific tsunami source points in the hazard calculation, the results of our model are not highly accurate. If we would consider specific tsunami source points, the results would increase in accuracy but would also increase the complexity, computing time and loose the general factor. Our general hazard results are still in good agreement with the existent hazard maps for parts of the Special Region of Yogyakarta (Bantul and Parangtritis) ( [14]; [15]). In advantage of our model is a relatively simplistic approach thus a short computing time is needed, that still delivers adequate results.

Conclusion
Tsunami risk zones were calculated for the Special Region of Yogyakarta following a GIS-based workflow. Two approaches were made for risk calculation where the results of the first approach show risk zones based on the environmental terrain conditions, and the results of the second approach include the exposure factor. For the final risk calculation, eight parameters were taken into account: geology, terrain conditions, land use, slope, elevation, tsunami inundation, critical infrastructures, and exposure of people. It is shown for both approaches that most areas of the Special Region of Yogyakarta are located in low or very low tsunami risk areas. Nevertheless, on the southern coastal area, the districts of Kulon Progo and Bantul are situated in high tsunami risk zones susceptible to be directly hit and most severely affected by a tsunami. Since the Special Region of Yogyakarta is prone to tsunami hazards, it is very important to know which areas are situated in risk zones and create risk reduction and preventative measures on this basis. Consideration should also be given to building high tsunami evacuation towers near the coast as well as to expand the early tsunami warning system.
In general, the method is highly dependent on the input source data and their quality. Influencing parameters like the land use must be documented accurately. The results should be taken as approximate and only valid for the present-day due to the changing current environmental conditions including climate change and sea level rise which change the land use and erosion/sedimentation rates. Furthermore, the results show a general risk estimation which is not based on a specific tsunami event.
Summarizing it can be used as a fast efficient and easy method.