Flood risk assessment for Davao Oriental in the Philippines using geographic information system‐based multi‐criteria analysis and the maximum entropy model

The assessments of flood‐prone areas and flood risk due to pluvial flooding for Davao Oriental on Mindanao Island in the Philippines were carried out by the analytic hierarchy process (AHP) and maximum entropy (Maxent) models using multiple criteria such as slope, elevation, soil type, rainfall, drainage density, distance to the main channel, and population density. Flood records from 70 survey points were obtained and used to verify the model results. The criteria weights of the top three important factors in the AHP are rainfall (42%), slope (23%), and elevation (15%), whereas those in the Maxent model are elevation (36%), rainfall (23%), and soil (19%). The verification results show that the accuracies of the AHP and Maxent model are 81 and 95.6%, respectively, indicating that both approaches are reliable in flood hazard and risk assessments. Approximately 22% of the total area and approximately 30% of the total population of Davao Oriental are classified as high risk of pluvial flooding in the current situation by the AHP method. This study shows a broad‐scale high‐level data‐driven screening method that can be used to help identify potential hot spots for pluvial flooding for which more detailed numerical modelling studies should be undertaken.


| 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. Human activities such as urbanisation and the growth of settlements and assets in flooding areas likewise contribute to the increasing impacts of floods (Danumah et al., 2016;Gigovi c, Pamučar, Baji c, & Drobnjak, 2017). Phenomena that commonly contribute to flooding are the limited capacity of river channels, human settlements in low-lying areas, and rapid growth of human settlements without upgrading the drainage infrastructures.
Effective modelling of floods will help in proper flood risk management planning and provide various insights into addressing the hazard and disaster problems (Otieno, 2004). Therefore, it is necessary to establish and implement a systematic approach to define areas that may have disastrous flood events. Hazard mapping is a vital component of flood risk management. A flood hazard map shows the extent of the water level in a flooded area. Another critical element in flood risk management is flood risk mapping, which shows the potential risk to the population, economy, and environment due to flooding (Gigovi c et al., 2017).
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; 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 (Gurenko & Lester, 2004). The primary causes of the disasters in this area are typhoons and flooding because of their frequent occurrence and their magnitude of impact to society. Approximately 20 typhoons per year approach or make landfall, and this is the highest frequency of typhoon events in the world (Otieno, 2004).
Mindanao in the Philippines is known to be rather typhoon-free. However, disasters due to tropical storm Sendong (international name, Washi) and Super Typhoon Pablo (international name, Bopha) have changed this such that Mindanao is now considered to be in one of the typhoon paths. Sendong and Pablo made landfalls in 2011 and 2012, respectively. These two deadly storms killed more than 1,000 people and led to 100 missing persons because of heavy rains that triggered the rivers to overflow and that caused flash floods and landslides (National Disaster Risk Reduction andManagement Council, 2012a, 2012b).
Although flooding is a severe hazard in the province of Davao Oriental in southeastern Mindanao Island due to typhoons and heavy rainfall, insufficient attention has been paid to flood hazard assessment. Recent scientific works undertaken in Davao Oriental by Lagmay, Racoma, Aracan, Alconis-Ayco, and Saddi (2017) and Ross, Santiago, and Lagmay (2015) were concentrated on pluvial flood mapping in major river basins using 2D flood simulations and vulnerability assessment against pluvial flooding and coastal flooding due to storm surge during the passage of Typhoon Pablo. Ross et al. (2015) focused on the east coast municipalities, but Typhoon Pablo affected the entire province. All the floodingrelated studies in Davao Oriental are beneficial to providing a spatial presentation of the distribution of the flooded areas. However, no studies have yet undertaken the evaluation and flood risk mapping of the entire province of Davao Oriental using multiple datasets.
The abovementioned methods successfully provided a flood risk assessment, but every method faced some limitations. In hydrological modelling, the preparation and calibration of parameters are time-consuming (Hong, Panahi, et al., 2018) and have high computer resource requirements. Quantitative modelling, such as statistical and data-driven methods, increases the subjectivity and uncertainty because they require an expert to select the flood conditioning factors (Cabrera & Lee, 2019). Nonlinear machine learning algorithms may lead to poor projections due to the large and inconsistent value ranges in the datasets in data-scarce areas (Hong, Panahi, et al., 2018;Tien Bui et al., 2016).
Thus, the objective of this study is to identify floodprone areas and to assess flood risk due to pluvial flooding in a data-scarce area, Davao Oriental in Mindanao Island, Philippines, based on a broad-scale high-level data-driven screening method to help identify potential hot spots for pluvial flooding for which more detailed numerical modelling studies have to be undertaken. For this purpose, GIS-based processes for the spatial assessment of the pluvial flood-prone hazard map were conducted by using the AHP in a quantitative approach and maximum entropy (Maxent) model in machine learning with multiple datasets including rainfall, slope, elevation, drainage density, soil type, and distance to the main channel.
The study area is described in section 2. In section 3, the datasets and criteria selection are described. The multi-criteria decision analysis (MCDA) methodology by the AHP and Maxent model used in this study and validation process are presented in section 4. The results and discussion are presented in sections 5 and 6, respectively. Finally, the conclusions are described in section 7.

| STUDY AREA
Davao Oriental is a province in the Davao Region, Philippines. Davao Oriental is the easternmost province of the country and is located between 6 20 0 and 7 10 0 N latitude and 125 0 0 and 126 20 0 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, includes four municipalities and one city: San Isidro, Governor Generoso, Banaybanay, Lupon, and Mati City. Davao Oriental covers an area of approximately 5,679 km 2 . The population is approximately 558,958 (Philippine Statistics Authority, 2015) with a population density of 98/km 2 .
The topographic condition of Davao Oriental is characterised 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 site 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 ha or 26% of the total land area of the Davao Region. 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 southwest, with approximately 3% of the total coastline of the country.

| MATERIALS
The multiple datasets considered in this paper were slope, elevation, soil type, rainfall, drainage density, and distance to the main channel and were defined as criteria or factual evidence. Figure 2 presents the flowchart of MCDA for pluvial flood risk assessment.

| Rainfall
The observed rainfall data are obtained from the Hinatuan and DOST-RXI Stations ( Figure 1d). The Hinatuan Station has observed daily rainfall data from 1973 to 2015, and 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 and 33.4 km, respectively ( Figure 1d). Therefore, utilising 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 1). The NCDC and GPCC rainfall data are global daily precipitation with spatial grid resolutions of 0.5 × 0.5 and 1.0 × 1.0 . The rainfall data were extracted within the boundary of Davao Oriental. Additionally, the nearest grid points were compared to the Hinatuan and DOST-RXI stations to evaluate the reliability of the data (Figure 1c,d). Figure 3 reveals that the NCDC rainfall data show better agreements with the observed rainfalls at the two stations. Thus, the NCDC rainfall data are used in the analysis of the study.

| Digital elevation model
Topography is arguably the key factor for the estimation of flood extent (Hawker, Bates, Neal, & Rougier, 2018). Therefore, generally speaking, the quality of flood predictions increases with higher-resolution digital elevation model (DEMs). Higher-resolution DEMs are more important when modelling urban environments with buildings to be captured (Fewtrell, Bates, Horritt, & Hunter, 2008). However, resolution can be less important for rural environments. For example, Savage, Bates, Freer, Neal, and Aronica (2016) showed that running simulations finer than 50 m had little performance gain without incurring additional unnecessary computational cost.
The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM2) is one of the most widely used DEM datasets. The overall vertical accuracy of GDEM2 is 17 m at the 95% confidence level, whereas that of GDEM1 is approximately 20 m at the 95% confidence level. ASTER provides improved spatial resolution from 120 m in GDEM1 to 71-82 m in GDEM2, provides reduced horizontal displacement from 0.95 pixels in GDEM1 to 0.23 pixels in GDEM2, and GDEM2 provides better water body coverage and detection that that of GDEM1 (Tachikawa et al., 2011). The DEM dataset has been applied in many fields, such as soil erosion, topography, geomorphology, hydrology, and flood hazard maps (Chatterjee, Huong & Nagasawa, 2014;Othman, Jafri, Lim, & Tan, 2011;Ramani & Balasubramaniam, 2013;Reddy, Kumar, Sahu, & Singh, 2017). It has a 30 m spatial resolution in the GeoTIFF image format with decimal degrees and WGS84 datum and is used for the study area, Davao Oriental, as shown in Figure 1d.

| Administration boundary
The administrative boundaries of Davao Oriental include provincial, 10 municipal and 183 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 data were provided as shapefiles from the global administrative areas and Philippine GIS organisation (Table 1). These shapefiles are in decimal degrees and have a WGS84 datum. The data were then projected to the Universal Transverse Mercator (UTM) coordinate system zone 51 N. 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 T A B L E 1 Summary of the dataset used in this study: rainfall, DEM, administration map, soil type, population, and socio-economic data the total provincial population, followed by the municipalities of Baganga and Governor Generoso with 10.1 and 9.9% of the total provincial population, 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 (Philippine Statistics Authority, 2010). Because of the geographic location of Davao Oriental, the primary sources of income are 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 two sets of input parameters based on USDA-NRCS (1986). 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.

| Criteria selection
There are many criteria affecting flood hazard identification and modelling, and they vary from one study area to another. This paper used a composite pluvial flood risk map index based on seven criteria: rainfall, slope, elevation, drainage density, distance to the main channel, soil type, and population density. Although Davao Oriental is a coastal province, the distance to the coast is not considered because this study focuses only on pluvial flooding due to heavy rainfall. For a vulnerability map, only F I G U R E 3 Comparisons of the observed annual rainfall at the (a) Hinatuan and (b) DOST-RXI stations with the NCDC and GPCC rainfall data at DO63 and DO20 points, respectively. The annual maximums of 5-day continuous rainfall at DO63 and DO20 are also presented in green. The locations of DO20 and DO63 are 7.083 N, 125.95 E, and 8.00 N, 126.33 E, respectively population density is considered. Due to the lack of statistical records, property, infrastructure, and other critical criteria are not considered. These criteria were selected based on various case studies (Danumah et al., 2016;Di Baldassarre, Castellarin, Montanari, & Brath, 2009;Elkhrachy, 2015;Wang, Li, Tang, & Zeng, 2011;Yahaya, 2010) and based on the available data in the study area.
Rainfall (Figure 4a) The 5-day continuous maximum rainfall in every year was determined at all data points in Davao Oriental for 1991-2016 from the NCDC rainfall data ( Figure 3). Then, the average of the 5-day continuous maximum rainfall over 25 years was used in the analysis. Finally, a spatial interpolation was performed to obtain the spatial rainfall pattern for the study area by using the kriging interpolation method, as in other case studies (Basconcillo et al., 2016;Kazakis, Kougias, & Patsialis, 2015;Shahabi, Salari, Ahmad, Bin, & Mohammandi, 2016;Zhao, Wang, Cheng, Liu, & He, 2015).

Slope (Figure 4b)
The slope is the rate of change of the surface in the horizontal and vertical directions from the centre of the grid. The slope is measured in units of degrees and is one of the crucial elements in flooding. The risk of flooding may increase as the slope increases. The slope is a reliable indicator of flood susceptibility (Bapalu, 2006;Islam & Sado, 2000).
Elevation (Figure 4c) The elevation of an area is also a major factor in flooding. Low elevation is a good indicator of areas with a high potential for flood accumulation. Flat areas at low elevations may accumulate water more quickly than areas at higher elevations with steeper slopes (Kazakis et al., 2015).   Falaky, 2016;Nyarko, 2002). Generally, runoff from intense rainfall is likely to be faster and greater in clay soils than in sand (The University Corporation for Atmospheric Research, 2010). Additionally, rainfall runoff from intense rainfall is likely to be faster and greater in loam than in sand.
Drainage density (Figure 4e) Drainage density is the length of all channels within the basin divided by the area of the basin (Elkhrachy, 2015). A dense drainage network is a good indicator of flow accumulation pathways and areas with a high potential for flooding (Islam & Sado, 2000).
Distance to the main channel ( Figure 4e) Pluvial flooding occurs when extreme rain saturates the drainage systems, and the surplus water cannot be absorbed. Thus, the 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 (Gigovi c et al., 2017;Islam & Sado, 2000).
Once the six criteria were defined, the next step was to build the spatial database. Each criterion was converted into raster data with a 30 × 30 m grid resolution. The ASTER GDEM V2 data were registered and projected to the UTM coordinate system zone 51 N. 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 method. Likewise, in the Maxent modelling, all criteria are considered environmental layers, and the survey of the flooding events is the sample layer. The sample and environmental layers are combined to predict the pluvial flood-prone areas.

| METHODOLOGY
In this paper, the methodology is based on GIS-based spatial assessment processes for flood hazard and risk assessment by using the AHP and the Maxent model.

| AHP modelling process
The multi-criteria analysis was used to analyse a series of alternatives or objectives by ranking them from the most preferable to the least preferred using a structured approach (Papaioannou, Vasiliades, & Loukas, 2015). The AHP introduced by Saaty (1980) is a multi-criteria decision-making approach, a decision support tool, and it is used to solve complex decision problems. The AHP uses a hierarchical approach to represent a problem, and this hierarchy can be summarised as three main levels: goals, criteria, and indicators. The evaluation of indicators and their weights must be determined according to their importance using paired comparisons. The AHP consists of six steps (Danumah et al., 2016;Saaty, 1980). First, a complex, unstructured problem is decomposed into its component factors in a hierarchy of goal, criteria and indicators. Second, data are collected from experts or decision-makers corresponding to the hierarchy structure in the pairwise comparison of indicators. Third, the pairwise comparison of the six indicators generated at the previous step are organised into a square matrix, the pairwise comparison matrix. Fourth, the pairwise comparison matrix from the previous step is normalised, and then the relative weights of each criterion are calculated. Fifth, the priority variables are determined by the synthesised judgements from the normalised comparison matrix. Finally, the consistency of the assessments and judgements from the normalised comparison matrix is evaluated. In the following, the key steps are described in detail.

| Pairwise comparison
The key step in the AHP is to make a pairwise comparison between each criterion based on the scales by Saaty (1980). The survey was conducted among the four stakeholders in the study area who are involved in the field of natural disasters in the city and provincial planning. The four local experts provided their judgements of the relative importance of one indicator against another. Then, the average results of the survey were used to determine which criterion has high importance and vice versa. The results of the comparison were described regarding integer values from 1 to 9, where a higher number means that the chosen criterion is considered more important than the other criterion used in the comparison. The results of the pairwise comparison matrix are shown in Table 2.

| Normalisation
This step is the process of normalising the matrix by adding the numbers in each column. Each entry in the column is then divided by the column sum to yield its normalised score as described in Equation (1). The sum of each column should be 1. Finally, the priority vector (PV) is computed by dividing the sum of the normalised column of the matrix by the number of criteria used, n, as shown in Equation (2). Table 3 shows the normalised matrix in this study.
where C ij is the value of a criterion in the pairwise comparison matrix, X ij is the normalised score, and PV ij is the priority vector of a criterion.

| Consistency analysis
There are three steps to determine the consistency ratio. First, the consistency measure is calculated by multiplying the pairwise matrix by the PV, and then the result is divided by the weighted sum vector with the criterion weights. Second, the consistency index (CI) is calculated, as described in Equation (3). Last, the consistency ratio (CR) is computed by Equation (4), where λ max is the sum of the consistency measure (CM) divided by n and set to 6.126. Since the number of criteria in this study is 6, the random index (RI) is 1.24, as given in the random index table of Saaty (1980).
The key component of the AHP is the calculation of the CR. If the CR exceeds 0.1, the set of judgements may be too inconsistent to be reliable. If the CR is less than 0.1, then the comparison matrix can be considered to have an acceptable consistency. A CR of 0 means that the judgements are entirely consistent.
In Table 3, the column PV contains the relative importance weights for each criterion. From the input values in the pairwise comparison and the weights calculation, the CR was found to be 0.02. Regarding the subjectivity issue raised in the AHP, CR sensitivity tests were conducted to reduce the subjectivity problem (Cabrera & Lee, 2019). The CR result indicates a reasonable level of coherency in the pairwise comparison.
Then, the hazard index (HI) was used to consider the rate of probability and was calculated based on where St, Sl, Dd, Dc, E, and R represent the soil type, slope, drainage, distance to the main channel, elevation, and rainfall, respectively.

| Maximum entropy model
The entropy algorithm to describe the uncertainty or information content of a random event quantitatively is described by Shannon (1948) as where H is the information entropy, p i is the occurrence probability of a random event, and C is a positive constant. The entropy H is a function of p i , and thus, under given experimental conditions, a distribution exists that maximises H. This distribution has a dominant probability and is the most common distribution, so it is called the "most probable distribution." To summarise, the principle of maximum entropy is to choose the distribution under given restrictions when the entropy is the maximum value in all possible compatible distributions (Feng & Hong, 2009). Based on the principle of maximum entropy and the Lagrange undetermined multiplicator method, the distribution when the entropy is maximum can be obtained (Feng & Hong, 2009). Let random variable x be x 1 , x 2 , …, x n , and the corresponding probabilities are p 1 , p 2 , …, p n . They could satisfy The mean value F k of several known functions f k (x i ) is given as To discover the distribution when entropy is the maximum value under the constraint condition of Equations (7) and (8), undetermined multiplicators α and β k are introduced to form a new function: Using the inequality lnx ≤ x − 1, Equation (9) is changed into: To have H be the maximum value, the above formula must be an equation; then, p i should satisfy the following equation: Using Equation (7), Equation (10) To obtain the value of β k , substitute Equation (11) into constraint Equation (8) to obtain In Equation (12), both F k and f k (x i ) are known, while the real unknowns are m values of β(β 1 , β 2 , …, β m ). The M equations could obtain m β values; thus, we could have the value of p i when entropy is the maximum value (Feng & Hong, 2009). The above calculation formula obtained from discrete conditions could also be used in the calculation process in continuous conditions (Feng & Hong, 2009). All criteria and data related to the locations of the flooding events are used as inputs in Maxent (Steven, Miroslav, & Schapire, 2017) and in ArcGIS software for mapping the potential pluvial flood-susceptible area. In the Maxent modelling, 70% of the flood point records (49 points) were randomly selected and used for the training set, and the remaining 30% (21 points) were used for the validation set in assessing the Maxent result.

| Validation of the AHP method
For the validation of the resulting pluvial flood-prone areas by the AHP method, a GPS-based field survey of local people was carried out along the east coast of Davao Oriental to investigate true ground points flooded by historical flooding events. From the field survey, coordinates of 70 true ground points were collected and used for performance evaluation (see Figure 1d for their locations).
The performance evaluation for the AHP method was performed 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 Equations (13)-(17), respectively (Feng, Liu, & Gong, 2015; where P, N, TP, TN, FP, and FN are condition positive, condition negative, true positive, true negative, false positive, and false negative, respectively. In the accuracy assessment of the pluvial flood-prone map, a confusion matrix was calculated based on the field survey of ground points for historical flooding events. However, the field survey was carried out along the east coast of Davao Oriental, and the Davao Gulf municipalities were not included in the survey. Thus, the points do not cover all barangays in the study area (Figure 1d). To solve this limitation, the elevations of 70 points were extracted from the DEM. The elevations of the 70 points were in the range between 10 and 20 m above mean sea level. Then, the points were spatially interpolated, and the areas below 20 m elevations where all points are located were used to determine the pluvial flood and non-flood areas in Davao Oriental. The result is shown in Figure 1c. Finally, the accuracy assessment was performed using the confusion matrix.

| Validation of the Maxent model
As described, the Maxent model was trained with 70% of the field survey points, and the result was validated with the remaining 30% of the survey points. The accuracy assessment for the Maxent model in the validation was performed using the area under the curve (AUC). The AUC is a popular, comprehensive quantitative method of accuracy assessment that can be used to evaluate the prediction and success rates (Tehrany et al., 2017).
The process of validation was executed by a comparison of known data on flooding with the probability map of acquired pluvial flooding using AUC (Bui, Pradhan, Lofman, Revhaug, & Dick, 2012). The values AUC = 1 and AUC = 0.5 demonstrate a perfect classification and a classification by chance, respectively. Some studies have used AUC in evaluating accuracy (Bui et al., 2012;Pourghasemi, Pradhan, Gokceoglu, Pourghasemi, & Gokceoglu, 2012;Tehrany et al., 2017). The technique involves dividing the probability map into equal area groups and ranking these values in a minimum-tomaximum hierarchy. The percentage of pluvial flooding occurrence for each probability category is determined from the success and prediction curves. These curves are formed by plotting percentages of pluvial floodsusceptible areas from highest to lowest on the x-axis and the percentages of pluvial flood events on the y-axis. A steeper curve indicates a greater number of pluvial flood events falling into categories of greater susceptibility (Tehrany et al., 2017).

| Vulnerability map
Vulnerability expresses the level of inability to resist a hazard or to respond when disaster has occurred. 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 (Danumah et al., 2016).
The susceptibility and exposure issues include several factors, such as the age and health of the population, socio-economic activities, the quality of buildings and their location with respect to flooding. In this study, the only indicator used for the vulnerability to pluvial flooding is the population density. The data were recorded in the attribute table and available at the barangay level. The barangay boundary is a shapefile and converted into a raster with a 30 × 30 m grid resolution. The vulnerability map obtained from the population density, as shown in Figure 4f, is divided into 10 classifications. Barangays with high populations are mostly proximal to coastal areas.

| AHP flood-prone hazard map
The six criteria for pluvial flood hazards are presented in Figure 4. The resulting hazard maps for the AHP are shown in Figure 4g. The AHP flood-prone map depicts in the ground point accuracy assessment that they have a good accuracy rate of 81%, as shown in Table 4. In the overall assessment, the ACC was 63%, the TNR was 81%, and the TPR was 43%, as shown in Table 5. Therefore, the analysis reveals that the AHP has high accuracy in predicting flood-prone areas.
Finally, in the AHP, the very low and low classes cover 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, rainfall (42.2%) and slope (22.5%) are the most significant causative factors of pluvial flood occurrence. In the flood-prone area, 64 out of 183 barangays (towns) and 30% (168,034/558,958) of the population are at a high risk of pluvial flooding in the current situation, as shown in Table 6. Additionally, these areas have high occurrences of rainfall, high drainage densities, low elevations, and low slopes, and they are proximal to river channels and shorelines.

| Maxent flood-prone hazard map
The resulting pluvial flood hazard map by Maxent modelling is presented in Figure 4h with the same five classifications as those of the AHP flood hazard map. In the simulation of the flooding using the Maxent model, the six criteria together with 70% of the survey ground points of the historical events (see Figure 1c) were used to generate the flooding.
The percent contribution of every criterion is 37. 7, 22.4, 20.8, 14.1, 4.3, and 0.7% for the elevation, rainfall, soil, slope, drainage density, and distance to the main channel, respectively. The results show a very high accuracy with an AUC of 0.965. The classifications of the map show approximately 90.39, 4.36, 2.72, 1.4, and 1.13% spatial distribution of flooding for very low (VL), low (L), moderate (M), high (H), and very high (VH), respectively. The VL,L,M,H,and VH indicate 20,40,60,80, and 100% of the chances of flooding. Although the classifications of the M, H, and VH are approximately only 6%, these areas are on the coastline and riverside, which are prone to the risk of flooding. Additionally, those areas are in the east coast municipalities of Davao Oriental.
The Maxent model highly emphasises elevation, soil, slope, and rainfall. Moreover, the AUC is 96.5% compared to the AHP accuracy of 63%, as described in the T A B L E 4 Accuracy assessment using the ground truth points from the field survey

No.
Flood-prone class Verification points previous section. The high AUC close to 1 indicates high accuracy. The crucial part in the analysis or modelling of the flood-prone area using Maxent is the dataset of the historical events of flooding. Due to the limited field survey for the historical flood events along the east coast of Davao Oriental, the information for the municipalities in the Davao Gulf are not considered, resulting in the high value of 90.39% for very low classification over the areas without survey points ( Figure 4h). Therefore, a direct comparison between the two methods for the whole study area is yet to be made.

| Pluvial flood risk map
Because of the restriction of the Maxent model due to the limited survey points, the pluvial flood risk map was evaluated using the AHP flood hazard map. The pluvial flood risk map was generated by a weighted overlay of the hazard and vulnerability maps with equal weights with five classifications ( Figure 5). The very low, low, and moderate classes cover 0.63, 60.64, and 35.35% of the total area, respectively. The categories of high and very high are estimated at 3.34 and 0.05% of the total area, respectively. These areas are barangays with high population densities and are mostly in low-lying areas near the shoreline and rivers. The moderate and very low classes are areas at high elevations with smaller populations. Thirty-one out of one-hundred eighty-three (31/183) barangays are at high risk of pluvial flooding. Only one barangay has a very high risk of pluvial flooding. Therefore, that barangay requires the implementation of a strong intervention policy that would decrease the number of casualties during pluvial flooding.

| DISCUSSION
Interestingly, the weights for criteria in the AHP method are 42% for rainfall, 23% for slope, 15% for elevation, 10% for distance to the main channel, 6% for drainage, and 4% for soil type, whereas the contribution rates for criteria in the Maxent model are 37.7% for elevation, 22.4% for rainfall, 20.8% for soil type, 14.1% for slope, 4.3% for drainage, and 0.7% for distance to the main channel. In this pluvial flood hazard and risk assessment, the highest weight is given to the rainfall, and the elevation is ranked third in weight in the AHP method. The soil type is the least important factor in the AHP. However, elevation is the most important criterion, and rainfall is the second most important criterion in the Maxent model. The soil type becomes very important in Maxent, with a contribution rate of 14.1%. Overall, the relative importance of each criterion is very different, varying with the method applied.
In the AHP, data from different sources with different resolutions were factors of bias during processing and analysis. Conversely, the addition of weights reduced the bias and uncertainty in the result. However, the AHP, as one of the MCDA methods, has some limitations due to the subjectivity in choosing the weights of criteria from the arbitrary judgements of experts (Danumah et al., 2016;Papaioannou et al., 2015). Fortunately, this subjectivity problem can be reduced by the consistency ratio test of the judgements. To make a coherent judgement, Saaty (1980) provided a consistency ratio threshold that must be less than 0.10 (10%).
Unfortunately, the range of the CR between 0 and 9.9% is a rough estimation. The HI is coherent if the CR is 9.9% or the CR is 0%, but the weights of the indicators vary in every CR generation. In Cabrera and Lee (2019), the sensitivity of HI was tested in three different CR scenarios: (a) the lowest CR is 2.0%, (b) the highest CR is 8.6%, and (c) the middle CR is 5.8%. Their results demonstrate that the average change in the ranking of flood hazard classifications according to the CR values varies from −0.77% for the H classification to 1.21% for the M classification, and the overall average change is 0.02%. Therefore, it is found that the changes in the CR did not significantly affect the results of the hazard map, and the changes are acceptable. The consistency ratio of this T A B L E 6 Flood risk assessment represented by the affected population and the number of affected barangays in each municipality by the AHP method study is 0.02 such that the judgements and the HI can be considered coherent. The Maxent model learns by combining the sample layer (the layer for the flood points from the field survey) and the environmental layers (rainfall, slope, elevation, drainage density, and distance to the main channel) and shows the contribution rates of the environmental layers to pluvial flooding. In this study, we found that the two models can complement each other. For instance, in an observation data-scarce area, the Maxent model can be applied to determine the contribution rates of the criteria used. Then, those contribution rates can be utilised as the weights for criteria in the AHP method. By applying this combined Maxent and AHP approach, the subjectivity in the weighting criteria due to the judgements of experts can be reduced in the AHP while assessing the pluvial flood hazard and risk with improved accuracy in the areas with a lack of observations and historical records.

| CONCLUSIONS
Flooding is a natural hazard that poses a risk to both people and physical structures within the affected areas. Flood risk occurrence can be better understood through the use of a spatial assessment. Thus, in this study, the flood-prone areas and flood risk due to pluvial flooding for Davao Oriental in Mindanao of the Philippines are F I G U R E 5 Pluvial flood risk map by the AHP method assessed by a GIS-based MCDA based on the AHP method and Maxent model with seven criteria, such as elevation, rainfall, soil type, slope, drainage density, distance to the main channel, and population density.
The 55 and 22% of the total area of Davao Oriental are classified into moderate and high classes of pluvial flooding, respectively, and approximately 30% of the total population appears to be at high risk of pluvial flooding in the current situation by the AHP method. The rainfall (42.2%) and slope (22.5%) are the most significant causative factors of pluvial flood occurrence in the AHP. In the Maxent model, elevation (37.7%) and rainfall (22.4%) are the most important factors.
In the validation of pluvial flood-prone area assessments, the Maxent model results show better accuracy over the limited area with historical flood records available, particularly in addressing the relationship between pluvial flood occurrences and the environmental factors (criteria), compared to that from the AHP. In contrast, the AHP flood-prone hazard map is appropriate in datascarce areas such as Davao Oriental of this study due to the limited historical records of flooding events.
In most cases, developing countries lack data measurements, which will lead to inappropriateness in predicting future events such as flooding. Successful flood risk assessment requires an understanding of the quantity and quality of observed events. The AHP and Maxent approaches and their combination for pluvial flood hazard and risk assessment can reduce the subjectivity and uncertainty in selecting and weighting criteria due to the judgements of experts in the AHP while improving the accuracies in the assessment results for data-scarce areas in developing countries. Finally, this study, based on a broad-scale high-level data-driven screening method, can help identify potential hot spots for pluvial flooding, for which more detailed numerical modelling studies should be undertaken. It cannot be used in place of hydrodynamic modelling.