Flood-Prone Area Assessment Using GIS-Based Multi-Criteria Analysis: A Case Study in Davao Oriental, Philippines

: Flooding is one of the major destructive natural disasters in Davao Oriental, Philippines, and results primarily from a high incidence of typhoons and heavy rainfalls. The main objective of this study was to identify ﬂood-prone risk areas by mapping them based on the integration of multiple indicators, including rainfall, slope, elevation, drainage density, soil type, distance to the main channel and population density. For this purpose, a GIS-based ﬂood risk spatial assessment was conducted by using analytic hierarchy process (AHP), weights by rank (WR) and ratio weighting (RW) frameworks to determine the relative importance of each indicator against another in the province of Davao Oriental. The resulting ﬂood-prone areas by the three methods are validated by comparing with the estimated ﬂood map based on ground truthing points from a ﬁeld survey. The comparison results show that AHP is the most appropriate method among them to assess ﬂood hazard. The result of the AHP ﬂood risk map shows that 95.99% (5451.27 km 2 ) of Davao Oriental is under low and moderate ﬂood risk. The high and very high ﬂood risk area covers approximately 3.39% (192.52 km 2 ) of the province, primarily in the coastal areas. Thirty-one out of the one hundred eighty-three (31 / 183) barangays (towns) are at a high to very high risk of ﬂooding at current climate, calling for the immediate attention of decision-makers to develop mitigation strategies for the future occurrence of ﬂooding in Davao Oriental.


Introduction
Flooding is considered one of the most serious and widespread natural hazards due to its devastating effects that endanger lives and cause property damage in the affected areas [1][2][3][4][5][6]. The changing climate behavior of extreme rainfalls and typhoons are the primary contributors to this problem. Human activities such as urbanization and the growth of settlements and assets in flooding areas likewise contribute to the increasing impacts of floods [1][2][3][4].
Phenomena commonly contributed to flooding are the limited capacity of river channels [7], human settlements in low-lying areas [8] and rapid growth of human settlements without upgrading the drainage infrastructures [1]. The problems related to flooding have increased extensively, and there is a need for efficient modeling to help mitigate the worst effects of flood disasters [9]. Thus, effective modeling will help in proper flood risk management planning and provides various insights into addressing the hazard and disaster problems. It is necessary to establish and implement a systematic approach to define areas that may have disastrous flood events.
The abovementioned methods successfully provided flood susceptibility assessment, but every method faced some limitations. In hydrological modeling, the preparation and calibration of parameters are time-consuming [26], the quality of DEM as input will affect the performance of the model [21,23,27], and it needs high requirements of computer resources. The most common and challenging problem in hydrological modeling during the analysis is also the scarcity or lack of hydro-meteorological data. In particular, data issues are critical in ungauged catchments. The area proximity, ordinary least squares and regression analysis are among the commonly used methods to solve this problem. However, as stated by Haddad et al. [31], the use of regression analysis has limitations, which involves the trading space for time. Yet, Haddad et al. [31] compared the two regression-based methods that employ a Bayesian generalized least squares framework: (i.e., quantile regression technique and parameter regression technique) to understand more the contributing factors in flooding. In quantitative modeling like statistical and data-driven methods increased the subjectivity as it requires an expert to select the flood conditioning factors. In the non-linear machine learning algorithms, it may lead to poor projections due to the large and inconsistent value ranges in the datasets [26,32]. In this study, the quantitative approach, AHP, was used for flood risk assessment due to data scarcity in the study area. Furthermore, other quantitative approaches such as WR, and RW were also applied for flood risk assessment, and their results were compared with that from AHP. These approaches have the capability to include multiple datasets and can be used on limited computer resources. Using multiple data sources potentially also allows for more reliable disaster risk reduction plan, especially in the context of flood management and mitigation. For multi-data sources, the multi-criteria analysis approach has become widely used to solve complex phenomena and to assess flood risk [1,4,30,33].
AHP is a popular and widely used method developed by Saaty [34] to solve a wide range of multi-criteria decision-making problems, combined with the use of a pairwise comparison matrix that calculates the weights of every criterion considered. Recently, researchers have given considerable attention to using AHP in flood risk assessment [1,4,13,35,36]. In those case studies, it has been shown that AHP can assess and map flood risk areas with reliable accuracy. However, the result is based on expert opinions and thus may be subjected to intellectual limitations due to uncertainty and subjectivity.
In Danumah et al. [1], AHP was used to integrate several elements under two criteria, hazard and vulnerability, for flood risk assessment and mapping. The slope, drainage density, type of soil, and isohyet were the elements used to generate their hazard map. The land use, population density, and drainage system were the elements used to create their vulnerability map. Danumah et al. [1] evidenced the reliability and the powerful role of geoinformation techniques in natural disaster assessment that requires the contribution of multi-source data. In Elkhrachy [2], a flash flood map was generated by using satellite images and the AHP was used to determine the relative impact or weight of flood-causative factors: the runoff, soil type, surface slope, surface roughness, drainage density, distance to the main channel and land use. Gigović et al. [4] conducted a study of the hazard zone mapping of urban flood-prone areas using a GIS multi-criteria methodology by considering six factors relevant to urban flooding: the height, slope, distance to the sewage network, distance from the water surface, water table and land use. Their results indicated that the scenario in which the methodology was used provides the highest level of compatibility with historical flooding data. The case studies of Danumah et al. [1], Gigović et al. [4] and Elkhrachy [2] were, however, limited to show the resulting flood hazard map or flood-prone area without validations of their resulting maps with ground trothing points or historical flood events. In this study, the resulting flood-prone area by AHP is validated by comparing with ground trothing points of historical flood events from a field survey. Moreover, the sensitivity test of consistency ratio is conducted to investigate the effect of varying ratio on flood-prone area.

Case Study Area
In the Philippines, the occurrence of natural hazards and disasters is frequent due to the physical environment of the country, which faces the Pacific Ocean where catastrophic typhoons originate from. Furthermore, the Philippines is located in a part of the Pacific Ring of Fire. The Philippines is vulnerable to typhoons, floods, earthquakes, storm surges and tsunamis [37]. The primary causes of the disasters in this area are typhoons and flooding because of their frequent occurrence and their magnitude of impact on society. Approximately 20 typhoons per year approach or make landfall, and this is the highest frequency of typhoon events in the entire world [9].
On the other hand, Mindanao Island, as described, is known to be rather typhoon free regions with less flood and storm surge risk in the Philippines. However, disasters due to tropical storm Sendong (international name, Washi) in 2011 and Super Typhoon Pablo (international name, Bopha) in 2012 changed that; Mindanao is now considered to be in one of the typhoon paths. These two deadly storms killed more than 1000 persons and led to approximately 100 missing persons because of heavy rains that triggered the rivers to overflow and that caused flash floods and landslides [38,39]. Davao Oriental is located in the southeast of Mindanao Island, and disasters due to floods and storm surge in this region are anticipated to be more frequent under global warming. Therefore, Davao Oriental is considered for the study area. the disasters in this area are typhoons and flooding because of their frequent occurrence and their magnitude of impact on society. Approximately 20 typhoons per year approach or make landfall, and this is the highest frequency of typhoon events in the entire world [9]. On the other hand, Mindanao Island, as described, is known to be rather typhoon free regions with less flood and storm surge risk in the Philippines. However, disasters due to tropical storm Sendong (international name, Washi) in 2011 and Super Typhoon Pablo (international name, Bopha) in 2012 changed that; Mindanao is now considered to be in one of the typhoon paths. These two deadly storms killed more than 1000 persons and led to approximately 100 missing persons because of heavy rains that triggered the rivers to overflow and that caused flash floods and landslides [38,39]. Davao Oriental is located in the southeast of Mindanao Island, and disasters due to floods and storm surge in this region are anticipated to be more frequent under global warming. Therefore, Davao Oriental is considered for the study area.
Davao Oriental is a province of Davao Region, Philippines. Davao Oriental is the easternmost province of the country and is located between 6°20' and 7°10' N latitude and 125°0' and 126°20' E longitude ( Figure 1). The province is composed of 183 barangays (towns) and two congressional districts. District 1, also known as the East Coast, consists of five municipalities; Tarragona, Manay, Baganga, Cateel and Boston. District 2, also called Davao Gulf, and includes four municipalities and one city: San Isidro, Governor Generoso, Banaybanay, Lupon and Mati City. Davao Oriental covers an area of approximately 5679 km 2 . The population is approximately 558,958 with a population density of 98 km 2 .  Davao Oriental is a province of Davao Region, Philippines. Davao Oriental is the easternmost province of the country and is located between 6 • 20' and 7 • 10' N latitude and 125 • 0' and 126 • 20' E longitude ( Figure 1). The province is composed of 183 barangays (towns) and two congressional districts. District 1, also known as the East Coast, consists of five municipalities; Tarragona, Manay, Baganga, Cateel and Boston. District 2, also called Davao Gulf, and includes four municipalities and one city: San Isidro, Governor Generoso, Banaybanay, Lupon and Mati City. Davao Oriental covers an area of approximately 5679 km 2 . The population is approximately 558,958 with a population density of 98 km 2 .
The topographic condition of Davao Oriental is characterized by a widespread chain of mountain ranges with an uneven distribution of plateaus, swamps and lowlands. The Mount Hamiguitan range is the newly enlisted UNESCO Heritage located at the administrative boundaries between the municipality of San Isidro, Governor Generoso and Mati City. This province occupies the largest land area of the provinces of Region XI (Davao Region): approximately 516,446 hectares or 26% of the total land area of Davao Region.
Davao Oriental is unique in that it is the only province in the country where all municipalities have a coastline. The coastline of the province measures 513.2 km from the municipality of Boston in the northern part of the province to the municipality of Banaybanay in the southwestern part of the province; it is the longest section of coastline in the country and approximately 3% of the total coastline of the country.

Rainfall
The observed rainfall data are obtained from the Hinatuan and DOST-RXI Stations ( Figure 2). The Hinatuan Station has observed daily rainfall data from 1973 to 2015. On the other hand, the DOST-RXI Station has observed daily rainfall data from 2006 to 2015. These two stations are far from the boundary of Davao Oriental. The distances of the Hinatuan and DOST-RXI Stations from the nearest border of the province are 43.9 km and 33.4 km, respectively, as shown in Figure 1D. Therefore, utilizing the rainfall data from these two stations will give an unreliable result because of their geographic locations.
To address the geographic limitation of the weather stations, rainfall data were obtained from the National Climatic Data Center (NCDC) and the Global Precipitation Climatology Centre (GPCC; see Table S1 in Supplementary Material 1). The NCDC rainfall data are global daily precipitations with a spatial grid resolution of 0.5 • latitude × 0.5 • longitude. In contrast, the GPCC rainfall data are global daily precipitations on a regular grid with a spatial resolution of 1.0 • latitude × 1.0 • longitude. The rainfall data were extracted within the boundary of Davao Oriental. Additionally, the nearest grid points were compared to the Hinatuan Station and DOST-RXI Station to evaluate the reliability of the data, as shown in Figure 1C,D. Figure 2 reveals that the NCDC rainfall data show good agreement with the observed rainfalls at the two stations. The GPCC rainfall, however, depicts lower values at Hinatuan while the rainfall at DOST-RXI shows good agreement with the observed and NCDC rainfalls. The rainfall data were interpolated as most of the related studies of flood hazard mapping used an interpolation method to create the rainfall distribution map [1,5,15,40].
The Hinatuan Station has observed daily rainfall data from 1973 to 2015. On the other hand, the DOST-RXI Station has observed daily rainfall data from 2006 to 2015. These two stations are far from the boundary of Davao Oriental. The distances of the Hinatuan and DOST-RXI Stations from the nearest border of the province are 43.9 km and 33.4 km, respectively, as shown in Figure 1D. Therefore, utilizing the rainfall data from these two stations will give an unreliable result because of their geographic locations.

Digital Elevation Model
The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) is one of the most widely used DEM datasets. ASTER GDEM has been applied in many fields, such as soil erosion, topography, geomorphology and hydrology [41,42]. ASTER GDEM is also widely used in developing flood hazard maps to extract drainage networks, flow accumulations and directions, basin boundaries, watershed boundaries, slopes and elevations [43][44][45][46].
The first version of ASTER GDEM, released in June 2009, was generated using stereo-pair images collected by the ASTER instrument onboard the Terra satellite. ASTER GDEM coverage spans from 83 • north latitude to 83 • south, encompassing 99% of the landmass of Earth [47]. The latest version of the ASTER GDEM V2 data was released on October 17, 2011, as shown in Figure 1D for the study area, Davao Oriental. The improved ASTER GDEM V2 dataset adds 260,000 additional stereo-pairs, improving coverage and reducing the occurrence of artifacts. The refined production algorithm provides improved spatial resolution, increased horizontal and vertical accuracy and superior water body coverage and detection [47]. Therefore, the ASTER GDEM V2 dataset is of better quality than the first version. ASTER GDEM V2 has a 30 m spatial resolution in the GeoTIFF image format with decimal degrees and WGS84 datum.

Administration Boundary
The Administrative boundaries of Davao Oriental include provincial, municipal and barangay boundaries. Davao Oriental is the easternmost province of the country. On the west side of Davao Oriental is the province of Compostela Valley, and the provinces of Surigao del Sur and Agusan del Sur are to the north. The Philippine Sea, part of the Pacific Ocean, is to the east of Davao Oriental. The administrative boundaries include ten municipal boundaries and 183 barangay boundaries. The administrative boundaries were provided as geographical information system (GIS) data from the global administrative areas and Philippine GIS organization, as indicated in Table S1. These GIS data are in decimal degrees and have a WGS84 datum. The data were then projected to the UTM coordinate system zone 51 N.

Population and Socio-Economic Data
Based on the 2015 census of population and housing [48] as shown in Table 1, the province of Davao Oriental has a total population of 558,958 in 2015. The 2015 census includes 41,340 more persons than were counted in the 2010 census, which determined a total population of 517,618 persons. This increase in the population from 2010 to 2015 translates into an average annual population growth rate of 1.47%. Mati City has the highest population of all municipalities, with 25.3% of the total provincial population. The municipality of Lupon is the second largest, with 11.8% of the total provincial population, followed by the municipalities of Baganga and Governor Generoso, with 10.1% and 9.9%, respectively. The rest of the municipalities contribute 43% of the total provincial population. Visayan migrants inhabit the province. Ethnic groups include the Mandaya, Mansaka, Manobo and Kalagan. The native languages spoken in the province are Bisaya, Mandaya, Dabawenyo and Chavacano. Due to the geographic location of Davao Oriental, the primary source of income is farming and fishing. The coastline of the province plays a vital role in the fishing activities of its fisherfolk.

Soil Type
The Davao Oriental soil cover is mainly loam and sandy clay loam, and a section of rough mountainous land has an unidentified soil type. The area of Davao Oriental is classified into four sets of input parameters based on United States Department of Agriculture Natural Resources Conservation Service (USDA-NRCS) [49]. The Camasan sandy clay loam and undifferentiated mountain soil are classified as sandy clay loam. The San Manuel silty clay loam and San Miguel clay loam are classified as clay loam. The Malalag loam is classified as loam, and the Bolinao clay is classified as clay.

Methodology
In this paper, the methodology is based on a GIS-based spatial assessment process for flood hazard and conducted by using multi-criteria analysis concepts and an AHP, WR and RW frameworks. This approach used the spatial data management capabilities of GIS and the flexibility of MCDA to combine factual evidence with value-based information. MCDA is a systematic procedure for designing, evaluating and selecting decision alternatives on the basis of conflicting criteria. The motivation of integrating GIS and MCDA stems from the need to make the geographic information technology more relevant for analyzing decisions. Spatial or GIS-based MCDA can be thought of as a collection method for transforming and combining geographic data and preferences (value judgments) to obtain information for decision-making. The factual evidences considered in this paper were slope, elevation, soil type, rainfall, drainage density and distance to the main channel, and defined as indicators. The value-based information approach, following Saaty [34], uses an expert decision to identify which indicators are the most crucial in flood hazard map. Figure 3 presents the flowchart of the methodology applied in this paper.

Indicators and Criteria
An important step of this analysis is the criteria selection for evaluating the flood risk map. The criteria considered are two, flood hazard and vulnerability. The barangay population density was used as the indicator in the vulnerability assessment. There are many indicators affecting flood hazard identification and modeling, and they vary from one study area to another. For instance, urban flood modeling is incredibly complex compared to rural flood modeling due to the interactions with manmade structures such as buildings, roads, channels, tunnels and underground structures. This paper used a composite flood hazard map index based on following six indicators. These indicators were selected based on various case studies [1,2,15,16,33] and based on the available data in the study area.
(1) Rainfall: At any location, the chance of flooding increases as the amount of rain increases. Higher rainfall intensity can result in more runoff because the ground cannot quickly absorb the water. In this study, the 5-day continuous maximum rainfall in every year was determined for 1991-2016 as shown in Figure 2. Then, the average of the 5-day continuous maximum rainfall over the 25 years was used in the analysis. Due to the limited number of weather stations in the study area, data from the NCDC was used in this study (Table S1). A spatial interpolation was performed to obtain the data for the study area by using the Kriging interpolation method, as in other case studies [5,[50][51][52]. The spatial distribution of the 5-day continuous rainfall is illustrated in Figure 4A.
(2) Slope: The slope is one of the crucial elements in flooding. The danger of flooding increases as the slope decrease. The slope is a reliable indicator of flood susceptibility [53]. The slope is

Indicators and Criteria
An important step of this analysis is the criteria selection for evaluating the flood risk map. The criteria considered are two, flood hazard and vulnerability. The barangay population density was used as the indicator in the vulnerability assessment. There are many indicators affecting flood hazard identification and modeling, and they vary from one study area to another. For instance, urban flood modeling is incredibly complex compared to rural flood modeling due to the interactions with manmade structures such as buildings, roads, channels, tunnels and underground structures. This paper used a composite flood hazard map index based on following six indicators. These indicators were selected based on various case studies [1,2,15,16,33] and based on the available data in the study area.
(1) Rainfall: At any location, the chance of flooding increases as the amount of rain increases. Higher rainfall intensity can result in more runoff because the ground cannot quickly absorb the water. In this study, the 5-day continuous maximum rainfall in every year was determined for 1991-2016 as shown in Figure 2. Then, the average of the 5-day continuous maximum rainfall over the 25 years was used in the analysis. Due to the limited number of weather stations in the study area, data from the NCDC was used in this study (Table S1). A spatial interpolation was performed to obtain the data for the study area by using the Kriging interpolation method, as in other case studies [5,[50][51][52]. The spatial distribution of the 5-day continuous rainfall is illustrated in Figure 4A.
Once the six indicators were defined, the next step was to build the spatial database. Each indicator was converted into raster data with a 30 m × 30 m grid resolution. The ASTER GDEM data were registered and projected to the UTM coordinate system zone 51N. The slope and elevation were obtained using the 3D Analyst algorithm based on a DEM. The drainage density and distance to the main channel were obtained using Arc Hydro, which is a set of data models that operate within ArcGIS to support geospatial and temporal data analyses. All data were integrated into the GIS environment using the AHP, WR and RW methods. Finally, the weighted overlay was used to calculate the flood hazard, which was then combined with the flood vulnerability to assess a flood risk in the study area.

Analytic Hierarchy Process
The AHP is a systematic approach developed in the 1970s to give decision-making based on experience, intuition and heuristics derived from sound mathematical principles [57]. It helps structure the decision-maker's thoughts and can help in organizing the problem in a way that is (2) Slope: The slope is one of the crucial elements in flooding. The danger of flooding increases as the slope decrease. The slope is a reliable indicator of flood susceptibility [53]. The slope is presented in Figure 4B.
(3) Elevation: In contrast, the elevation of an area is a major factor in flooding. Low elevation is a good indicator of areas with a high potential for flood accumulation. Water flows from higher to lower elevations; therefore, slope influences the amount of surface runoff and infiltration [5]. Flat areas at low elevations may flood more quickly than areas at higher elevations with steeper slopes [5]. The elevation is appeared in Figure 4C.
(4) Soil: Soil type is a significant factor in determining the water holding and infiltration characteristics of an area and consequently affects flood susceptibility [54,55]. Generally, runoff from intense rainfall is likely to be faster and greater in clay soils than in sand [56]. Additionally, rain runoff from intense rainfall is likely to be faster and greater in loam than in sand. The soil type of Davao Oriental is shown in Figure 4D.
(5) Drainage density: Drainage density is the length of all channels within the basin divided by the area of the basin [2]. A dense drainage network is a good indicator of flow accumulation pathways and of areas with a high potential for flooding [53]. The drainage density and distance to the main channel are exhibited in Figure 4E. (6) Distance to the main channel: Distance to the main channel significantly impacts flood mapping. Areas located close to the main channel and flow accumulation path are more likely to flood [4,53].
Once the six indicators were defined, the next step was to build the spatial database. Each indicator was converted into raster data with a 30 m × 30 m grid resolution. The ASTER GDEM data were registered and projected to the UTM coordinate system zone 51N. The slope and elevation were obtained using the 3D Analyst algorithm based on a DEM. The drainage density and distance to the main channel were obtained using Arc Hydro, which is a set of data models that operate within ArcGIS to support geospatial and temporal data analyses. All data were integrated into the GIS environment using the AHP, WR and RW methods. Finally, the weighted overlay was used to calculate the flood hazard, which was then combined with the flood vulnerability to assess a flood risk in the study area.

Analytic Hierarchy Process
The AHP is a systematic approach developed in the 1970s to give decision-making based on experience, intuition and heuristics derived from sound mathematical principles [57]. It helps structure the decision-maker's thoughts and can help in organizing the problem in a way that is simple to follow and analyze. Basically, the AHP helps in structuring the complexity, measurement and synthesis of rankings [58]. Therefore, these features make it suitable for a wide variety of applications, including alternative selection, resource allocation, forecasting, business process re-engineering, quality function development, benchmarking, public policy decisions, healthcare, hazard mapping, risk analysis and many more [57].
The AHP provides a means of decomposing the problem into a hierarchy of sub-problems, which can be more easily comprehended and subjectively evaluated. The subjective evaluations are converted into numerical values and processed to rank each alternative (indicator in this study) on a numerical scale. It can be explained in the following steps: Step 1: The problem is decomposed into a hierarchy of goal, criteria and indicators. Structuring the decision problem as a hierarchy is fundamental and the most important part in the process of the AHP. The hierarchy illustrates a relationship between elements of one level with those of the level immediately below. Figure 3 shows a tree structure representing an inverted hierarchy structure. At the root of the inverted hierarchy is the goal (to produce a flood risk map) of the problem being studied and analyzed. The leaf nodes are the six indicators to be compared.
Step 2: Data are collected from experts or decision-makers corresponding to the hierarchy structure, in the pairwise comparison of indicators based on a qualitative scale between 1 (equally important) and 9 (extremely important) proposed by Saaty [34] (see , Table S2 for the Saaty scale in the Supplementary Material). For each pairing within each criterion, a higher number means that the chose indicator is considered more important than the other indicator used in the comparison. Four local experts from the province of Davao Oriental in the field of natural disaster provided their judgments of the relative importance of one indicator against another (Table S3).
Step 3: The pairwise comparison of the six indicators generated at step 2 above are organized into a square matrix, the pairwise comparison matrix ( Table 2). The diagonal elements of the matrix are 1. The indicator in the ith row is more important that the indicator in jth column if the value of element (i,j) is larger than 1; otherwise the indicator in the jth column is more important than in the ith row. The (j,i) element of the matrix is the reciprocal of the (i,j) element. Step 4: The pairwise comparison matrix from step 3 is then normalized by Equation (1). Then, the priority vector (PV) is derived by dividing the sum of the normalized column of the matrix by the number of indicators, n, as shown in Equation (2).
where C ij is the value of an indicator in the pairwise comparison matrix, X ij is the normalized score and PV ij is the priority vector. The PVs give the relative importance of the indicators being compared and represent the weights of the indicators. Table 3 depicts the normalized comparison matrix and the PVs. Step 5: The consistency of the normalized comparison matrix is evaluated as described below. Since the numeric values of the matrix are derived from the subjective judgments of experts, it is impossible to avoid some inconsistency in the final matrix of their judgments [59]. Then, the question is how much inconsistency is acceptable. For this purpose, the AHP calculates a consistency ratio (CR) comparing the consistency index (CI) of the final normalized comparison matrix with that of a random-like matrix (RI). Saaty [34] provides the calculated RI value for matrices of different sizes (Table S4).
Firstly, use the PV as weights of each column of the pairwise comparison matrix as shown in Table S5. Secondly, multiply each value in the first column of the pairwise comparison matrix in Table  S5 by the PV of first indicator as shown in the first column of Table S6; continue this process for all the columns. Table S6 shows the resulting matrix. Thirdly, add the values in each row to obtain the weighted sum as shown in Table S7. Fourthly, divide the weighted sums by the corresponding PVs of each indicator to obtain the consistency measure (CM) as shown in Table S8. Fifthly, calculate the average of CMs, λ max . Lastly, calculate the CI and CR by the followings.
The λ max was 6.121 obtained from Table S8 and then the CI was 0.024 with n = 6. Therefore, by dividing the CI by RI (1.24 for n = 6 from Saaty), the CR was obtained as 0.02. Saaty suggests the value of CR should be less than 1 to have an acceptable consistency. Since the CR was 0.02 in this case, it can be assumed that the comparison matrix was reasonably consistent so that we could continue the decision-making process using AHP.
Step 6: By multiplying the six indicators by the corresponding weights and aggregating them, the hazard index (HI) is obtained by: where St, Sl, Dd, Dc, E and R represent the six indicators, soil type, slope, drainage, distance to main channel, elevation and rainfall, respectively.

Weights by Rank
In this method, the weight of indicator is arranging through a rank order. In that way, every indicator is ranked as per decision-maker preference. The decision-makers were the same four local experts considered in the AHP analysis. The most important is 1, the next is 2 and so on. The ranking of the indicators from the experts are shown in Table S9. Once ranking is established for a set of indicators, the numerical weights are calculated using Equation (6).
where W i is the normalized weight for ith indicator, n is the total number of indicators under consideration (j = 1, 2, . . . , n) and r j is the rank position of jth indicator. Each indicator is weighted (n − r i + 1) and then normalized by the sum of all weights (n − r j + 1). Therefore, the ranking method estimated weight should be considered as an approximation [60]. The results are given in Table 4. The average rank (AR) values are the average results from the four experts in the study area. Then, the resulting weight (W) in Table 4 was used to calculate the HI as shown in Equation (7). Table 4. The results of the weights by rank (WR) and ratio weighting (RW) derived from the ranking of the indicators from the same local experts considered in the AHP method. The W and AW are the resulting weights of each indicator used in the calculations of hazard index in Equations (7) and (10) (refer to Tables S9 and S10 in the Supplementary Material for derivation procedure).

Ratio Weighting
In this method, it also uses the pairwise comparisons to establish the relative importance among indicators. It is also likely that the pairwise comparisons will be inconsistent. Therefore, this approach requires more information to compensate for the inconsistencies of judgments. Efficient techniques to retrieve weights from pairwise comparison data include the simplified prioritization method as follows [61]: Step 1: A decision-maker assesses importance (weight) ratios between indicators using pairwise comparisons. In this study, the same local experts in Davao Oriental ranked the indicators. Table S10 shows the comparison matrix based on the expert's judgments among the indicators.
Step 2: Compute the geometric mean of each row of the comparison matrix in Table S10, and then normalize the resulting numbers. The geometric mean and the numerical weights are calculated using Equations (8) and (9). Table 4 depicts the resulting geometric mean and weights of the judgments.
where GM i is the geometric mean for ith indicator, n is the total number of indicators and W i is the normalized weight for ith indicator. The HI is computed by using the average weight (AW) in Table 4 as in Equation (10).
Finally, the HI obtained by the three models, AHP, WR and RW, was further processed using a weighted overlay analysis in ArcGIS. The results of weighted overlay analysis were classified using equal interval into four levels such as very low, low, moderate and high.

Validation
Super Typhoon Bopha was considered as one of the world's damaging storms in 2012. Bopha entered the Philippine Area of Responsibility (PAR) as a rapidly intensified Category 5 Super Typhoon and made its landfall at Baganga, Davao Oriental, of 4 December 2012. It has reached an average speed of 185 kph and gusts reaching 210 kph. Typhoon Pablo, as what Typhoon Bopha is called in the Philippines, was the most powerful storm to have hit the island of Mindanao, southern Philippines, in more than 100 years of recorded storms [60]. Its torrential rains generated massive debris flow in the Mayo River watershed in the Andap village in New Bataan municipality, causing some areas to be buried under a rubble as thick as 9 m and killing 566 people [60]. Moreover, it was also revealed that a storm surge with a height of approximately 1 m has reached a distance of 500 m inland away from the coastline [12]. Unfortunately, the government has no accounts in recording the flooded areas during the said event.
Due to the hazards caused by storm surge, storm surge prediction and forecasting currently is actually progressing. In fact, in a study conducted by Ross, et al. [60] regarding the projected susceptibility to storm surge modeling using historical records from Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA), it showed that there was a storm surge event as high as 3.66 m. The said storm surge happened in the city of Bislig, located about 60 km from the northern coastal boundary of the Davao Oriental province. Hence, to address such scenarios, this study was conducted using a GPS-based field survey to determine the flooded areas during Typhoon Bopha.
In order to validate the resulting flood-prone areas by the AHP, WR and RW methods, a GPS-based field survey to local people was carried out along the east coast of Davao Oriental to investigate ground true flooded points using historical flooding events. This survey was carried out to determine the flooded area during the Typhoon Bopha. It is to note that the field survey was also not carried out in the municipalities located near the Davao Gulf thus, the points do not cover all barangays specified in the study area. Hence, only 70 ground truthing points were collected and were used for performance evaluation (see Figure 1 for the survey result).
The performance evaluation was done using the accuracy assessment of the flood classification. One of the commonly used methods is to apply a confusion matrix or error matrix. This method can be used to compute several assessment elements such as overall accuracy (ACC), true positive rate (TPR), true negative rate (TNR), false positive rate (FPR) and false negative rate (FNR) using the Equations (11)-(15), respectively.
where P, N, TP, TN, FP and FN are condition positive, condition negative, true positive, true negative, false positive and false negative, respectively. To evaluate the performance, an accuracy assessment was performed using the confusion matrix, and it was calculated based on the 70 ground points collected from the field survey. Elevations were extracted from these truthing points, and it showed that it is in the range between 10 m and 20 m above mean sea level. After spatially interpolating the points, it was identified that all areas in Davao Oriental that belongs to the below 20 m elevation were considered as susceptible flood-prone areas.
Based on the survey, a resulting map was created (as seen in Figure 1C) and was overlaid with the prediction maps obtained by the AHP, RW and WR methods to compare and assess the accuracy of the prediction maps. In the accuracy assessment, the TP and FP are the numbers of pixels correctly classified and the numbers of pixels mistakenly classified, respectively. With this, the results showed that the number of correctly predicted pixels by the model (TP) are plotted to the number of incorrect predicted pixels (FP), as illustrated in Table 5 and Figure 5 , indicating that the AHP-based flood area is truly flooded in the ground truthing-based result. Table 5. Confusion matrix to calculate the accuracy (ACC), true positive rate (TPR), true negative rate (TNR) and false positive rate (FPR) for performance evaluation of three methods, AHP, WR and RW.  flood area that are shown in Figure 1C. In the overall assessment, AHP had superior results in ACC, TNR and TPR as shown in Figure 5A. Moreover, Figure 5B shows the accuracy of the positive prediction (TPR) compared to the wrong prediction (FPR). The upper part of the red line in the Figure  5B means the model was accurate, otherwise not accurate. Therefore, AHP shows higher accuracy compared to the other two methods. In overall, the analysis reveals that AHP was much more accurate and reliable for flood risk analysis in this study.  Data scarcity and sparsity are the main problems solved in this approach. In this paper, data scarcity means limited historical flood events, or the recorded floods are limited as compared to the amount needed in the validation of the model (see Figure 1D). Moreover, as shown in Figure 1D, the survey data of flood are not evenly distributed (i.e., sparsity) to the study area. In this study, the flood data used was derived from the flooding caused by typhoon Bopha and did not cover all flood events happened in Davao Oriental. This approach has uncertainty and limitation, same with other hydrological modeling. Fortunately, the GIS-based flood validation approach is a good start for the researcher to conduct flood susceptibility mapping in data-scarce and sparse region.

Vulnerability Map
Vulnerability expresses the level of inability to resist a hazard or to respond when disaster has occurred [62]. For example, people who live in low-lying areas are more vulnerable to floods than people who live at higher elevations. Moreover, vulnerability is the most crucial component of flood risk because vulnerability determines if exposure to a hazard constitutes a threat. Flood vulnerability mapping is the process of determining the degree of susceptibility and exposure of a given place to flooding. The susceptibility and exposure issues include several factors, such as the age and health of the population, the socio-economic activities, the quality of buildings and their location with respect to flooding. In this study, due to the limited socio-economic data from the provincial office of the Philippine Statistics Authority in Davao Oriental, the only indicator used for the vulnerability map is the population density from barangay (town). Thus, the vulnerability map obtained from the barangay population density. Vulnerability map is divided into five classifications (shown in Figure 4F). Barangays with high populations are mostly proximal to coastal areas and riverside.

Hazard Map
Hazard is the function of scale and probability of natural disaster to cause harm or loss in a specified location in a period of time. Thus, it reflects the natural properties of the disaster [36], and the physical and statistical aspects of the flooding process [63]. Flood hazard maps contain information about the probability and/or magnitude of an event [63], the return period of the flood, the extent, and depth of inundation [64]. Due to the exclusion of the hydrological modeling of this study, the depth of inundation, and the return period of the flood were not included. Thus, in this study, the flood hazard map only covered the extent or the possibility of the spatial distribution of flood in the study area as shown in Figure 4G. Figure 4G shows that the hazard map was a combination of maps: rainfall ( Figure 4A), slope and elevation ( Figure 4B,C), soil ( Figure 4D) and drainage density and distance to the main channel ( Figure 4E). Each map had different weights that, in this paper, the rainfall was emphasized more. The rainfall map ( Figure 4A) shows that each east coast municipality had a higher incidence of rainfall compared to that of the Davao Gulf municipalities. The slope and elevation maps ( Figure 4B,C) indicate that most of the areas in Davao Oriental were at low elevations, less than 217 m and have slopes up to 16 degrees. All maps were classified based on equal interval classification.
The resulting flood hazard maps by the AHP, WR and RW methods are illustrated in Figure 6. Table 5 shows the confusion matrix of the three models. The data in Table 5 was the number of pixels from the models (i.e., predicted) in high classification, and the observed data from the flood and non-flood area that are shown in Figure 1C. In the overall assessment, AHP had superior results in ACC, TNR and TPR as shown in Figure 5A. Moreover, Figure 5B shows the accuracy of the positive prediction (TPR) compared to the wrong prediction (FPR). The upper part of the red line in the Figure 5B means the model was accurate, otherwise not accurate. Therefore, AHP shows higher accuracy compared to the other two methods. In overall, the analysis reveals that AHP was much more accurate and reliable for flood risk analysis in this study.
Finally, in AHP, the very low and low classes covered approximately 1% and 22% of the total area of Davao Oriental, respectively. These are areas with high elevations, high slopes and low drainage densities. The categories of moderate, and high classes were estimated to cover 55%, and 22% of the total area of Davao Oriental, respectively, and are mostly in the east coast municipalities (Boston, Cateel and Baganga). For the high hazard risk areas of 22.15% of Davao Oriental, the rainfall (41.8%) and slope (23.8%) were the most significant causative factors of flood occurrence. In the flood-prone area, 64 out of 183 barangays (towns), and 30% (168,034/558,958) of the population were at a high of flooding at the current situation as shown in Table 6. Additionally, these were areas that had high occurrences of rainfall, high drainage densities, low elevations and low slopes, and they were proximal to river channels and shorelines. In the Davao Gulf municipalities, most of the areas of Mati City and Lupon belonged to the moderate hazard class due to their geophysical characteristics. Lupon had one of the widest rivers of all the municipalities in the province. On the other hand, Mati City was a low-lying city, especially in the city centre where the urbanized area was located. Finally, in AHP, the very low and low classes covered approximately 1% and 22% of the total area of Davao Oriental, respectively. These are areas with high elevations, high slopes and low drainage densities. The categories of moderate, and high classes were estimated to cover 55%, and 22% of the total area of Davao Oriental, respectively, and are mostly in the east coast municipalities (Boston, Cateel and Baganga). For the high hazard risk areas of 22.15% of Davao Oriental, the rainfall (41.8%) and slope (23.8%) were the most significant causative factors of flood occurrence. In the floodprone area, 64 out of 183 barangays (towns), and 30% (168,034/558,958) of the population were at a high of flooding at the current situation as shown in Table 6. Additionally, these were areas that had high occurrences of rainfall, high drainage densities, low elevations and low slopes, and they were proximal to river channels and shorelines. In the Davao Gulf municipalities, most of the areas of Mati City and Lupon belonged to the moderate hazard class due to their geophysical characteristics. Lupon had one of the widest rivers of all the municipalities in the province. On the other hand, Mati City was a low-lying city, especially in the city centre where the urbanized area was located.

Risk Map
In the study of Zhou et al. [36], the risk is the product of hazard and vulnerability as defined by the United Nations, Department of Humanitarian Affairs [65,66] and World Meteorological

Risk Map
In the study of Zhou et al. [36], the risk is the product of hazard and vulnerability as defined by the United Nations, Department of Humanitarian Affairs [65,66] and World Meteorological Organization [67]. A study by Thieken et al. [63] also defined flood risk as to the product of hazard and vulnerability.
The flood risk assessment has been widely used in understanding flood disaster evaluation, warning and awareness [36]. It is essential in management and decision making. In this study, flood risk map is a combination of hazard map and vulnerability map. The result of the flood risk map was generated by a weighted overlay of the hazard and vulnerability maps with equal weights. The flood risk map was classified into five classes ( Figure 4H). The very low, low and moderate classes covered 0.63%, 60.64% and 35.35% of the total area, respectively. The categories of high and very high were estimated at 3.34% and 0.05% of the total area, respectively. These areas were barangays with high population densities and were mostly in low-lying areas near the shoreline and rivers. On the other hand, the moderate and very low classes were areas at high elevations with smaller populations. Thirty-one out of one hundred eighty-three (31/183) barangays were at high risk of flooding. Only 1 barangay had a very high risk of flooding at the current climate, calling for the immediate attention of the government to develop strategies or flood protection for the future occurrence of flooding within the province of Davao Oriental.

Effects of Varying Consistency Ration (CR) on Flood Hazard Map
Integrating multiple indicators (data sources) in the multi-criteria analysis presents a real advantage, and results show that the AHP approach allowed a better understanding of all the indicator contributions in the flood assessment because a weight was given to each indicator. However, data from different sources with different resolutions were factors of bias during processing and analysis. On the other hand, the addition of weights reduced the bias and uncertainty in the result.
However, the AHP method shows some limitations due to the subjectivity in choosing the value of the indicator weighting from the arbitrary judgments of experts [1,6,13]. The CR is the crucial element when generating HI for the multi-criteria analysis using AHP method. According to Saaty [34], this subjectivity problem was reduced by the CR test. The CR threshold must be less than 0.10 or 10% to make a coherent judgment. The lower the CR, the better, and if CR is 0, it means the HI is perfect. Unfortunately, the range of CR between 0 and 9.9% is a rough estimation. The HI is coherent if the CR is 9.9% or CR is 0%, but the weights of the indicators vary in every CR generation.
Therefore, the sensitivity of HI was tested in different three CR scenarios; (1) the lowest CR was 2.0%, (2) the highest CR was 8.6% and (3) the middle CR was 5.8%. Figure 7 shows the flood hazard maps from the sensitivity analysis in different CR scenarios. Table 7 shows the weights of every indicator according to the ratio of indicators. Table 8 shows the rank of every classification in different CR. The result illustrates that the ranking of classifications is the same in the three CR scenarios. Average change according to classification was also ranging from −0.77% to 1.21%, and the overall average change was only 0.02%. Therefore, it was found that the changes of the CR did not significantly affect the result of the hazard map, and the changes were acceptable. Thus, the lowest criteria ratio was more appropriate to use as according to Saaty [34] so that the lower CR made the  The result illustrates that the ranking of classifications is the same in the three CR scenarios. Average change according to classification was also ranging from −0.77% to 1.21%, and the overall average change was only 0.02%. Therefore, it was found that the changes of the CR did not significantly affect the result of the hazard map, and the changes were acceptable. Thus, the lowest criteria ratio was more appropriate to use as according to Saaty [34] so that the lower CR made the judgment more coherent. Therefore, the CR with 0.02 was used for to obtain the flood risk map.

Conclusions
Flooding is a natural hazard that poses a risk to both population and physical structures within the affected areas. Flood risk occurrence can be better understood with the use of a spatial assessment. Thus, in this study, the main objective was to identify flood-prone and categorize risk areas to produce a flood risk assessment based on several factors. For this purpose, a GIS-based flood risk assessment was conducted by using the concepts of AHP, WR and RW frameworks in the province of Davao Oriental, and determines which model has a superior result.
The analysis was based on the integration of multiple indicators including rainfall, slope, elevation, drainage density, soil type and distance to the main channel. The AHP model depicts a superior result in the computation of relative importance of indicators compared to WR and RW methods. The result of the AHP model shows that Davao Oriental is under high flood-prone areas approximately 22.15% of the province, and these are in low-lying areas near the coastline and rivers. The multi-criteria analysis allowed the integration of the several indicators under two criteria for flood risk assessment: hazard and vulnerability. The results showed that the province of Davao Oriental is moderately exposed to the risk of flooding, which requires immediate attention from the decision-makers to develop strategies for the future occurrence of flooding within the province of Davao Oriental. An effective strategy should be drawn up at the barangay level since each barangay has its own physical and social characteristics.
Immediate strategies like installing weather measurement instruments in at least one station per municipality will give us precise weather forecasting especially during extreme events like drought, heavy rainfalls, storms and typhoon. Since the national government has no specific flood-risk mitigation plan, the municipal government in coordination with the office of the disaster risk reduction and management council can design a flood-risk education program in every barangay. The flood-risk education program is an education campaign on how to act before, during and after the event of flooding.
Additionally, the results may be improved by adding land use/land cover through image analysis methods that use high-resolution images such as Landsat 8+. Hydrologic modeling in 2D/3D for efficient analysis in determining the depth of flood [19], the extent of flood and the damage estimate using damage curve of the flood is also recommended as shown in literatures [20,21,23,24,27].
Lastly for future works, in the selection of methods, the subjectivity of the quantitative approaches such as AHP, WR and RW can be improved using data mining and machine learning techniques.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/11/11/2203/s1, The Tables S1-S10 can be found in the Supplementary Material 1. The Supplementary Material 2 is the Excel file demonstrating the AHP procedure to calculate the weighting factors used in this study. The Supplementary Material 3 includes the field survey result.
Author Contributions: J.S.C. developed the concept of this study under the supervision of H.S.L. The analysis was carried out by both authors. Both authors drafted the first version of the manuscript and worked on improving and finalizing the manuscript.
Funding: This research is partly supported the Grant-in-Aid for Scientific Research (17K06577) from JSPS, Japan.