Preliminary identification of earthquake triggered multi-hazard and risk in Pleret Sub-District (Yogyakarta, Indonesia)

ABSTRACT Yogyakarta is one of the large cities in Central Java, located on Java Island, Indonesia. The city, and the Pleret sub-district, where the study has taken place, is prone to earthquake hazards, because it is close to several seismically active zones, such as the Sunda Megathrust and the active fault known as the Opak Fault. Since a devastating earthquake of 2006, the population of the Pleret sub-district has increased significantly. Thus, the housing demand has increased, and so is the pace of low-cost housing that does not meet earthquake-safety requirements, and furthermore are often located on unstable slopes. The local alluvial material covering a jigsaw of unstable blocks and complex slope is conditions that can amplify the negative impacts of earthquakes. Within this context, this study is aiming to assess the multi-hazards and risks of earthquakes and related secondary hazards such as ground liquefaction, and coseismic landslides. To achieve this, we used geographic information systems and remote sensing methods supplemented with outcrop study and existing seismic data to derive shear-strain parameters. The results have revealed the presence of numerous uncharted active faults with movements visible from imagery and outcrops. show that the middle part of the study area has a complex geological structure, indicated by many unchartered faults in the outcrops. Using this newly mapped blocks combined with shear strain data, we reassessed the collapse probability of buildings that reach level >0.75 near the Opak River, in central Pleret sub-district. Classifying the buildings and from population distribution, we could determine that the highest risk was during nighttime as the buildings susceptible to fall are predominantly housing buildings. The secondary hazards follow a slightly different distribution with a concentration of risks in the West.


Introduction
Cascading hazards have received relatively little attention in disaster and risk studies, as the majority of hazard and disaster risk studies have concentrated on one single hazard. The reality of scientific funding and projects often leads researchers to pay less attention to the possibility that one incident can cascade to other secondary hazards. Consequently, the interactions between hazards are still relatively ignored (Budimir, Atkinson, and Lewis 2014), especially because as scientists we are often trained in one discipline and it is still problematic to work across boundaries. As a particular area may be exposed to more than one type of hazard, even if each hazard can lead to a given type of disaster with different magnitudes (Westen 2011), each single incident can trigger secondary hazards that cannot be encapsulated by the magnitude approach from the original hazard (Ren and Liu 2013). An example of this such disaster chain-reaction occurred in Beichuan, China, due to the 2008 Wenchuan earthquake. The Wenchuan earthquake generated a complex disaster chain that caused important damage to towns in the Beichuan county. The earthquake had a shallow focal depth of about 19 km, and the magnitude of the earthquake was 7.9-8.0 Ms. It caused ~90% of the buildings in the area to collapse. They had also triggered coseismic landslides upstream the township, generating in turn a dammed lake, which collapsed with the heavy rainfall a month after the earthquake. The water from the Tangjiashan Lake mixed with the rubbles of the earthquake and washed away the ruins of the collapsed buildings. Fortunately, during the event, the Beichuan county town was off-limit due to the expected cascading effects triggered by the earthquake.
Another similar situation occurred in Christchurch, New Zealand between September 2010 and February 2011, Christchurch -the second largest city in New Zealand -experienced two large series of earthquakes. The first earthquake (6.3 Ms) took the lives of 185 people, making this event the second deadliest disaster that ever occurred in New Zealand. Amidst numerous aftershocks, a second earthquake (7.3 Ms), just underneath the city, caused significant damage in the central city of Christchurch. Afterward, in March 2014, Christchurch experienced more flooding due to the impact of the Canterbury Earthquake and the general lowering of the ground level (Gomez et al. 2019). Some researchers concluded that this extraordinary flood occurred not just due to the heavy and prolonged rainfall but due to ground deformation, liquefaction, subsidence, narrowing of channels, and uplifting of river beds after the earthquake.
Coastal plains are therefore particularly at risk, especially for earthquake-prone island coastal plains of New Zealand, Japan or Indonesia. In Indonesia, the study location is situated approximately 300 km north of the Java Mega thrust. The study area, in Central Java (the west flank of Baturagung Escarpment) is particularly prone to earthquakes. Having a total subduction segment of 840 km, the Java Mega thrust can potentially generate earthquakes with maximum recordedmagnitudes up to 8.1 Mw (Figure 1) (Irsyam et al. 2012). Link to this activity, one of the main active fault of Java, the Opak Fault, passes through the research area along the west flank of the Baturagung Escarpment. According to (PUSGEN 2017), which is the National Center of the Earthquake studies, with a total fault length of approximately 45 km and a slip rate of about 2.4 mm per year, the Opak fault can potentially generate shallow earthquakes with a maximum magnitude of 6.2 Mw. According to the same agency, such an earthquake would have a devastating impact on the study area and the southeast part of Yogyakarta City.
Administratively, the west flank of Baturagung consists of several sub-districts, including Kretek, Pundong, Imogiri, Dlingo, Pleret, and Piyungan. As a pilot study, the present work has focused on the Pleret sub-district, Bantul Regency, of Yogyakarta Province. Pleret sub-district is located approximately 10 km southeast of Yogyakarta Central City. Located near to subduction zone and over the active Opak Fault makes this area very prone to earthquake hazards. Furthermore, high seismicity is expected to be amplified due to unconsolidated river and marine sediments, which could lead to liquefaction in the center lowland, while to the west, amplified movement is expected due to the unstable slope of the Baturagung Escarpment.
Probabilistic Seismic Hazard Analysis (PSHA) (Ashadi and Kaka 2015), Deterministic Seismic Hazard Analysis (DSHA) (Chopra et al. 2012), liquefaction (Civico et al. 2015), ground motion, landslide (Hadmoko et al. 2010), and tsunami studies have been developed worldwide. However, earthquakes can trigger other related hazards that increase the impact on society (Marano et al. 2010) and studies of multi-hazards and risks in this area are still limited. Thus, the study of earthquake triggers and other related secondary hazards is needed in Pleret Sub-District, Southeast Yogyakarta Province.
The primary objective of this study was to construct a multi-hazard risk assessment in Pleret Sub-District using remote sensing and a Geographic Information System (GIS). The primary objective can be divided into several sub-studies with the objectives of (1) identifying the potential area of coseismic landslide, (2) conducting an outcrop study in order to better understand the fault configuration, (3) identifying the liquefaction zonation, (4) assessing the vulnerability level of some elements at risk, and (5) describing the multihazards and risks in the study area.

Overview of the study area
Pleret sub-district is located 10 km southeast part of Yogyakarta City, Indonesia. The study area is the part of Bantul Graben. This graben was formed by two convergent normal faults, which generated two horst zones in the east and the west and a graben zone in the middle. Two normal faults, Progo and Opak Faults, formed the Bantul Graben. Based on the difference in gravity, Progo fault is located near the Progo River in the west, and Opak Fault is located in the border area of the Southern Mountain (Baturagung Escarpment) to the east (Nurwihastuti et al. 2014). The study area is located in the transition zone between a flat area and the escarpment zone in the Eastern part of Bantul Graben ( Figure  2). Saputra et al. (2018) stated that Pleret sub-district is located in the Eastern Horts of Bantul Graben, which has experienced step faults due to the complexity of specific geological structure.
Geographically, the study area can be divided into three major zones: east, middle, and west. The west and the middle part zones are dominated by extensive flat areas, while the Eastern zone is dominated by undulating and mountainous areas. The extensive flat area mainly consists of Young Volcanic Deposits of Merapi Volcano (Qmi) and few of them consist of alluvium (from the denudational process in the Eastern mountainous areas), which is located in the foothills, near the undulating areas. The Eastern part consists of mainly tertiary deposits of Semilir Formation (Tmse). Only a few areas, including the summit of the Baturagung Escarpment, belong to the Nglanggran Formation (Tmn). Tmse formed approximately 27.82-23.03 million years ago. Tmse consists of mainly interbedded layers of breccia pumice, tuffbreccia, dacite tuff, andesite tuff, and tuffaceous clay, whereas the Tmn was deposited in early Miocene approximately 23.03 to 11.608 million years ago. Tmn was deposited in parallel on the top of Tmse. Tmn is mainly distributed in the summit of Baturagung Escarpment in the Eastern part of the study area. Both formations were generated from the eruption of a nearby ancient volcano (Bronto et al. 2008;Winarti and Hartono 2015;Pandita, Sukartono, and Isjudarto 2016). Bronto et al. (2018) found that both Tmse and Tmn can be classified as pyroclastic density flow from the ancient volcano that might be located in the east part of the study area.
In terms of geomorphology, three major groups of landforms can be identified in the study area: structural, fluvial, and denudation. The structural landforms are indicated by the existence of the Baturagung Escarpment in the east part. The intensive denudation process occurred in the middle to upper slope of the escarpment. Some triangle facets formed due to the strong erosion processes in the middle to upper slope of the escarpment. The fluvial process occurs in the middle part of the study area along the Opak River and in the narrow flat area near the undulating area and Baturagung Escarpment. The fluvial process along the Opak River and narrow flat area produced an extensive alluvial plain and a colluvial plain, respectively. The alluvial plain consists of Qmi, whereas the colluvium consists of Qa (denudation material from Semilir and Nglanggran Formation). Figure 3 depicts the geomorphological aspect of the study area.

Scope of the analysis
Situated in a complex geological and geomorphological zone, the study area is prone to several hazards, with frequent and strong seismic activity. Pleret Sub-District is located approximately 250 km north of the active Sunda Megathrust and intersected by an active inland fault namely Opak Fault. The middle part of the study area is prone to soil amplification because this area is dominated by the dense quaternary material from the Merapi Volcano. The abundance of shallow groundwater in this area can also potentially generate liquefaction. Furthermore, the mountainous area in the east part of the study area is also prone to slope stability issues and coseismic landslides due to the type of lithology and intensive erosion. Based on these considerations, the multi-risk study several types of hazards such as earthquake, soil amplification, liquefaction, and coseismic landslide aspects.

Data
The study mainly used the historic earthquake data from the United States Geological Survey (USGS) between 1900 and 2019 to support the probabilistic seismic hazard analysis and liquefaction analysis. The fieldwork data of the outcrop study (Saputra et al. 2018) and seismic vulnerability index (kg) from micro tremor measurement (Daryono 2011) were used in the earthquake and liquefaction analysis, respectively. The direct measurement data of groundwater were obtained from household well measurements. Some boreholes and geoelectric data were used to support the analysis of groundwater condition in the middle part of the study area. The rainfall data from surrounding weather stations, a 1:25,000 topographic map, and detailed geology map based on remote sensing analysis were used to determine the coseismic landslide-prone area (Saputra et al. 2015(Saputra et al. , 2016. Furthermore, detailed land use data and building damage were used to analyze the vulnerability level. The land use data were generated using the visual interpretation of the latest QuickBird image (Saputra et al. 2017). The data for the building damage caused by the 2006 Yogyakarta earthquake were obtained from previous research conducted by Kerle (2010). The list of data used in this study is provided in Table 1.

Peak ground acceleration (PGA)
As much as 4,593 earthquake data from USGS were used to obtain the PGA of 320 observation points which were distributed entire study area. Kanai (1966) (Douglas 2019) attenuation was used to generate the peak ground acceleration. Kanai attenuation was used because this formula considers the fundamental period of the site, which is closely related to the local site effect that can amplify earthquake shaking. We calculated the Peak Ground Acceleration (PGA) in the study area as follows: whereα is peak ground acceleration, α1 is the first constant value (5), α2 is the second constant value (0.61), T is the fundamental period of the site, M is the earthquake magnitude, R is the hypocenter (km), and P and Q are the values from equations 2 and 3, respectively.
Thus if the equations 2 and 3 are substituted in equation 1, the Peak Ground Acceleration (PGA) can be calculated using equation 4 below.

The outcrop study
The outcrop study was conducted to get better understanding of geological condition in study area as only geological map of Yogyakarta (scale 1:100,000) was available in study area. The outcrop study was divided into three main stage: pre-fieldwork analysis, fieldwork activities, and post-fieldwork analysis. The prefieldwork analysis included geological data (lithological and geological structure) extraction from mainly Yogyakarta Geology map scale 1:100,000 and Landsat 8 interpretation (Saputra et al. 2018). The visual interpretation of QuickBird imagery was conducted in the pre-fieldwork analysis of the outcrop study. The primary purpose of the interpretation was to identify the location of the outcrop and determine the location for outcrop observation. The aims of the fieldwork activities were to characterize the outcrop (identify the lithofacies and qualitative grain size of each rock layer), identify the micro fault, and to record the threedimensional (3D) surface model of the outcrop using the structure from motion technique. The brief workflow of the outcrop study stage is provided in Figure 4.

Coseismic landslide assessment
The method proposed by Mora and Vahrson (1999) was adopted to generate the coseismic landslide susceptibility (Saputra et al. 2016). There were two main parameters, site characteristic factors and trigger factors, that could cause coseismic landslide occurrence. The site characteristics consist of physical parameters which are related to slope stability analysis such as relief, geology, and soil humidity. Meanwhile the trigger factors consist of seismic and rainfall intensity. As shown in Table 1, other spatial data, such as relief and lithology, were derived from contour and geology maps, respectively. The contour map resolution was about 12.5 m. The geology map scale was 1:100,000. Saputra et al. (2016) interpreted some remote sensing data and completed fieldwork observations to increase the scale of the geology map and provided more detailed geology and lithology information. The other parameters, such as natural humidity of soil and rainfall intensity, were derived from monthly average rainfall and annual rainfall intensity, respectively. The detailed work flow of the coseismic landslide assessment that was conducted by Saputra et al. (2016) is provided in Figure 5.

The liquefaction assessment
Liquefaction is a secondary earthquake hazard that often causes the worst damage to cities around the world. In 1999, a 7.6 Mw earthquake occurred in Chi-Chi and Taiwan, China. This earthquake generated massive liquefaction throughout Taiwan, China and especially in Nantou, Wufeng, and Yuanlin prefectures (Chu et al. 2004). In New Zealand, two strong sequential earthquakes in September 2010 and December 2011 in Christchurch caused massive liquefaction across Christchurch (Morgenroth, Hughes, and Curbinovski 2016). On 28 September 2018, a strong earthquake struck Palu City, Central Sulawesi, Indonesia. The Palu-Koro strike-slip fault generated inland earthquake with a magnitude of 7.4 on the Richter scale at an earthquake depth of 10 km. This earthquake mainly affected Palu Central City. A 1.5-m tsunami occurred in Palu and Donggala coastal areas. Massive liquefaction also occurred in some parts of Palu, causing severe damage.
Liquefaction is closely related to geological and geomorphological characteristics, as liquefaction often reoccurs in the same locations (Youd and Perkins 1987;Yasuda and Tohno 1988;Tatsuoka et al. 1980 Yogyakarta, Indonesia earthquake, the liquefaction occurred in a flat area, with shallow ground water table, not too far from the earthquake source, and has a high seismic vulnerability index (Kg). Thus, an area with these same characteristics are highly susceptible to liquefaction. Using a similar approach (Yasuda and Tohno 1988;Tatsuoka et al. 1980), the liquefaction assessment was conducted in the study area. We also analyzed the liquefaction from the share-wave velocity and seismic vulnerability index point of view. By knowing the share-wave velocity and seismic index vulnerability, the sites exposed to liquefaction risk can be identified (Daryono 2011;Yasuda and Tohno 1988). The higher the value of ground sharestrain, the more easily the ground deformation occurs, leading to a crack, liquefaction, and coseismic landslide (Huang and Tseng 2002;Ishihara 1982). The relationship of the shear-strain value and the potential liquefaction is outlined in Table 2. Based on this concept, we combined Probabilistic Seismic Hazard Assessment (PSHA) and the Ishihara concept to analyze the potential liquefaction area in the study area. Thus, the first step of this analysis was to determine the Peak Ground Acceleration (PGA) based on the historical earthquake data. The PGA value was used to assess the ground shear strain based on the earthquake vulnerability index (Kg). The ground shear strain value in the study area was calculated using following equation (Daryono 2011;Ishihara 1982): The complete workflow used to assess the potential liquefaction area in the study area is provided in Figure 6.

Multi-vulnerability assessment
The multi-vulnerability assessment used in this study followed the previous research conducted by Saputra et al. (2017). We applied logistic regression analysis to predict the damage level of the residential buildings based on the 27 May 2006 Yogyakarta earthquake damage data. The main data used in that research were earthquake damage data, QuickBird images, geologic maps, and building footprint data. The earthquake damage data were obtained from Kerle (2010), which included building damage data (low, medium, and collapsed) in impacted area (Jetis, Pleret, Imogiri, and Bantul Sub-Districts). The QuickBird imagery was used to obtain the land use data in more detail (the year 2012, scale 1:25,000) based on the modified Anderson system 2002. Saputra et al. (2017) used the geology map of Yogyakarta to extract additional data such as the lithology, type of material, and distance from the epicenter of Yogyakarta earthquake (27 May 2006). In general, there were two main elements at risk used in the analysis: population and building collapse probability. The population data were obtained from the local statistics agency (Badan Pusat Statistik (BPS) or Statistics Indonesia), and the building collapse probability was generated from the logistic regression analysis. We followed four steps to conduct the multivulnerability assessment. The first step was to extract the land use data via the visual interpretation of QuickBird imagery. The second step was to conduct the probability of building collapse. The third step involved creating a population distribution model and the last step involved combining the building collapse probability and population distribution into multi-vulnerability analysis. The steps of the multivulnerability analysis are provided in Figure 7.
We used 15 parameters to calculate the probability of building collapse based on the 2006 earthquake damage data. The parameters consist of dependent variable (Y), which refers to the building damage, and independent X variables (X1-X15) ( Table 3). We applied logistic regression based on the binary model (0 and 1 for No or Yes, respectively). The input data (the spatial building damage data) were converted into the binary system. For instance, the building ID 505 has the characteristics of wood structure, asbestos or zinc roof material, located within 8 km of the earthquake epicenter, has a Semilir geological condition, and experienced moderate damage, was converted into binary code as shown in Table 4.

Multi-hazard and risk assessment
The approach to generate the multi-hazard and risk assessment followed the general concept of risk assessment and reduction. Risk is the function of hazard and vulnerability, which can be expressed as follows:  Figure 6. The potential liquefaction assessment.
where R is risk, H is hazard, and V is vulnerability. For the multi-hazard and risk scheme of this study, we referred to Alkema et al. (2018) who considered hazards as the combination of ground motion, soil amplification, liquefaction, and coseismic landslide. In terms of vulnerability, we focused only on building block (building collapse probability) and the population (the population distribution). The multi-hazard  and risk assessment of earthquakes and other related secondary hazards applied in this study are provided in Figure 8. Thus, based on Figure 8, equation 5 was modified to accommodate the multi-hazard aspect in this study. Equation 6 shows how the multi-hazard and risk were generated in this study: where P is PGA value, SA is soil amplification, L is liquefaction, CL is a coseismic landslide, BV is building vulnerability, and PV is population vulnerability.

Peak ground acceleration (PGA) results
Based on the Kanai attenuation, we used the historical data of over 3,481 earthquakes with a magnitude greater than 5 on the Richter Scale, which occurred between 1900 and 2017, to determine that the peak ground acceleration of study area ranges from 531.04 to 967.66 cm/s 2 . Theoretically, the higher the PGA value in a particular area, the higher the degree of damage probability when an earthquake occurs (Yamazaki and Matsuoka 2018; Walter et al. 2008). Based on this attenuation, we found higher PGA values in the middle part of study area, extending to northeast and southeast along the Opak River, which is closely associated with the location of the Opak Fault. Lower PGA values were found in the northwest part in the Yogyakarta City direction. This distribution of PGA values is inversely related to the dominant period of the soil. The direct microtremor measurements that were recorded by Daryono (2011) explained that the highest (0.84) predominant period of soil in the study area is located in the northwest area. The predominant soil period and PGA are provided in Figure 9(a,b), respectively.

Results of outcrop study
Based on the outcrop study that was conducted by Saputra et al. (2018), the study area has a complex geological structure condition. We found significant evidence of fault displacement including the great normal fault in this area -the Opak Fault. Saputra et al. (2018) found the same rock of the Nglanggran Formation (Tmn is described as a volcanic breccia, lava flow containing breccia agglomerate rock and tuff). The main location of the Nglanggran formation is at the summit of Baturagung Escarpment. However, based on the field observation, the authors found an isolated hill of Nglanggran Formation located in the center of study, separated by approximately 4.24 km from the main location of the Nglanggran Formation. Saputra et al. (2018) found that the middle parts of the study area are more vulnerable to ground amplification. At least 30 fault displacements were found in the middle part of the study area, with a maximum displacement of 2.39 m. Most of the faults are typically normal faults, and only a few of them are strike-slip faults. The direction of most of the micro faults is similar to the direction of the Opak Fault. The outcrop study conducted by Saputra et al. (2018) revealed that the Segoroyoso Village, Srumbung Sub Village, the middle part of Bawuran Village, the middle part of Pleret Village, and the middle part of Wonolelo Village are vulnerable due to the complex geological structure and ground amplification. The map of the fault evidence derived from the outcrop study and an example of maximum fault displacement in Srumbung Sub-Village are provided in Figures 10 and 11, respectively.

The groundwater condition
Based on the integrated direct measurement of household wells during the 2015 rainy season, some geo- Table 4. An example of the binary model for the variables. Building ID 505 = wood structure, asbestos roof, within 8 km, and Semilir Formation. electricity observations, especially fluvial plains and local people drilling activities in South Segoroyoso, the study area is dominated by shallow groundwater areas (less than 10 m). This abundance of groundwater is concentrated in the middle part of the study area (Wonokromo, Pleret, middle Segoroyoso, middle Bawuran, and middle Wonolelo). Geomorphologically, this area is dominated by fluvial and colluvial plains, which consists of Alluvium from the Young Merapi Volcanic Deposit (Qmi) and Colluvium from the denudational material of Semilir and Nglangran mountainous area (Qa), respectively ( Figure 12). In line with the household wells observations, the geo-electricity observation showed that the study area has abundant groundwater. Shallow groundwater can be found at depths of 2.5-15 m in the colluvial plain in the narrow plain in the east part of the study area. From the borehole data, we determined a similar variation in groundwater depth. In alluvial and colluvial plains, shallow groundwater can be found at depths of 1-15 m. This shallow aquifer is known as a confined aquifer. Deep aquifers can be found at depths of approximately 80-110 m.

The liquefaction
Based on the ground-shear strain and the potential liquefaction by using Ishihara method, we determined that the areas more prone to liquefaction are the middle part of study area, especially in Kerto, Keputren, and Kanggotan Villages ( Figure 13). This area has the highest ground-  shear strain of 9,221.43 × 10 −6 , found on the alluvial plain of young Merapi volcanic deposits with a sediment thickness of 130 m and groundwater depth of around 2.31 m.

Coseismic landslide
Based on previous studies (Saputra et al. 2015(Saputra et al. , 2016, the study area could be classified into four Coseismic Landslide Susceptibility Levels (CLSLs): negligible, low, moderate, and medium. The negligible zone is the safest area in the study area. This zone is mainly distributed in the west to the middle part of the study area. This zone is located in a flat to gentle slope area and in the alluvial plain zone, colluvial plain, and natural levee of the Opak River. The low CLSL zone is associated with the narrow plain that is located in the border area between the flat and mountainous area in the east part study area. The moderate zones are mainly located in the middle slope of Baturagung Escarpment, which features weakly to strongly eroded denudational hills of Semilir Formation. The medium zone -the most unstable areas -are located along the upper slope of Baturagung Escarpment, which consists of Semilir Formation. This result is in line with that reported by Samodra et al. (2016), who found that the middle slope a greater probability of rock fall occurrence in the Sewu Mountainous Area, which is in the west part of Yogyakarta. The complete coseismic landslide susceptibility map is shown in Figure 14.

Multi-vulnerability results
The first element at risk used in this study was residential buildings. The vulnerability of residential buildings was determined from the logistic regression. There were two datasets used in this study. The first dataset was the building damage data due to the Yogyakarta earthquake. These data include some  information about attributes such as location (x and y), the owner, the building structure, the roof material, geology, the distance from the epicenter, and the level of damage. These data were obtained from field observations right after the disaster occurrence. The second dataset was the building footprint data of existing buildings after rehabilitation and reconstruction process. These data were generated from a visual interpretation process. These data included the same attributes such as location (x and y), the type of buildings structure, roof materials, geology, and the distance from the epicenter.
The other examined element at risk was population density. We used the dasymetric technique to determine the population density per unit land use. Several scenarios, such as the population of each  land use unit in daytime and nighttime both on weekdays and holidays, were applied to determine the real population condition. The land use units were generated from visual interpretation based on the Anderson classification. The vulnerability results of each element at risk are explained below. The results of multi-vulnerability are provided in Figure 20.

Land use based on visual interpretation of Quick bird imagery
Based on the modified Anderson system, we found 29 land use units in the study area. This land use unit covers the classification level III, which is suitable for interpretation of QuickBird or medium-altitude data captured between 3100 and 12,400 m or 1:20,000 to 1:80,000 map scales. Based on this classification, we classified the study area into 29 land use units: abandoned mining sites, agricultural wetlands, canal, cemetery, cemetery on wetland, commercial strip development, educational institution, government centers, harvested cropland, health institution, inactive cropland, light industrial, not built up, open areas, other agricultural, other institutional (mosque), pastureland, poultry farm, residential high density, residential low density, residential medium density, rural single unit, road, shrub land, specialty farm, stone quarries, stream, traditional market, and wetlands.
The distribution and the total area of the land use units are provided in Figure 15.

Building unit collapse probability
There were three main results based on the logistic regression: the model fitness results, pseudo R 2 , and estimated parameter values. The model fitness shows that the models (binary coding) for both independent and dependent variables were significant, that the models fit, and can be used for further analysis (Table 5). The pseudo R 2 result provides information about how far the model can explain the results. Based on the logistic regression, the model can explain at least 33.20%, which was indicated by the Nagelkerke value (Table 6). This means 66.80% of the dependent variables cannot be explained from this model. Thus, some additional parameters need to be added in future research.
The estimated parameter values (Table 7) show that the threshold values of the damage categories are 1.529 and 2.426 for damage categories 1 and 2, respectively. This means the predicted response value (Y*i) were categorized as follows: (1) Damage category 1 (low damage) if Y*i ≤ 1.529 (2) Damage category 2 (moderate damage) if 1.529 < Y*i < 2.246 (3) Damage category 3 (high damage or collapsed) if Y*i ≥ 2.246 Figure 15. The land use map and its total area.
The Y*i value can be calculated from the formula that was derived from the value of each parameter in Table 7. The formula followed the general regression equation. Therefore, from Table  7, the general equation of damage category is expressed as: where the Y*i is the prediction value of damage category, X1 is wood structure, X2 is unreinforced masonry, X3 is reinforced masonry, X4 is asbestos and zinc, X5 is cement tile, X6 is clay tile roof, X7 is concrete slap roof, X8 is distance within 8 km of epicenter, X9 is distance between 8.1 and 10 km of epicenter; X10 is distance between 10.1 and 15 km of epicenter, X11 is distance greater than 15 km from epicenter, X12 is Semilir formation (Tmse), X13 is alluvium (Qa), X14 is young volcanic deposits of Merapi Volcano (Qmi), and X15 is Nglanggran Formation (Tmn). Based on the visual interpretation of 17,512 buildings, there were only 33 combinations of house characteristics. Each combination has specific building attributes and damage category. For instance, the combination number 5 (Table 8)  This means that combination of building attribute number 5 has damage category 3 (high damage or collapsed). The results of damage category calculation using equation 7 of all building unit entire study area can be seen in Figure 16.
The building damage probability of each building unit was obtained by applying the damage category resulted from equation 7 into equations 8, 9, and 10 below (Agresti 1990;Hosmer and Lemeshow 2000).

Probability of no damage Y1
Probability of high damage or collapsed Y3 ð Þ Thus, for example, combination number 5 (wood structure, clay tile roof, distance more than 15 km from the 2006 epicenter, and young volcanic deposit of Merapi Volcano) was categorized to damage category 3 or collapsed. By using equations 8,9, and 10, the building damage probability can be determined. For Based on the calculation above, the building combination number 5 has a predicted damage category of 3 (collapsed) with a probability of high damage or collapsed of 0.655. The safest building type in the study area is the building with the attributes of reinforced masonry structure, asbestos or zinc roof material, located between 10.1 and12 km from the earthquake source, and located above the Semilir Formation. This buildings unit have the probability of collapse only 0.07. The most vulnerable buildings are those with a combination of reinforced masonry structure, clay tile roof material, located between 8.1 and 10 km from the epicenter, and located above the young volcanic deposits of the Merapi Volcano. The probability of collapse of this type reached 0.84.
Building block probability of collapse was determined from converting the probability of collapse of a building unit to the block or land use scale. The average value was applied to convert the collapsed probability of building unit to the building block. The illustration of how to convert the collapse probability of building units to building block is provided in Figure 17. Based on the results, residential blocks located in the western part of the Opak Fault are rated as high probability of collapse between 0.60 and 0.73. The middle part of a study area, to the left and right of the Opak River, scored the highest values (very high) of building collapse probability, ranging from 0.74 to 0.84. The Eastern part tended to have very low to moderate probability of collapse (0.08-0.59) because this area is mountainous areas consisting of compact lithological characteristics. The distribution of building collapsed probability in the study area are provided in Figure 18.

Population distribution model
The distribution of the population in several scenarios (day time, night time, both on weekends and weekdays) was defined based on dasymetric analysis through the simple percentage of the occupation of the local people. Based on the data, the people of Pleret Sub-District are predominantly casual workers, students, unemployed, entrepreneurs, and farmers, at 17.92%, 20.51%, 18.86%, 20.16%, and 23.99% of the population, respectively. Thus, based on the dasymetric analysis, the population distribution on each block of land use was estimated. For instance, during work hours, the population in Wonokromo Village, the densest village in the study area, tends to distribute into three major land use units: settlements, schools, and commercial areas. The population can be estimated as 27.7% staying in commercial areas, 16.65% in schools, and 16.47% in settlement areas. The rest of the population distributes into other land use units. Similar to this condition, the population in other villages, such as Pleret, Segoroyoso, Bawuran, and Wonolelo, mainly distribute into commercial areas, school, settlements, and agricultural areas. The population distribution in all villages in the study areas showed no significant difference among the scenarios used. The population distribution under several scenarios is provided in Figures 19 and 20.

Multi-risk results
The multi-hazard zonation was generated by applying the summation function of all hazard parameters (PGA, the proximity of the faults, liquefaction, and coseismic landslide). The summation function was determined based on the consideration of the selected hazards being secondary earthquake hazards that can amplify the earthquake impact. Based on the analysis, we found that the index of multi-hazard in the study area is between 0.46 and 0.85, with minimum and maximum values of 0 and 1, respectively. Based on this analysis, the study area is dominated by the multi-hazard value 0.61, which can be classified as a moderate multi-hazard zone. Most of this zone is located in the flat area in the middle part of the study area. These areas are located near Opak Fault, have a complex geological structure, and soft surface sediment of alluvium or colluvium. In terms of PGA, this zone has PGA values around Figure 17. The example of converting the building unit collapse probability to block unit. 676.59-822.12 gal and low ground shear strain, which means vulnerable to ground motion and liquefaction occurrence. The total area of this zone is 39.12% of the study area ( Figure 21).
The multi-risk value was obtained by adding the multi-temporal information of the population distribution. The results show that both in daytime or night time on weekdays, the middle part has the highest multi-risk value because this area has a higher multihazard index value, high population density, and most of the buildings inside the block have a higher probability of collapse. The holidays scenarios produced different results. Based on the holiday scenarios, during both daytime and night time, the population tend to distribute within the residential blocks. For example, Segoroyoso Village has a slightly higher multi-hazard risk assessment value at night. The complete multi-risk value picture is provided in Figure 22.
The multi-risk models produce slightly different results between weekday and holiday multi-risk index and produce distinctive results for daytime and night time. On daytime, both on weekdays and holidays, the high multi-risk index is mainly distributed in the southern part of Wonokromo village, in the middle part of Pleret village, and some areas of Wonolelo village. During the night time both on weekdays and  holidays, high multi-risk index values are mainly distributed in most of the residential areas, as almost all people stay in residential area at night both on holidays and weekdays. This model shows (Figure 22, left picture) that the south part of Wonokromo village, the middle part of Pleret village, the middle part of Bawuran, Segoroyoso, and Wonolelo villages scored high multi-risk index values.  Multi-hazards zonation in study area. Hazards information: PGA was derived based on Kanai attenuation by using USGS earthquake historical data . Liquefaction was generated based on the ground shear strain value of USGS earthquake historical data . Soil amplification was assessed based on the proximity to active faults which was generated from outcrop study. Coseismic landslide resulted from the coseismic landslide model grade 2 proposed by Mora and Vahrson (1999) and the multi-hazard index was generated from raster overlay function with benefits standardization method.

Conclusions and recommendations
Based on the results above, the multi-risk index of the study area for earthquakes and other secondary hazards differs, being influenced by PGA, ground motion, liquefaction, coseismic landslide, population distribution, and building collapse probability. Based on the CLSL analysis, the study area includes negligible, low, moderate, and medium zones. The negligible and low CLSL areas are mainly distributed in the extensive flat region in the west part of the study area. The more vulnerable areas are mainly distributed in the mountainous area in the middle part of the study area. This study also shows that the middle slope of Baturagung Escarpment is more prone to coseismic landslide occurrence. These findings are in line with the landslide study that was conducted in the western part of Yogyakarta (Walter et al. 2008).
The study area has a complicated geological structure. Based on the lineament analysis, we found various lineament features of valleys, ridges, river, sudden changes of river direction, ridge lines, scarp face, and straight drainage segments. All these lineaments are associated with the existence of faults (Samodra et al. 2016;Gannouni and Gabtni 2015). Based on the earthquake fracture displacements recorded on the outcrop surface, the middle part of the study area, around the north to the middle part of Segoroyoso and the middle part of Pleret and Bawuran, has a complex fault configuration as indicated by the dense area of faults and highest displacement. Additionally, based on the lithostratigraphic analysis, the outcrop deposits are ignimbrite type that originated from the ancient volcano (tertiary) in the east part of study area.
Based on the shear strain analysis that was derived from the PGA analysis, the study area, and especially the extensive flat area in from middle to west parts, is dominated by shallow groundwater up to 10 m deep with high ground shear-strain values. This means this area is vulnerable to liquefaction occurrence.
In term of vulnerability analysis, we concluded that, statistically, buildings with reinforced masonry structure (RM), clay tile roof material, built above Young Volcanic Deposits of Merapi Volcano (Qmi), and 8.1-10 km from the epicenter of the 2006 earthquake have a higher probability of collapse. In terms of population distribution, there is no significant difference in population density on each land use unit between holidays and weekdays. However, we found that the population density on each land use unit is different between daytime and night time. Thus, population density affects the multi-risk index, which demonstrated a different pattern in daytime and night time.
The multi-hazard risk results show that residential houses that are located in the center part of study area, including the center of Wonokromo, Pleret, Bawuran, and Wonolelo, and the north part of Segoroyoso Villages, have a high multi-hazard and risk index (>0.60) both during daytime and night time on both weekdays and holidays. This means this area is more prone to earthquake hazards and other related secondary hazards.
This study shows the general characteristic of multi-hazard and risk in study area. This study is a preliminary assessment before conducting more detailed investigation. This study has some limitations and we have recommendations for future research. Longer-term earthquake occurrence data are needed to accommodate the longer reoccurrence period of larger earthquake. The amplification due to the micro fault was only based on the proximity analysis; thus, further detailed investigations are required to calculate the size the amplification for complex geological structures when earthquakes occur from surrounding sources. A good inventory of landslide and coseismic landslide is needed to improve the coseismic landslide hazard result. Additionally, several sitespecific observations need to be recorded to support the coseismic landslide data. A deterministic seismic hazard analysis model needs to be created to obtain a better understanding of the type of fault and mechanism that might trigger coseismic landslides when the earthquakes occur. Some limitation of liquefaction analysis was also found in this study such as this study was still used Ishihara model to approach the liquefaction which is only suitable for preliminary assessment. We need to conduct some Standard Penetration Tests (SPT) and cone Penetration Tests (CPT) for verification and for detail site investigation. Pore-water pressure, seasonal effects on groundwater, and storm and heavy rainfall need to be considered as parameters in future research. The statistical analysis of building collapse vulnerability needs to be improved by adding more parameters as independent variables, such as soil type, the shape of buildings, and the age of buildings. In terms of multi-hazard and risk analysis, the cascading effect of earthquakes and other related secondary hazards need to be examined to better explain multi-hazard and risk aspects. The effects of earthquakes also need to be considered on the karst landform located in the southeast part of the study area. Strong earthquakes can cause sinkhole collapse, which leads to rock fall and block gliding. Lastly, the further analysis such as loss estimation is required to support the risk evaluation, management, and reduction process.