Comprehensive Assessment of Flood Hazard, Vulnerability, and Flood Risk at the Household Level in a Municipality Area: A Case Study of Nan Province, Thailand

: Estimating ﬂood hazard, vulnerability, and ﬂood risk at the household level in the past did not fully consider all relevant parameters. The main objective of this study is to improve this drawback by developing a new comprehensive and systematic methodology considering all relevant parameters and their weighting factors. This new methodology is applied to a case study of ﬂood inundation in a municipal area of Nan City in the Upper Nan River Basin in Thailand. Field and questionnaire surveys were carried out to collect pertinent data for input into the new methodology for estimating ﬂood hazard, vulnerability, and risk. Designed ﬂoods for various return periods were predicted using ﬂood simulation models for assessing ﬂood risk. The ﬂood risk maps constructed for the return periods of 10–500 years show a substantial increase in ﬂood risk with the return periods. The results are consistent with past ﬂood damages, which were signiﬁcant near and along the riverbanks where ground elevation is low, population density is high, and the number of household properties are high. In conclusion, this new comprehensive methodology yielded realistic results and can be used further to assess the effectiveness of various proposed ﬂood mitigation measures.


Introduction
Floods are the most frequent type of disasters. They can cause widespread devastation, loss of life, and damage to personal property and public infrastructure. Between 1998 and 2017, floods affected more than 2 billion people, or 29% of the world population [1]. Floods cause more than USD 40 billion in damage worldwide annually [2]. The number of deaths due to the fact of significant floods worldwide from 1906 to 2014 was reported by Ritchie and Roser [3] (updated in 2020) as shown in Figure 1. In 1931, deaths were as high as 3.7 million in China, while in 1974, there were 28,700 deaths in Bangladesh [4].
People who live in floodplains or lack information on flood warning and flood hazard awareness are the most vulnerable to floods. Floods are expected to increase in frequency and intensity due to the anticipated increase in extreme precipitation as a result of continued climate change [5].
Many people tend to live in floodplain areas to satisfy needs for water and transportation. However, living in floodplain areas increases the risk to their assets, household, and life security during flooding [6][7][8]. A flood hazard, vulnerability, and risk assessment is necessary to determine the probability that a flood of a particular intensity may occur over an extended period of time [9]. A risk assessment aims to estimate this flooding probability from over a period of years to decades to support flood mitigation planning and management activities [10]. Extreme flood events can be more harmful to people living in the floodplain areas [11]. Another important consequence of floods is damage to physical infrastructures such as bridges, main roads, public buildings, and service facilities [12][13][14]. This type of structural damage is often defined as the amount of money needed to repair or rebuild. On the other hand, damages to populations and their household properties are related to socioeconomic conditions [15,16]. From flood hazard and flood vulnerability, flood risk can be determined, and flood risk maps can be constructed for various return periods for flood mitigation planning and management. A flood risk map can be helpful for land-use planning and management, reducing flood damage and saving human lives, agricultural products, and private and public properties.
Previous studies [17][18][19][20][21][22][23] considered flood hazard as a function of flood depth only or flood depth and velocity. This study considered all parameters including flood depth, velocity, inundation duration, and their weights. For the total vulnerability at the household level, this study considered all parameters including the sensitivity of the population, adaptive capacity, exposure of household properties, and their weights. These parameters have not been fully considered in the past.
The main aim of this study was to develop a comprehensive and systematic methodology for determining flood risks considering all hazard and vulnerability parameters and their weights using flood simulation models, questionnaire surveys, and the analytic hierarchy approach (AHP) [24,25]. The computational framework is given in Section 3 on the research structure and methodology.
The comprehensive and systematic methodology developed in this study is the main contribution, as it will be an essential foundation for properly assessing flood hazard, vulnerability, and flood risk in the future.
Hence, the objectives of this study were as follows: (a) To simulate floods in the Upper Nan River and its floodplain in the municipal area of Nan Province; (b) To develop a comprehensive and systematic methodology to determine flood hazard, flood vulnerability, and flood risk in a municipal area at the household level; (c) To apply the developed methodology to assess flood hazard, flood vulnerability, and flood risk at the household level in the Nan Municipality area considering floods of various return periods; Extreme flood events can be more harmful to people living in the floodplain areas [11]. Another important consequence of floods is damage to physical infrastructures such as bridges, main roads, public buildings, and service facilities [12][13][14]. This type of structural damage is often defined as the amount of money needed to repair or rebuild. On the other hand, damages to populations and their household properties are related to socioeconomic conditions [15,16]. From flood hazard and flood vulnerability, flood risk can be determined, and flood risk maps can be constructed for various return periods for flood mitigation planning and management. A flood risk map can be helpful for land-use planning and management, reducing flood damage and saving human lives, agricultural products, and private and public properties.
Previous studies [17][18][19][20][21][22][23] considered flood hazard as a function of flood depth only or flood depth and velocity. This study considered all parameters including flood depth, velocity, inundation duration, and their weights. For the total vulnerability at the household level, this study considered all parameters including the sensitivity of the population, adaptive capacity, exposure of household properties, and their weights. These parameters have not been fully considered in the past.
The main aim of this study was to develop a comprehensive and systematic methodology for determining flood risks considering all hazard and vulnerability parameters and their weights using flood simulation models, questionnaire surveys, and the analytic hierarchy approach (AHP) [24,25]. The computational framework is given in Section 3 on the research structure and methodology.
The comprehensive and systematic methodology developed in this study is the main contribution, as it will be an essential foundation for properly assessing flood hazard, vulnerability, and flood risk in the future.
Hence, the objectives of this study were as follows: (a) To simulate floods in the Upper Nan River and its floodplain in the municipal area of Nan Province; (b) To develop a comprehensive and systematic methodology to determine flood hazard, flood vulnerability, and flood risk in a municipal area at the household level; (c) To apply the developed methodology to assess flood hazard, flood vulnerability, and flood risk at the household level in the Nan Municipality area considering floods of various return periods; (d) To analyze and discuss the consequences of flood hazard, vulnerability, and flood risk on local residents and the physical environment.

Study Area
The Upper Nan River Basin in Northern Thailand has a drainage area of 8709.50 km 2 [26]. Its main drainage way is the Nan River. It has an annual average temperature of 26.3 • C, relative humidity of 75.9%, wind speed of 0.9 knots, evaporation of 1457.4 mm, and rainfall of 1371 mm.
The basin drainage area from the river gauging station (N1) at Nan Province up to its headwater boundary of the river basin has an area of 5663 km 2 [27]. Another river gauging station upstream of station N1 is the Tha Wang Pha station (N64). The basin study area is subdivided into five sub-basins, namely, Upper Nan, Nam Yao (W), Nam Yao (E), Nam Samun, and Nam Nan part-2 as shown in Figure 2 [28]. (d) To analyze and discuss the consequences of flood hazard, vulnerability, and flood risk on local residents and the physical environment.

Study Area
The Upper Nan River Basin in Northern Thailand has a drainage area of 8709.50 km 2 [26]. Its main drainage way is the Nan River. It has an annual average temperature of 26.3 °C, relative humidity of 75.9%, wind speed of 0.9 knots, evaporation of 1457.4 mm, and rainfall of 1371 mm.
The basin drainage area from the river gauging station (N1) at Nan Province up to its headwater boundary of the river basin has an area of 5663 km 2 [27]. Another river gauging station upstream of station N1 is the Tha Wang Pha station (N64). The basin study area is subdivided into five sub-basins, namely, Upper Nan, Nam Yao (W), Nam Yao (E), Nam Samun, and Nam Nan part-2 as shown in Figure 2 [28]. These five sub-basins contribute runoff to the Nan River which flows from north to south through Nan's municipal area of 7.6 km 2 on the right bank of the Nan River as shown in Figure 3. The municipal area is an urban area that has 31 villages that are mainly residential. The Nan River flows through the municipal area for a distance of 7 km. In this river reach, overbank flow occurs frequently during high flood periods. The municipal area has flooded often in the past. Large flood inundations occurred in 2006 [29] and 2011 due to the fact of heavy tropical storms [30]. Considerable flood damages and casualties were reported in both years. These five sub-basins contribute runoff to the Nan River which flows from north to south through Nan's municipal area of 7.6 km 2 on the right bank of the Nan River as shown in Figure 3. The municipal area is an urban area that has 31 villages that are mainly residential. The Nan River flows through the municipal area for a distance of 7 km. In this river reach, overbank flow occurs frequently during high flood periods. The municipal area has flooded often in the past. Large flood inundations occurred in 2006 [29] and 2011 due to the fact of heavy tropical storms [30]. Considerable flood damages and casualties were reported in both years.

Research Structure
The research structure and methodology of this study is well described in Figure 4. The research structure has two parallel paths: Path A for flood hazard estimation and Path B for vulnerability estimation. The flow of work in Path A is from Block A.1 to Block A.4 as shown in the figure. The flow of work in Path B is from Block B.1 to Block B.9. The output of Block A.1, or flood hazard, and the output of Block B.9, or vulnerability, are multiplied to obtain the flood risk as shown in Block C near the bottom of the figure. The flood risk maps are constructed, as shown in Block D, based on the computed flood risk in Block C. From here onward, flood risk analysis and assessment were carried out to obtain a risk reference for further evaluation of risk reduction of proposed flood control

Research Structure
The research structure and methodology of this study is well described in Figure 4. The research structure has two parallel paths: Path A for flood hazard estimation and Path B for vulnerability estimation. The flow of work in Path A is from Block A.1 to Block A.4 as shown in the figure. The flow of work in Path B is from Block B.1 to Block B.9. The output of Block A.1, or flood hazard, and the output of Block B.9, or vulnerability, are multiplied to obtain the flood risk as shown in Block C near the bottom of the figure. The flood risk maps are constructed, as shown in Block D, based on the computed flood risk in Block C. From here onward, flood risk analysis and assessment were carried out to obtain a risk reference for further evaluation of risk reduction of proposed flood control measures. The methodology for the computational flood hazard, vulnerability, and flood risk analysis is described in detail in the following section.

Determination of Flood Hazard
The flood hazard index (FHI) is expressed by Equation (1) in which FHIT is the flood duration index, FHID is the flood depth index, and FHIV is the flood flow velocity index. The symbols α, β, and µ are their weighting factors, respectively.
The values for FHIT, FHID, and FHIV are determined according to the specified ranges of the computed flood duration, flood depth, and flood velocity in the study area as shown in Table 1. The computed flood duration is classified into five ranges: very low, low, medium, high, and very high as shown in Table 1. The corresponding FHIT of each class is specified as 1, 2, 3, 4, and 5, respectively. The same method is applied to the FHID and FHIV. The sum of the weights α, β, and µ is equal 1. The weights α, β, and µ are determined using the analytic hierarchy process (AHP) of multi-criteria decision analysis [24,25,[31][32][33] based on field and questionnaire surveys. The details of determining α, β, and µ by using AHP is explained in the Appendix A. Finally, the hazard level of FHI is specified as very low when FHI ≤ 1, low when 1 < FHI ≤ 2, medium when 2 < FHI ≤ 3, high when 3 < FHI ≤ 4, and very high when 4 < FHI ≤ 5. The total flood damage vulnerability index (FVI) for each village is calculated as the sum of the weighted flood vulnerability index of the population, FVI pop , and the weighted flood vulnerability index of households, FVI hh . The total vulnerability index (FVI) of each village is calculated by the following equations: The weights w pop and w hh are determined using AHP. The level of total vulnerability FVI is specified as very low when FVI ≤ 1, low when 1 < FVI ≤ 2, medium when 2 < FVI ≤ 3, high when 3 < FVI ≤ 4, and very high when 4 < FVI ≤ 5.
In determining FVI pop , population vulnerability, VI pop , is represented by population density (person/km 2 ) and is classified into five ranges in ascending order as given in Table 2. These five ranges are represented by the integer numbers 1, 2, 3, 4, and 5 for the index FVI pop : very low, low, medium, high, and very high respectively. The household vulnerability index, FVI hh , is determined by the value of the household vulnerability. VI hh and its ranges according to Table 2. The VI hh of the household samples is expressed as a function of three major contributing factors, namely: sensitivity F1, adaptive capacity F2, and exposure F3 according to the following equation: where w 1 , w 2 , and w 3 are the weights of F1, F2, and F3, respectively. The major contributing factors F1, F2, and F3 are defined in [24,25]. Equation (5) was developed in this study based on the original equation in [34][35][36] by including the weights w 1 , w 2 , and w 3 in Equation (5). The weighting factors w 1 , w 2 , and w 3 were introduced to normalize the relative importance of F1, F2, and F3. The values of w 1 , w 2 , and w 3 were determined by the AHP using data from collected samples of questionnaires and field surveys. Equation (6) requires that the sum of the weighting factors w 1 , w 2 and w 3 is equal to one. Each major contributing factor is divided into various components in which each component is further classified into sub-classes. An impact score is assigned to each class, as shown in Table 3, according to the influence or impact of the class on its component.  Table 3. Cont. (5) and (6) Weights W i of F i in Equations (5) and (6) Components C i of Factor F i in Equation (7) & Definitions

Major Con-tributing Factors F i in Equations
Classification of Component C i in Equation (7)  θ 1 = 0.408 Area at shorter distance to river is more exposed to flooding For example, the sensitivity F1 is further classified into 7 contributing components C 1 -7 in which C 1 = family size, C 2 = gender, C 3 = health, C 4 = land use type, C 5 = household damage, C 6 = public property damages, and C 7 = ownership of household. The sensitivity F1 is computed as: where C i is the ith contributing component of the sensitivity factor F1, and n = 7 is the total number of the component C i of F1.
Each contributing component C i of F1 is computed as: where the weight θ i of C i is determined by AHP based on the collected samples from questionnaires and field surveys; m is the number of classes of the collected samples for the component C i , arranged from the highest to the lowest significance of vulnerability; Q j is the number of samples of class j as a percentage of the total collected samples of all classes; K j is the assigned impact score of the class j between 0% and 100%. A score of 100% is assigned to the class that has the highest impact on vulnerability. More details regarding the calculation are given in the next section on the computational procedure. The component C 1 for F1 is calculated as follows: where Q j is the number of family samples in the jth class with the assigned impact score K j , and θ 1 is the weight of C 1 to be determined by AHP. The same overall procedure described above is applied to the other two major contributing factors, namely: adaptation F2 and exposure F3. For the adaptation capacity F2, there are seven components as shown in Table 3, namely: education level of householders, their employment, household income, saving deposit, insurance for flood events, land price, and notification of flood events.
The exposure F3 in Table 3 is composed of six major components, namely: the distance from the river, ground elevation, the inundation depth in 2011, the flood velocity in 2011, the duration of inundation in 2011, and the number of flooding events in 2011. The year 2011 was the only year that flood inundation data in the study area were available. The exposure factor F3 was calculated based on the sum of its weighted components in which the weights were determined using AHP.
After obtaining the values of the three major contributing factors, namely, sensitivity F1, adaptation F2, and exposure F3, they were substituted into Equation (5) to calculate household vulnerability, VI hh . The computed household vulnerability, VI hh , of each village was fitted into the five classified ranges: very low, low, medium, high, and very high as shown in Table 2. Then, the value of FVI hh corresponding to VI hh was represented by an integer number from 1 to 5 according to Table 2.
Finally, the vulnerability indices of population density, FVI pop , and of household FVI hh were substituted into Equation (3) to compute the total vulnerability, VI.
Each village area occupies many grids of 30 by 30 m, which are used in calculating the flood hazard. As the grid size is relatively much smaller than the area of each village, all the grids in the same village are considered to have the same value of the total vulnerability, FVI.

Determination of Flood Risk
By using the same grids for hazard and vulnerability, the flood risk index, FRI, for each grid in the study area was computed as the product of the FHI and the total flood damage vulnerability index, FVI. The flood risk index (FRI) is calculated as: Equation (11) does not explicitly include exposure. This is different from the equation given in [16], which includes exposure. In this study, exposure data were collected from the questionnaires and field surveys as shown in Table 3. This exposure term, F3, was considered in Equation (5) when calculating household vulnerability, V hh , and the total vulnerability, FVI. The computed values of the FRI were classified into five ranges corresponding to very low for 1 < FRI ≤ 5, low for 5 < FRI ≤10, medium for 10 < FRI ≤ 15, high for 15 < FRI ≤ 20, and very high for 20 < FRI ≤ 25.

Hydrological Data
For the hydrological model simulation and the hydrodynamic model simulation, the hourly hydro-meteorological data from four stations (i.e., 331201, 331301, 331401, and 331402), as shown in Figure 2, and the daily river discharge and water level data at two stations at N1 (Nan) and N64 (Tha Wang Pha) were collected from the regional offices of the Thai Meteorological Department and of the Royal Irrigation Department (RID), Thailand. The river cross-sections, floodplain geometry, land topography, and DEM were collected from the Nan Provincial Office, the Regional Office of the Royal Irrigation Department (RID), and from the Geo-Informatics and Space Technology Development Agency (GISTDA).

Vulnerability Data
The following data were collected for the analysis: (a) land-use maps in the Nan River Basin from the Land Development Department (LDD); (b) population data and household flood damage and public infrastructure damage date in Nan Municipality from the Nan Provincial Administration; (c) socioeconomic and household data were collected from field and questionnaire surveys in the study area. The number of household samples collected for the questionnaire survey was estimated using Yamane's equation [37]: where n is the number of household samples; N is the total number of households in the study area, which was 6310; e is the allowable error taken equal to 10%. The value of n was calculated and found to be approximately 100.

Computation of Flood Hazard
The HEC-HMS rainfall-runoff model and the HEC-RAS flood routing model are connected as shown in Figure 5 in computing flood discharges and depths. The computational procedure is as follows: Figure 5. Configuration of the HEC-HMS and HEC-RAS models for the upstream sub-basins of the Nan River and the Nan municipal area.
(a) Rainfall-runoff computation: The HEC-HMS rainfall-runoff model [38] was applied to compute the runoff hydrograph using hourly rainfall input at four stations in the Upper Nan Basin. The hourly rainfalls at the four stations were averaged over the basin area using the Thiessen polygon method. The river basin was divided into seven sub-basins in which the hourly average rainfalls were used in each sub-basin.
The computed runoff was used as the upstream boundary condition of the HEC-RAS flood routing model [39]. The HEC-HMS model requires a digital elevation model (DEM), soil and land-use maps, soil characteristics, and input rainfall hyetographs. The HEC-Geo HMS model, which is an extension of HEC-HMS, prepares raster layers of delineated sub-basins and river network systems for exporting to HEC-HMS as base maps. By inputting rainfall data, land cover, and soil maps to HEC-HMS, the model computes daily runoff hydrographs for each sub-basin. The HEC-HMS model was calibrated and verified against the observed daily discharges at station N64 at Tha Wang Pha and at station N1 at Muang Nan. The calibration period was during the wet period from June to December 2006-2011, and the verification period was from June to December 2012-2017. In the model calibration, the model parameters, such as initial and maximum storages of canopy, SCS curve number, time of concentration, and lag of unit hydrograph, were assumed and adjusted by trial and error to obtain a satisfactory agreement between the observed and computed discharge hydrographs; (b) Flood routing computation: The 1D and 2D HEC-RAS flood routing model [39] for the Upper Nan River and its floodplain were used to route the runoff from the upstream station N1 along the Nan River to the downstream end station, which is 7 km downstream of station N1. The river passes through the municipal area, which is in the river floodplain. The geometrical inputs to HEC-RAS were the measured river cross-sections every 1.2 km and the floodplain topography from the digital elevation model with a 30 by 30 m resolution with a 1 m contour interval. The river crosssections from the field measurements and the floodplain topography from the DEM were merged using HEC-GEO RAS, an extension of HEC-RAS, to obtain the complete river and floodplain cross-sections [40]. This geometry data was input into HEC-RAS (a) Rainfall-runoff computation: The HEC-HMS rainfall-runoff model [38] was applied to compute the runoff hydrograph using hourly rainfall input at four stations in the Upper Nan Basin. The hourly rainfalls at the four stations were averaged over the basin area using the Thiessen polygon method. The river basin was divided into seven sub-basins in which the hourly average rainfalls were used in each sub-basin.
The computed runoff was used as the upstream boundary condition of the HEC-RAS flood routing model [39]. The HEC-HMS model requires a digital elevation model (DEM), soil and land-use maps, soil characteristics, and input rainfall hyetographs. The HEC-Geo HMS model, which is an extension of HEC-HMS, prepares raster layers of delineated sub-basins and river network systems for exporting to HEC-HMS as base maps. By inputting rainfall data, land cover, and soil maps to HEC-HMS, the model computes daily runoff hydrographs for each sub-basin. The HEC-HMS model was calibrated and verified against the observed daily discharges at station N64 at Tha Wang Pha and at station N1 at Muang Nan. The calibration period was during the wet period from June to December 2006-2011, and the verification period was from June to December 2012-2017. In the model calibration, the model parameters, such as initial and maximum storages of canopy, SCS curve number, time of concentration, and lag of unit hydrograph, were assumed and adjusted by trial and error to obtain a satisfactory agreement between the observed and computed discharge hydrographs; (b) Flood routing computation: The 1D and 2D HEC-RAS flood routing model [39] for the Upper Nan River and its floodplain were used to route the runoff from the upstream station N1 along the Nan River to the downstream end station, which is 7 km downstream of station N1. The river passes through the municipal area, which is in the river floodplain. The geometrical inputs to HEC-RAS were the measured river cross-sections every 1.2 km and the floodplain topography from the digital elevation model with a 30 by 30 m resolution with a 1 m contour interval. The river cross-sections from the field measurements and the floodplain topography from the DEM were merged using HEC-GEO RAS, an extension of HEC-RAS, to obtain the complete river and floodplain cross-sections [40]. This geometry data was input into HEC-RAS for flood flow simulation. In the HEC-RAS model, the 1D flow routing procedure was used to compute the 1D flow in the Nan River, and the 2D flow routing was used to compute the 2D flow in the floodplain; The hydrological inputs to HEC-RAS were the observed daily upstream discharge at station N1 as the upstream boundary condition. At the downstream end of the model, there was no river gauging station; thus, a depth-discharge relationship according to the Manning equation was used. The model was calibrated and verified by trial-and-error adjustment based on the values of the Manning roughness coefficient n.  Table 3. As shown in Figure 4, for Blocks A.1 to A.4, the computed FHI D , FHI V , and FHI T were substituted into Equation (1) to compute the FHI. The hazard weighting factors α for flood duration, β for flood depth, and µ for flood velocity in Equation (1) were determined using AHP [24,25].

Computation of Total Vulnerability
As shown in Figure 4, for Blocks B.1 to B.8, the total vulnerability of the study area was computed as a function of the vulnerability of the population and the vulnerability of households. According to Equations (3) and (4), the total vulnerability is equal to the sum of the weighted vulnerability of the population and the weighted vulnerability of households. The weights of the vulnerability of the population and of household vulnerability were determined using AHP based on the questionnaire surveys.
The vulnerability of the population of each village was computed from the population density (number of persons/km 2 ) of each village from census data. For the household vulnerability, VI hh , each major contributing factor of VI hh was composed of a number of components in which each component was further subdivided into classes with assigned impact scores. For example, the major contributing factor F1, or sensitivity, had seven components (C 1 -7 ), namely: family size, gender, health, land-use type, household damage, public damage, and household ownership as shown in Table 3. An example of calculating the component C 1 or family size of F1 in Table 3 for Phumin-Thali village is given here. The family size was subdivided into three classes: class 1 for households with more than five family members; class 2 for between three and five family members; class 3 for fewer than three family members. The impact scores of 100%, 67%, and 33% were assigned to each family class according to the information from local residents on vulnerability to their families. A score of 100% was assigned to class 1, because it is most sensitive to vulnerability. The number of collected household samples in Phumin-Thali village was six.
The component C 1 of the major contributing factor F1 was classified into three classes (m = 3). In class 1 (j = 1), the number of the families was equal to one. Hence, Q 1 was equal to 16.67% of the total collected samples of six, and the impact score K 1 was assigned to be 100%. In class 2 (j = 2), the number of families was equal to two or Q 2 = 33.33% and assigned K 2 = 67%. In class 3 (j = 3), the number of families was three or Q 3 = 50% and K 3 = 33%. The weighting factor θ 1 of the sensitivity factor F1 was determined by AHP and equal to 0.386. This is shown at the beginning of Table 3. Knowing Q 1 , Q 2 , Q 3 , K 1 , K 2 , K 3 , and θ 1 , the value of the component C 1 of F1 is equal to θ 1 (Q 1 K 1 + Q 2 K 2 + Q 3 K 3 )/100 = 21.41%.
The same procedure was repeated for the remaining classes 2 and 3 of the component C 1 of the factor F1 as shown in Table 3. In this way, the components C 1 -7 were calculated using Equation (8) and summed to obtain the major contributing factor F1 according to Equation (7).
The same overall procedure was repeated for the major contributing factors F2 and F3. Hence, the total vulnerability, VI, and its index, FVI, was determined from Equations (3) and (4).

Computation of Flood Risk
As shown in Figure 4, the FRI was computed as the product of the FHI and the FVI. The distribution of flood risk over the municipal area can be represented by a flood risk map of grids of 30 by 30 m resolution. The maps of flood hazard and risk of a spatial resolution of 30 × 30 m were constructed for 10, 50, 100, and 500 year return periods. The flood risk index of each grid for each return period was classified into five levels of equal intervals, namely: very low for 1 < FRI ≤ 5, low for 5 < FRI ≤ 10, medium for 10 < FRI ≤ 15, high for 15 < FRI ≤ 20, and very high for 20 < FRI ≤ 25.

Calibration and Verification of the HEC-HMS Rainfall-Runoff Model and HEC-RAS Flood Routing Model
The HEC-HMS rainfall-runoff model, which covers the Upper Nan River Basin area from its headwater to the model outlet at the gauging station N1, was calibrated and verified against the observed daily discharges at gauging stations N64 at Tha Wang Pha and N1 at Mueang Nan. The calibration period was during the wet period from June to December in 2006-2011, and the verification period was from June to December in 2012-2017. In the model's calibration, the model parameters, such as initial and maximum storages of canopy, SCS curve number, time of concentration, and lag of unit hydrograph, were assumed based on previous studies [21,40] and adjusted by trial and error to obtain a satisfactory agreement between the observed and computed discharge hydrographs ( Figure 6). The model's performance was evaluated using the following statistical parameters, namely: coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (NSE), percent bias (PBIAS), volume ratio (Vr), and normalized root mean square error (NRMSE). The results of the model's calibration at stations N64 and N1 were found to be satisfactory with coefficients of determination of 0.72 and 0.76; Nash-Sutcliffe efficiencies of 0.67 and 0.64; percent biases of 15.78 and 15.11; volume ratios of 0.95 and 0.89; normalized root mean square errors of 12.29% and 14.36%, respectively. The statistics showed acceptable agreement in both the model's calibration and verification.
The HEC-RAS flood routing model covers the Nan River reach from the river gauging station N1 to the downstream end station 7 km downstream and its floodplain in the municipal area. The model was calibrated by comparing the computed and observed daily water levels at station N1 during the peak flood periods in 2006-2011, while in the floodplain, there was no water level data. According to [41], the Manning roughness coefficient n of 0.03 was initially assumed for the river and its floodplain in the calibration. The assumed Manning roughness coefficients were adjusted by trial and error and were found to be 0.032 for the river and 0.035 for the floodplain. These Manning n values matched with the roughness condition of the river channel and its floodplain, which is moderately rough according to [9].  The HEC-RAS flood routing model covers the Nan River reach from the river gauging station N1 to the downstream end station 7 km downstream and its floodplain in the municipal area. The model was calibrated by comparing the computed and observed daily water levels at station N1 during the peak flood periods in 2006-2011, while in the floodplain, there was no water level data. According to [41], the Manning roughness coefficient n of 0.03 was initially assumed for the river and its floodplain in the calibration. The assumed Manning roughness coefficients were adjusted by trial and error and were found to be 0.032 for the river and 0.035 for the floodplain. These Manning n values matched with the roughness condition of the river channel and its floodplain, which is moderately rough according to [9].  The results of the calibration and verification of the HEC-RAS model for the river flow during flood peak periods are shown in Figure 7. As shown in the results, the performance statistics showed acceptable agreement for both the model's calibration and verification. The comparison of the computed and observed inundation depths of high-water The results of the calibration and verification of the HEC-RAS model for the river flow during flood peak periods are shown in Figure 7. As shown in the results, the performance statistics showed acceptable agreement for both the model's calibration and verification. The comparison of the computed and observed inundation depths of high-water marks in the floodplain at the three locations, namely, Phuang-Payom, Phumin-Thali, and Mueang Len, were found to be satisfactory with percentage errors of 6.7, 6.3, and 9.1, respectively. This assures the accuracy of the model simulation in the floodplain.
The HEC-RAS model was applied to calculate the flood duration, flood depth, and flood velocity in the municipality's floodplain area for the flood return periods of 10, 50, 100, and 500 years. It was found that the computed flood duration, T, varied from 74 to 151 h; the flood inundation depth, D, from 1.5 to 2.97 m; the flood flow velocity, V, from 0.34 to 0.81 ms −1 .
From the results of the HEC-RAS model, flood hazard maps of the municipal area for the return periods of 10, 50, 100, and 500 years were constructed as shown in Figures 8 and 9 respectively.

Flood Hazard
In Equation (1) Table 2. They were substituted into Equation (1) to compute the FHI.
The hazard weighting factors α for flood duration, β for flood depth, and µ for flood velocity in Equation (1) were determined by AHP [24,25] as described in the Appendix A. It was found that weight α was 0.63 for the flood duration, β was 0.26 for the flood depth, and µ was 0.11 for the flood velocity. The sum of the three weights was equal to unity according to Equation (2).
The values of the FHI of each grid was computed for floods in the return periods of 10, 50, 100, and 500 years. In the municipal area of 7.6 km 2 , the total flooded area for the 100 year return period was found to be 1,614,545 km 2 , or 21.24% of the study area. The remaining part of the study area had no flood inundation. According to the ranges of classification, the hazard index of the flooded area of 1,614,545 km 2 was classified into a very high hazard area of 307,638 km 2 , or 19.05% of the study area; a high hazard of 25,622 km 2 , or 1.59%; a medium hazard of 68,432 km 2 , or 4.32%; a low hazard of 1,014,760 km 2 , or 62.85%; a very low hazard of 198,094 km 2 , or 12.26%. The high and very hazardous areas were in Phumin-Thali village in the south of the study area, and another major flooding location was in the area of Suriyapong Camp in the north of the study area.

Total Flood Vulnerability
The total vulnerability index, FVI, was determined for each village using Equations (3) and (4). The computed population vulnerability index, FVI pop , and the computed household vulnerability index, FVI hh , were substituted into Equation (3). The weights w pop for population vulnerability and w hh for household vulnerability were determined using AHP based on the data collected from the Nan Municipality and from the questionnaire survey. The weights w pop and w hh were found to be 0.33 and 0.67, respectively. This shows that the household vulnerability, FVI hh , has much more influence on the total vulnerability, FVI, than the population vulnerability, FVI pop .
The computed total vulnerability index, VI, was further classified into five ranges of equal intervals corresponding to the index values FVI of 1, 2, 3, 4, and 5 and named very low, low, medium, high, and very high, respectively, as shown in Table 2. Figure 10 shows the vulnerability level of different villages in the municipal area.

Flood Risk
In calculating the flood risk, the municipal area was divided into grids of 30 by 30 m. In each village in the municipality, the hazard varied from one point to another depending on the flood conditions and topography. For the total vulnerability, it was considered to be the same everywhere, as the village was the smallest unit in the questionnaire survey. The computed FRI for all grids in the same village was classified into five levels with equal intervals, namely: very low for 1 < FRI ≤ 5, low for 5 < FRI ≤ 10, medium for 10 < FRI ≤ 15, high for 15 < FRI ≤ 20, and very high for 20 < FRI ≤ 25. The maps of the computed flood risks for the 100 and 500 year return periods are shown in Figures 11 and 12 respectively. As shown in Figure 10, it was found that the areas along the right bank (west bank) of the Nan River and in the southern part near Phumin-Thali village were highly vulnerable and highly sensitive to floods compared to the other parts of the municipality.
From the computed total vulnerability index, FVI, of each village, it is found that 10 villages, namely, Chiang Kang, Phumin-Thali, Pouang Phrayom, Phra Nate, Mueang Len, Tha Chang, Don Khwan, Huea Khaung, Phrayaphu, and Don Sri Saim, located near the Nan River and its tributary Chao Fa canal have very high flood vulnerability.

Flood Risk
In calculating the flood risk, the municipal area was divided into grids of 30 by 30 m. In each village in the municipality, the hazard varied from one point to another depending on the flood conditions and topography. For the total vulnerability, it was considered to be the same everywhere, as the village was the smallest unit in the questionnaire survey. The computed FRI for all grids in the same village was classified into five levels with equal intervals, namely: very low for 1 < FRI ≤ 5, low for 5 < FRI ≤ 10, medium for 10 < FRI ≤ 15, high for 15 < FRI ≤ 20, and very high for 20 < FRI ≤ 25. The maps of the computed flood risks for the 100 and 500 year return periods are shown in Figures 11 and 12 respectively.   From the maps, it was found that the risk level and inundation area increased significantly by four times for the return period from 100 to 500 years. In Figures 11 and 12 , the highest risk areas consistently appear in the most southern part of the municipal area, near to the wastewater treatment plant, where the Phumin-Thali and Phuang-Payom villages are situated. There, the ground elevation is low; flood depth and flood duration are large; the population density and number of household properties are high. For the 100 From the maps, it was found that the risk level and inundation area increased significantly by four times for the return period from 100 to 500 years. In Figures 11 and 12, the highest risk areas consistently appear in the most southern part of the municipal area, near to the wastewater treatment plant, where the Phumin-Thali and Phuang-Payom villages are situated. There, the ground elevation is low; flood depth and flood duration are large; the population density and number of household properties are high. For the 100 year return period, more than half of the two villages are projected to be inundated and heavily damaged.
These two villages have the highest total vulnerability and highest hazard. Other locations with the highest or at high risk are along the right bank of the Nan River, especially in the northeastern part of the municipal area around Mueang Len village and Suriyapong Camp. The medium-and low-risk areas are in the inner part of the municipality approximately 1 or 2 km away from the riverbank. In the central part of the municipal area and in the northwestern part, the flood risk is very low or none, due to the fact of their higher ground elevations.

Discussions
Globally, the costs of flood losses and corresponding adaptations are generally preliminary and subject to a number of assumptions. Uncertainties are largely due to the fact of projected climate changes and estimation of flood damages as well as methodological shortcomings [42]. Options for adaptation measures and assumptions of vulnerability and exposure under various scenarios of socioeconomic conditions and physical environments are important factors in changes of flood risks. An attempt to overcome these shortcomings in risk assessments was conducted in the present study by developing a comprehensive and systematic method for flood risk assessment at the household level and to demonstrate its application in a municipal area of Nan Province in Thailand.
The discussion of the results of this study is as follows: 1.
The performance of the HEC-HMS model was evaluated using the following statistical parameters, namely: R 2 , NSE, PBIAS, Vr, and NRMSE. The results of the model's calibration at stations N64 and N1 were found to be satisfactory as shown in the results, and. the performance statistics showed acceptable agreement in both the model's calibration and verification; 2.
The comparison of the computed and observed inundation depths of high-water marks in the floodplains at Phuang-Payom, Phumin-Thali, and Mueang Len villages were found to be satisfactory with percentage errors of 6.7, 6.3, and 9.1, respectively. This assures the accuracy of the model's simulation in the floodplain. Importantly, the results show the reliability of the HEC-RAS flood simulation model and the data used in the calculation; 3.
The flood hazard was calculated by using Equations (1) and (2) in which the weights α, β, and µ representing flood duration, depth, and velocity, respectively, were systematically determined by AHP. In the other previous studies [17,18,20], the weights α, β, and µ for flood duration, depth, and velocity, respectively, in Equation (1) were specified by the researchers according to their experiences without using AHP. These weights could be subjective or questionable. By using the AHP, the weights α, β, and µ of 0.63, 0.26, and 0.11, respectively, were obtained. In AHP, the relative importance of a flood duration, depth, and velocity of 1, 3, and 5 was assigned according to the results of field and questionnaire surveys. A sensitivity analysis was conducted to determine the effects of change on the relative importance of flood duration, depth and, velocity from the original values of 1, 3, and 5 to the new values of 1, 5, and 9, respectively. By using AHP, the new weights α, β, and µ were found to be 0.72, 0.21 and 0.06, respectively. In percent, the changes were +14.3% for weight α, −19.2% for β, and −45.4% for µ. The changes in weights α and β were considered to be small and acceptable. The change in weight µ of −45.4% was negative and moderate. Based on the field and questionnaire surveys in this study, it was revealed that the damaging effect of flow velocity in the municipal area was much smaller compared to flood duration and depth; therefore, the change in weight µ of −45.4% was considered non-significant. Hence, weights α, β, and µ of 0.63, 0.26, and 0.11, respectively, were considered reasonable; As shown in Figures 8 and 9, the flood hazard increased both spatially and in magnitude with the increase in the flood return periods. The 10 year flood hazard along the right bank of the Nan River was smaller than the 50 and 100 year floods, and it was much smaller compared to the 500 year flood. Significant flood hazard occurs in the Phumin-Thali and Phuang-Payom villages in the southern part of the municipality. The locations of the observed high-water marks are shown by the small red circles in Figure 10. It can be seen that the center of the municipal area has higher ground elevation than the surrounding area and, hence, it has less of a flood hazard. On the other hand, in the Mueang Len area in the northeastern part of the municipality, the hazard is significant when the flood magnitude is greater than the 100 year return period.

4.
For total flood vulnerability, the weights W pop and W hh in Equation (3) were found to be 0.33 and 0.67, respectively. This shows that FVI hh had much more influence on FVI than FVI pop . For household vulnerability, FVI hh , the relative importance of the major contributing factors, namely, sensitivity F1, adaptation F 2 , and exposure F 3 were set to be 1, 3, and 5, respectively. By using AHP, the weights of w 1 of sensitivity F1, w 2 of adaptation F2, and w 3 of exposure F3 were found to be 0.63, 0.26, and 0.11, respectively. The same sensitivity analysis was conducted for the case of flood hazard, and it was found that the values of 0.63, 0.26, and 0.11 for weights w 1 , w 2 , and w 3 , respectively, were the most reasonable ones. In previous studies [21,23], the weights w 1 , w 2 , and w 3 were not determined by AHP but were assumed to be equal to one. Such an assumption could be incorrect, as the values of the sensitivity F1, the adaptation F2, and the exposure F3 were calculated on different bases, and they were not normalized. In this study, each major contributing factor was considered to have various contributing components Ci as shown in Table 3. The weight θi of each component Ci was determined by AHP based on the collected samples from the questionnaire surveys as given in Table 3; The distribution of the total vulnerability index, FVI, in the municipal area is shown in Figure 10. Depending on the sensitivity, adaptive capacity, and exposure of the households and population, the very high and highly vulnerable areas were found along the right bank of the Nan River, from upstream to downstream. These areas included Phumin-Thali, Phuang-Payom, and the Mueang Len villages. The medium vulnerable areas are in the central and western parts. Only two villages in the western rim of the municipal area, namely, Pha Mai and Don Swan, have very low vulnerability, because they have very low population densities and are located on higher ground elevations, far away from the river.

5.
Flood risks depend on flood hazard probabilities and vulnerabilities. Therefore, flood risks also change with flood probabilities or flood return periods. The study's results show that when the flood hazard changes, the flood risk also changes correspondingly; 6.
To mitigate flood problems in the municipal area, various flood control or mitigation measures can be proposed such as dredging of the Nan River channel and its tributaries, raising crest elevations of river flood control levees, or construction of flood bypass channel around the Nan municipal area. The effectiveness of each measure in reducing the flood hazard can be evaluated by using the hydrological model (HEC-HMS) and the hydrodynamic flood routing model (HEC-RAS). Hence, the changes in flood risk in Nan Municipality can be determined. More details can be found in [21].

Conclusions
The main results of the present study can be summarized as follows: 1.
The novelty of this study is the development of an advanced comprehensive and systematic methodology in determining flood hazard, total flood vulnerability, and flood risk at the household level in a municipal area. This is an important improvement over previous studies in which the flood and vulnerability parameters were not all considered. Moreover, these parameters were not systematically determined; 2.
The methodology was applied to a case study of a municipal area of 7.6 km 2 in Nan Province, Northern Thailand. The study area was located in the floodplain on the right bank of the Nan River; 3.
The HEC-HMS hydrological model and the HEC-RAS flood flow simulation model were applied to predict flood depths, velocities, and durations for the return periods of 10, 50, 100, and 500 years; 4.
The model computed results showed that significant flood hazards occur in Phumin-Thali and Phuang-Payom villages in the southern part of the municipal area and in Mueang Len village in the northeastern part of the area. These villages have low ground elevations, and they are near to Nan River. The central part of the municipal area had less of a flood hazard, as it has a higher ground elevation, and it is far from the river. The computed results were found to be consistent with the past flood situations during the field survey; 5.
The questionnaire survey in the municipal area reported that the flood duration had the most significant impact on households, while the flood depth and velocity had lesser impacts; 6.
In-depth analysis of the total vulnerability in the municipal area showed that the vulnerability of the population constituted one-third of the total vulnerability, while the household vulnerability constituted the remaining two-thirds; 7.
From the computed flood risks, flood risk maps were constructed for various return periods. The maps show that Phumin-Thali and Phuang-Payom villages, located near the right bank of the Nan River, are under very high risk, and more than half of the villages will be inundated and prone to high flood damages. This is consistent with past flood situations; 8.
The flood risk in the municipal area increased by approximately four times for the increase in the return period from 10 to 500 years; 9.
Overall, the methodology developed in this study yielded realistic results, and it should be applied to other study areas.
Author Contributions: The first author (T.T.) contributed 65%, and the second author (T.P.) contributed 35% of this research article. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
Most of the data and information used in this study such as river basin characteristics, river geometry, field survey and questionnaire survey can be obtained from Reference [21]. Other additional data can be obtained from References [22,[26][27][28][29][30].

Acknowledgments:
The authors express their acknowledgements to the Royal Irrigation Department, the Thai Meteorological Department, the Land Development Department, and the Nan Provincial Administration for their full cooperation in providing useful data in this study. Acknowledgments are also extended to the Asian Institute of Technology for the support in this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The AHP procedure starts from a pairwise comparison with a rating scale from 1 to 9. Generally, and in this study, scale 1 represents the most important flood parameter and scale 9 represents the least important. From the results of field and questionnaire surveys, the effect of flood duration was reported to have the highest impact on people's living conditions in the municipal study area and, hence, it was considered to be the most important. On the other hand, the flood depth and flood velocity were found to have less and lesser effects, respectively. A pairwise comparison of the three flood parameters was conducted, and the results were formulated in the form of a 3 × 3 judgment matrix as shown in Table A1. The diagonal elements of the matrix are equal to 1, as they are compared to themselves. In the first row, the matrix element values for flood duration, depth, and velocity were judged, according to their relative importance, to be 1, 3, and 5, respectively. According to the AHP, the successive elements of column 1 are the reciprocal of the row elements that is 1, 1/3, and 1/5 for flood duration, flood depth, and flood velocity, respectively. Next, from the second diagonal element of the 3 × 3 matrix, the element values of the row from the second diagonal element to the right-hand side were assigned to be 1 and 3 according to the relative importance of depth and velocity. For the column below the second diagonal element, the element values are the reciprocal of the mentioned row elements, that is 1 and 1/3, respectively.
The judgment matrix was then normalized. For example, in the first column of the matrix in Table A2, each element of column one was normalized by dividing the value of each element by the sum of all elements of the first column, i.e., 1 + 1/3 + 1/5 = 23/15. Hence, the normalized elements of column one in Table A2 are 15/23, 5/23, and 3/23. The same procedure was repeated for the remaining two columns. Then weight α for flood duration was calculated as the average of the normalized elements of the first row that is (15/23 + 9/13 + 5/9)/3 = 0.63. By the same procedure, weight β for flood depth and weight µ for flood velocity were calculated and equal to 0.26 and 0.11, respectively, as shown in Table A2. The sum of the three weights α + β + µ = 1.00 as given by Equation (2).
The judgment matrix was then normalized. For example, in the first column of the matrix in Table A2, each element of column one was normalized by dividing the value of each element by the sum of all elements of the first column, i.e., 1 + 1/3 + 1/5 = 23/15. Hence, the normalized elements of column one in Table A2 are 15/23, 5/23, and 3/23. The same procedure was repeated for the remaining two columns. Then weight α for flood duration was calculated as the average of the normalized elements of the first row, that is (15/23 + 9/13 + 5/9)/3 = 0.63. By the same procedure, weight β for flood depth and weight µ for flood velocity were calculated and equal to 0.26 and 0.11, respectively, as shown in Table A2. The sum of the three weights α + β + µ = 1.00 as given by Equation (2).
The judgment matrix was then normalized. For example, in the first column of the matrix in Table A2, each element of column one was normalized by dividing the value of each element by the sum of all elements of the first column, i.e., 1 + 1/3 + 1/5 = 23/15. Hence, the normalized elements of column one in Table A2 are 15/23, 5/23, and 3/23. The same procedure was repeated for the remaining two columns. Then weight α for flood duration was calculated as the average of the normalized elements of the first row, that is (15/23 + 9/13 + 5/9]/3 = 0.63. By the same procedure, weight β for flood depth and weight µ for flood velocity were calculated and equal to 0.26 and 0.11, respectively, as shown in Table A2. The sum of the three weights α + β + µ = 1.00 as given by Equation (2). The consistency of the calculated weights α, β, and µ was determined using consistency indicators such as principal eigenvalue (λ max ), consistency ratios (CR), and the consistency index (CI).
According to [24,25], the principal eigenvalue should be equal to the number of flood hazard parameters, which was three for satisfactory weights. The principal eigenvalue λ max , the CR, and the CI were computed using the following equations: where i is the row number of the judgment matrix in Table A1, j is the column number, λ max is the principal eigenvalue, n is the number of rows or columns of the square matrix, a ij is the element of judgment matrix, W i is the normalized weight of row i in Table A2, CR is the consistency ratio, CI is the consistency index, and RI is the random inconsistency index according to [24,25]. The principal eigenvalue (λ max ) computed using Equation (A1) was equal to (23/15 × 0.63 + 13/3 × 0.26 + 9 × 0.11) = 3.083. This was approximately equal to three, which is the number of flood hazard parameters. The consistency index and consistency ratio should be approximately zero for perfect weights. The CI was computed by Equation (A2) as: where n is the number of elements of the eigenvector, which was three. The CR = CI/RI = 0.0415/RI, where RI = 0.59 for n = 3 according to [24,25]. Therefore, CR = 0.0415/0.59 = 0.07.
It was found that both CI and CR were small numbers and were near to zero and, hence, they were very acceptable weight estimations.