An Efficient Modeling Approach for Probabilistic Assessments of Present-Day and Future Fluvial Flooding

Risk-informed flood risk management requires a comprehensive and quantitative risk assessment, which often demands multiple (thousands of) river and flood model simulations. Performing such a large number of model simulations is a challenge, especially for large, complex river systems (e.g., Mekong) due to the associated computational and resource demands. This article presents an efficient probabilistic modeling approach that combines a simplified 1D hydrodynamic model for the entire Mekong Delta with a detailed 1D/2D coupled model and demonstrates its application at Can Tho city in the Mekong Delta. Probabilistic flood-hazard maps, ranging from 0.5 to 100 year return period events, are obtained for the urban center of Can Tho city under different future scenarios taking into account the impact of climate change forcing (river flow, sea-level rise, storm surge) and land subsidence. Results obtained under present conditions show that more than 12% of the study area is inundated by the present-day 100 year return period of water level. Future projections show that, if the present rate of land subsidence continues, by 2050 (under both RCP 4.5 and RCP 8.5 climate scenarios), the 0.5 and 100 year return period flood extents will increase by around 15- and 8-fold, respectively, relative to the present-day flood extent. However, without land subsidence, the projected increases in the 0.5 and 100 year return period flood extents by 2050 (under RCP 4.5 and RCP 8.5) are limited to between a doubling to tripling of the present-day flood extent. Therefore, adaptation measures that can reduce the rate of land subsidence (e.g., limiting groundwater extraction), would substantially mitigate future flood hazards in the study area. A combination of restricted groundwater extraction and the construction of a new and more efficient urban drainage network would facilitate even further reductions in the flood hazard. The projected 15-fold increase in flood extent projected by 2050 for the twice per year (0.5 year return period) flood event implies that the “do nothing” management approach is not a feasible option for Can Tho.


INTRODUCTION
Flooding is one of the most frequently occurring and damaging natural disasters in the world (Hirabayashi et al., 2013;Kundzewicz et al., 2014;Arnell and Gosling, 2016;Alfieri et al., 2017;Forzieri et al., 2017;Mora et al., 2018). Coastal and delta areas are among the most flood hazard-prone areas of the world (Nicholls, 2004;Nicholls et al., 2007;Wong et al., 2014;Neumann et al., 2015;Pasquier et al., 2019), and flood intensity and frequency are already increasing, especially in coastal and delta cities, due to changes in river flows, downstream sealevel and local changes in rainfall and land use (Merz et al., 2010;Balica et al., 2012;Huong and Pathirana, 2013;Chen et al., 2018). Coastal and delta regions are also among the most densely populated areas of the world, which have experienced a rapid expansion of settlements, urbanization, infrastructure, economic activities, and tourism especially over the last 5 decades or so (Small and Nicholls, 2003;Valiela, 2006;Ranasinghe and Jongejan, 2018).
Flood risk assessment is performed to facilitate the implementation of risk-informed measures aimed at minimizing (or mitigating) flood damage (Hall et al., 2003;Meyer et al., 2009;Bureau Reclamation, 2020;Mishra and Sinha, 2020). Flood hazard estimation, which computes the probability and intensity of a possible event (Pappenberger et al., 2012;Alfieri et al., 2013;De Moel et al., 2015) is the first step in flood risk assessment (Penning-Rowsell et al., 2005;De Moel et al., 2015;Foudi et al., 2015;Kvočka et al., 2016). However, flood risk management and planning decisions in many parts of the world have historically utilized flood hazard or risk maps associated with one or two pre-determined return period water levels. While this may have been sufficient in the past, the need to move from stationary to innovative time-dependent (non-stationary) flood hazard and risk modeling approaches that can account for the uncertainty, arising from anthropogenic and climatic induced stressors, is rapidly increasing (Mosavi et al., 2018).
Climate change is now recognized as a major global challenge in the twenty-first century and beyond. Future projections indicate that climate change will have implications on the trends of sea level rise due to global warming and extreme events (e.g., extreme precipitation and hurricanes), which may lead to increases of flooding in the future (Panagouliaas and Dimoub, 1997;Menzel et al., 2002;Nicholls and Cazenave, 2010;IPCC, 2013Prudhomme et al., 2013;Alfieri et al., 2015). Besides the challenges posed by climate change, the population in coastal and delta zones is projected to increase in all Shared Socioeconomic Pathways (SSPs) by 2050, with especially the population living in the low elevated coastal zone projected to exceed one billion (by 2050) in all SSPs (Merkens et al., 2016) and likely to reach 1.4 billion by 2060 under the high end growth assumption (Neumann et al., 2015). Increases in population inevitably increase water demand, which is often satisfied by excessive groundwater extraction, which, more often than not, leads to land subsidence, further exacerbating the flood hazard due to increased inundation levels and frequency. Lowering of the inundated areas as well as the river banks due to land subsidence can contribute to this increase. The combination of increased flood hazard and increased population/infrastructure will, in turn, lead to an increase in flood risk in these vulnerable areas (Nicholls et al., 2007;Lenderink and Van Meijgaard, 2008;Syvitski et al., 2009;Min et al., 2011;Balica et al., 2012;Rojas et al., 2013;Wong et al., 2014;Takagi et al., 2015). However, projections of these trends contain uncertainties. This particularly holds for deltas due to their multi-faceted and dynamic character (Nicholls et al., 2020).
The present study attempts to address the challenge of estimating non-stationary fluvial flood hazard in a way it can inform risk modeling via a computationally efficient modeling approach based on a simplified 1D hydrodynamic model for the entire Mekong Delta (area of 40,577 km 2 ) that is coupled with a detailed 1D/2D coupled model and demonstrates its application at the urban center of Can Tho city (Ninh Kieu district) in the Mekong Delta.

ESSENTIAL ELEMENTS FOR PROBABILISTIC ASSESSMENT OF FLOOD HAZARDS
Robust quantification of flood risk assessment typically requires multiple (thousands of) river and flood model simulations to derive probabilistic flood hazard assessments. This constitutes a major challenge, especially for large river systems, such as the Mekong, due to the associated computational and resource demands. Although the last decade has witnessed a great improvement in computational capabilities and models, traditional modeling approaches still pose significant challenges in terms of computational time required to obtain fully probabilistic flood hazard estimates (McMillan and Brasington, 2007;Neal et al., 2012). Therefore, developing computationally efficient modeling approaches to circumvent this particular bottleneck remains an important challenge.
Many studies have attempted to improve the computational performance of 1D and 2D hydraulic models in different ways. These attempts include methods such as: (i) Model simplification-incorporating simpler representations of the physical processes to reduce the complexity of the model structure, e.g., ignoring certain terms such as inertia to simplify 2D shallow water equation (Bates and De Roo, 2000;Seyoum Solomon et al., 2012); applying Cellular Automata approaches in 1D drainage networks (Austin et al., 2014), and 2D overland flows (Ghimire et al., 2013); using conceptual models (Wolfs and Willems, 2013;Teng et al., 2015). (ii) Detail reduction-reducing the level of detail and input details of the model, and/or simulation time-step, e.g., using simplified drainage networks (Davidsen et al., 2017), simplified river networks (Ngo et al., 2018), or lower resolution topographic data (Savage et al., 2016). (iii) Using superior computational resources: parallel computation in 1D (Burger et al., 2014) and 2D models (Leandro et al., 2014;Zhang et al., 2014), and Cloud computing (Glenis et al., 2013).
These approaches significantly reduce model simulation times (Davidsen et al., 2017;Ngo et al., 2018), especially simplified models. Simplifying the model usually reduces the expressive ability (due to reduction of the behavioral complexity), statistical precision and sometimes the statistical accuracy of the model. The appropriateness of a simplified model should be looked at in terms of the purpose of the intended model application.
According to the "fit for purpose model" concept, the choice of a model for flood simulation should be done such that the model that can provide predictions with an acceptable level of statistical accuracy at a reasonable computational time and cost (Wright and Esward, 2013;Haasnoot et al., 2014). Statistical methods including cross-validation (often used with machine learning models), computation of goodness-of-fit with a normalized statistic like NSE (Nash-Sutcliffe efficiency, typically used with hydrologic models Nash and Sutcliffe, 1970), are used to establish the appropriateness of a simplified model to represent results at a point of interest. In addition to the computational performance of the models, the number of simulations is also a factor influencing the computational cost, especially time consuming 2D simulations. Flood risk assessment in the modern context often calculates flood water levels for different return periods (e.g., 5, 10, 20, 50, and 100 year). This reduces the number of 2D simulations that need to be executed while still achieving probabilistic flood hazard distributions.
Availability of historical data is a factor that influences flood hazard and risk assessment, in terms of flood frequency analysis (Machado et al., 2015;Engeland et al., 2018). Flood frequency analysis is used to determine the peak of flood events corresponding to a specified return period (Bayliss and Reed, 2001) in order to design flood risk reduction measures. However, the length of observed data is not always sufficient for robust flood frequency analysis. However, synthetic data approaches, such as the synthetic streamflow generator developed by Giuliani et al. (2017), provide to the means to generate a large numbers of time series data based on limited historical data.
In summary, key features of an effective non-stationary fluvial flood hazard modeling approach may therefore include (a) model reduction-a substantially simplified 1D model of the river system which can complete a simulation in short time span [e.g., completing a 1 year simulation of river flow (with daily time step) in around a minute (Ngo et al., 2018)], (b) strategic use of limited input data (e.g., using a synthetic streamflow generator to generate a large number (e.g., 1,000) time series of river flow based on limited historical data , (c) reduction of the number of time-consuming 1D/2D coupled model runs required to achieve probabilistic flood hazard results via flood frequency analysis while also considering flood hydrograph patterns, (d) derivation of probabilistic flood hazard quantification.
Model reduction is achieved in this study by gradually reducing the complexity of the detailed 1D model of the river system, using trials with varying levels of complexity (Ngo et al., 2018). Here the main river artery is retained while systematically removing small and medium distributaries, based on their size and/or distance from the target area (in this case, Can Tho city). In doing this, nodes at essential locations in the river system (e.g., at branching points, places at which the direction of flow changes, locations at which there are large change in river crosssectional area) are retained, and nodes at less crucial positions are sequentially removed. After each stage of the model reduction, the model is run for a specific period (Ngo et al., 2018). The model simulated water levels are compared with the measured water level data at the area of interest to evaluate the accuracy of the model. If the model performance is still good enough, the process of reduction is continued until the simplest level of detail that still provides accurate predictions of water levels in the area of interest is obtained. The final model thus obtained is referred to hereon as the simplified model.
In this application, due to limited data availability, the synthetic streamflow generator developed by Giuliani et al. (2017) is used to generate 1,000 synthetic river flow time series. This synthetic streamflow generator uses the non-parametric method to re-sample flows from the recorded data, which combines the methods of Nowak et al. (2010) and Kirsch et al. (2012). First, Kirsch's method is used to generate flows on a monthly time step. Specifically, the monthly flows in the historical data are log-transformed and standardized, which are then randomly sampled to generate standard normal monthly synthetic flows. Subsequently, Nowak's method is used to disaggregate these monthly flows to daily flows. In detail, historical monthly flows are used to calculate the k nearest neighbors for each generated monthly flow. The k-nearest neighbors are then sorted from the closest to the furthest, and probabilistically selected for proportionally scaling flows in disaggregation. After a neighbor (i.e., one of the historical months) is selected, all historical daily flows from the selected neighbor are proportionally scaled to generate synthetic daily flows . This synthetic streamflow generator is available on Github (https://github.com/julianneq/ Kirsch-Nowak_Streamflow_Generator). Here, Matlab code was used to generate correlated synthetic historical data at multiple sites, while Python code was used to validate the synthetic historical data statistically (https://waterprogramming. wordpress.com/2017/08/29/open-source-streamflow-generatorpart-i-synthetic-generation/).

APPLICATION OF THE FLUVIAL FLOOD HAZARD MODELING APPROACH IN CAN THO, VIETNAM
Can Tho (population 1.6 million in 2016) is one of the five cities that are directly administered by the Central Government of Vietnam, and is located in the center of the Mekong Delta, next to the Hau River (Bassac River) ( Figure 1A). In the coming decades, Can Tho is expected to be a dynamic city not only in the Mekong Delta but also in the entire southern region of Vietnam (Huong and Pathirana, 2013;MDP, 2013). The city includes five urban districts (Ninh Kieu, Binh Thuy, Cai Rang, O Mon, and Thot Not) and four rural districts (Phong Dien, Thoi Lai, Vinh Thanh, Co Do) ( Figure 1B). Situated in a tropical monsoon climate region, Can Tho has two distinct seasons: rainy and dry season. The rainy season usually starts from May and lasts until the end of November, with rainfall during this season accounting for 90% of the annual rainfall (Can Tho City People's Committee, 2010), while the dry season lasts from December to April with low rainfall. Frequent flooding is a major problem in Can Tho (Huong and Pathirana, 2013;Ngo et al., 2018), at least 2-3 times per year, with inundation depths of up to 50 cm at some places in the city center. Although Can Tho city is some 80 km upstream from the ocean, it is impacted by tidal water level variations, thus making the city vulnerable to the effects of projected sea-level rise and storm surges. Additionally, Can Tho is already facing many challenges from rapid urbanization, population growth, and poor infrastructure (e.g., city flood drainage network)-all of which combine to increases the city's vulnerability to floods. Figure 2 shows the highest observed flood water levels in Can Tho since 2000 and official flood water level alarms in Can Tho.
This study focuses on a part of Ninh Kieu district (population 280,000 in 2019) which is the central area of Can Tho city ( Figure 1C). Ninh Kieu is the most developed district of Can Tho and home to many trade centers, urban areas, and residential areas.
To achieve the goal of producing probabilistic flood hazard estimates for the study area (for the present and 2050), the method adopted here broadly comprised four methodological steps: (a) Model reduction-a substantially simplified 1D model for the entire Mekong Delta developed and validated by Ngo et al. (2018), (b) strategic use of limited input data (seven years data of river flow) to generate 1,000 time series data of river flow at Kratie using the synthetic streamflow generator, (c) reduction of the number of 1D/2D coupled model runs required by flood frequency analysis while also considering the flood hydrograph patterns, and (d) derivation of probabilistic flood hazard quantification for the present and for 2050 under RCP 4.5 and 8.5, with and without land subsidence. The methodology adopted is summarized in Figure 3. Data sources used for this study are shown in Table 1 and the way in which these data were used here is explained in detail below.

Model Reduction
A simplified 1D model for the entire Mekong Delta, which is developed using the open-source dynamic rainfall-runoff simulation model SWMM (Storm Water Management Model, Rossman, 2015) is used to derive the river water levels at Can Tho. Here, first, a detailed SWMM model of the entire Mekong Delta was built based on an existing ISIS model, which comprised 575 nodes and 592 links ( Figure 4A). Using trial simulations with the varying levels of complexity, nodes and links were systematically removed to arrive at a much reduced SWMM model comprising 37 nodes and 40 links ( Figure 4B) which was able to accurately simulate 1-year of river water levels at Can Tho city in under 1 min. The model was calibrated with the year 2000 observed water level data and then validated with the years measured water level data at Chau Doc, Tan Chau and Can Tho gauging stations in 2001, 2002(Ngo et al., 2018. Years 2001 and 2002 are years that typical flood events occurred, while 2011 is a year that an extreme event occurred (on 27th October) (Figure 2). The simplified model was validated for the 2011 event with the goal of investigating the model's capability to simulate the extreme events, which may occur more often in the future due to Climate change (Panagouliaas and Dimoub, 1997;Menzel et al., 2002;Prudhomme et al., 2013;Alfieri et al., 2015). The general model performance ratings with respect to the October 2011 event based on two indicators NSE (Nash-Sutcliffe efficiency) and RMAE (Relative Mean Absolute Error) are 0.77 (very good) and 0.04 (excellent), respectively (Ngo et al., 2018).
The simplification of the model (e.g., removing small and medium distributaries, not considering the irrigation systems along the river) for the entire Mekong Delta leads to a degradation of the precision and accuracy of its results at locations far away from Can Tho, specifically at Chau Doc and Tan Chau (Ngo et al., 2018). Moreover, the simplified model is a 1D model, which cannot accurately simulate flood propagation, especially on floodplains, even though the floodplains have been included in cross-sections and assigned appropriated roughness coefficients. Therefore, the simplified 1D model is used here as a type of surrogate model for simulating water levels along the main rivers, bearing in mind that it may not provide reliable information on water levels and inundation dynamics in areas that are far away from Can Tho or located some distance from the main rivers on the floodplains.

Upstream Boundary Condition
In the application here, 7 years (2000)(2001)(2002)(2003)(2004)(2005)(2006) of observed flow data at the upper boundary location "Kratie" were used as forcing for the above described simplified 1D model of the Mekong River. As 7 years of data is not sufficient to derive probabilistic results, here a synthetic streamflow generator developed by Giuliani et al. (2017) was used to generate 1,000 synthetic flow time series (each 1 year long) for each scenario in Table 2 (current and future). While it is acknowledged that this approach does not fully replace the utility of long term observational data, due to the probabilistic nature of the streamflow generator, it nevertheless captures the statistical variation of upstream flows compared to simply using the 7 years of available data, which is important in flood hazard modeling.
For the future simulations that include the effects of climate change, here we used the annual river discharge projections of Hoang et al. (2016) under RCP 4.5 and 8.5 as the basis for generating the future river flow data. Hoang et al.'s (2016) projections indicate that annual river flow for RCP 4.5 and RCP 8.5 at Kratie in 2050 are expected to change between 3 to 8% and −7 to 11%, respectively, relative to the 1971-2000 baseline period adopted in that study. Future river flow data were generated by combining the flow of each year in 7 years of observed flows used here with a randomly selected % change of these projected changes in the flow corresponding to each RCP. Subsequently, the streamflow generator was used to generate 1,000 future riverflow time series (each 1 year long) corresponding to RCP 4.5 and RCP 8.5.

Downstream Boundary Condition
Thirty 6 years (1979-2014) of model simulated extreme sea level (tide + surge) data were extracted at the Mekong river mouths (Tran De, Ben Trai, and An Thuan) from the GTSR data set presented by Verlaan et al. (2016) to use as the downstream boundary condition of the 1D model.
For the future simulations, the above present-day extreme sea levels from the GTSR data were combined with the 2050 regional sea level rise projections presented by the Viet Nam Ministry of Nature Resource and Environment (MONRE, 2016) under RCP 4.5 and RCP 8.5 for the region containing the aforementioned river mouths. The projected regional sea level rise in the area by 2050 (relative to 1986-2005) for RCP 4.5 and RCP 8.5 are 13-32 cm and 16-35 cm, respectively. The simplified 1D model was then executed with 1 year time series of the above described boundary conditions to generate river water level at Can Tho. In all, 36,000 simulations were FIGURE 3 | The methodological framework adopted to derive probabilistic flood hazard estimates for the study area. This process was applied for the present and for 2050 (RCP 4.5 and RCP 8.5 together with 3 different local land subsidence scenarios).
undertaken corresponding with the total number of possible combinations of upstream (1,000 1-year time series of riverflow) and downstream (36 1-year time series of sea level) boundary conditions, resulting in 36,000 water level time series data (each 1-year long) at Can Tho per considered scenario (current and future).

Reduction of the Number of 1D/2D Model Runs
From the 36,000 1-year long water level time series at Can Tho, the maximum water level of each year in the 36,000 simulation years was extracted and used to fit an extreme value distribution (Gumbel or Type I). The water level at Can Tho corresponding to each return period for each scenario (current or future) was determined based on the fitted Gumbel distribution.
Apart from the maximum water levels, it is also important in flood hazard modeling to represent the shape of hydrograph around the peak water level. To achieve this, here we used a threshold-based method, as is commonly used in flood frequency analysis (Lang et al., 1999;Bezak et al., 2014). A threshold value of 2.15 m (i.e., maximum measured water level at Can Tho, see Figure 2) was adopted here and annual water level time series (of the full 36,000 series) that contain at least one water level value exceeding 2.15 m were identified.
For each peak flood value, 24 h long time series around each peak value (12 h earlier to 12 h later) were extracted. Due to there can be many different hydrograph shapes whose peak coincides with a design water level. Such time series each can have a unique shape, requiring, in principle, 1D/2D simulation run for each such time series. Since 1D/2D coupled model is the most expensive in terms of computational effort, this is a practically impossible task. This was circumvented by first, identifying a small number of representative flood hydrograph patterns and assuming that all flood hydrograph shapes with peak value H can adequately be represented by a hydrograph shape with peak value H can adequately be represented by a hydrograph shape F i,H (t) with unique hydrograph shape patterns i = {1,. . . , n}. For each pattern i, first the hydrograph shape having the maximum peak (called is typical flood hydrograph shape) within that pattern was selected as a representative pattern (f i (t)). F i,H (t) is calculated as follow: Where t max is the time at which the maximum value in representative pattern occurs.

Derivation of Probabilistic Flood Hazard Quantification
Several different flood parameters can be used to quantify the flood hazard, including inundation level, flow velocity, frequency of flooding, and flood duration, etc. (Ramsbottom et al., 2003;Ward et al., 2011;De Moel et al., 2015). Of these, inundation level (water depth) and flow velocity are considered the most important parameters (Penning-Rowsell et al., 1994;Wind et al., 1999;Merz et al., 2007;Kreibich et al., 2009). However, due to the relatively flat terrain combined with small inundation depths in Can Tho, the effect of the flow velocity is expected to be small compared to that of the flood inundation depth (Dinh et al., 2012). While in agricultural contexts as well as indirect damages (e.g., loss of livelihoods and nuisance), the duration plays an important role, in the context of urban property damage (e.g., buildings, furniture, and road), inundation depth is much more important than the duration. Hence, this study considers inundation levels as the main indicator of the flood hazard in the study area.
To simulate the effect of flood drainage network on flooding in the study area, here we used the detailed 1D urban model developed for the study area by Huong and Pathirana (2013). This model comprises 479 junctions, 612 conduits, 48 outfalls, and 303 sub-catchments (Figure 5).
The model parameters were calibrated for the flood event on the 17th of October 2016 based on the observed water depths (at 1-min measurement interval) in the manholes at Nguyen Van Cu and Tran Hung Dao streets (Figure 5). These are the only available observed flood water depths in the study area. Here, SWMM5-EA software (Pathirana, 2014), was used to calibrate the most uncertain parameters, i.e., Manning's roughness coefficient of conduits, Manning's roughness coefficient of the previous/impervious surfaces of the sub-catchments and the slope of the sub-catchments. The detailed 1D urban model performance was evaluated by using the NSE indicator, which determines the relative magnitude of the modeled variance compared to the observed data variance (Nash and Sutcliffe, 1970). NSE is computed as shown in Equation (1) in Appendix.
The calibrated 1D urban model was then used to establish an integrated 1D/2D model using PCSWMM software (http:// www.chiwater.com/Software/PCSWMM). This 1D/2D coupled model was validated based on inundation depths and flood extent as given in the report of Can Tho Water Supply and Drainage Construction Company for the same flood event on the 17th of October 2016.
Flood simulations were undertaken for water levels with return periods ranging from 0.5 to 100 year obtained from the flood frequency analysis described in section Strategic Use of Limited Input Data for each model scenario in Table 1. A 15 m resolution DEM developed by the Vietnam Institute of Meteorology, Hydrology, and Environment was used in all simulations. In the simulations that include land subsidence, an average subsidence rate of 1.6 cm/year for the entire Mekong Delta (Erban et al., 2014;Minderhoud et al., 2015) was used to consider the effect of land subsidence on the flood hazard at Ninh Kieu district. Being an urban area with generally less groundwater recharge and potentially more groundwater pumping combined with the contribution of other possible human-induced drivers (e.g., loading by infrastructure), it is very likely that Can Tho's subsidence rate is much higher than the average for the Mekong Delta. Although, in recent studies by Minderhoud et al., the average subsidence rate in the Mekong Delta was estimated to be 1.1 cm year −1 (Minderhoud et al., 2017) and 1.31 cm year −1 corresponding to B2 scenario (no-mitigation with a steady annual increase of 4% of the 2018 volume) (Minderhoud et al., 2020). However, this subsidence rate is likely to increase in the future due to increased groundwater demand (Minderhoud et al., 2017(Minderhoud et al., , 2020. Thus, this study used a subsidence rate of 1.6 cm year −1 in considering the effect of land subsidence on the flood hazard at Ninh Kieu district in the future. For all model scenarios indicated in Table 1, flood inundation maps were first developed for water level return periods ranging from 0.5 to 100 year (eight inundation maps for each model scenario). The flood maps thus obtained were added into ArcMap under shapefiles, which were then converted to raster files to extract the maximum inundation depths (during the flood events) corresponding to each return period by using the "Polygon to Raster" tool in Conversion tools of ArcToolbox. Subsequently, maximum inundation depths at all grid cells were aggregated to generate flood hazard maps for each return period.

Design Water Levels
Figures 6A-C show flood frequency curves at Can Tho for the present, and for 2050 corresponding to RCP 4.5 and RCP 8.5, obtained from the Gumbel distributions that were fitted to the modeled maximum water level at Can Tho based on the results of 1D hydrodynamic model described in Section Model Reduction. Table 3 shows the water level at Can Tho with return periods ranging from 0.5 to 100 year determined from the flood frequency analysis for each model scenario.

Model Calibration
The calibration results for the 1D urban drainage model for the flood event of October 2016 are shown in Figure 7. The

NSE indicator values between the measured and simulated (present conditions) water depths in the manholes at Nguyen
Van Cu and Tran Hung Dao streets is 0.75 (good) and 0.95 (very good), respectively. Thus, the 1D urban drainage model performance for the study area can be considered to be sufficiently accurate. The calibration of 2D results (inundation extent and depth) proved to be much more challenging due to the lack of observed data. An indirect visual approach was used to compare observed spot values of inundation heights with simulated depths. The validation results for the 1D/2D coupled model for the flood event of October 2016 are shown in Figure 8. Simulated flood extent and inundation depths at many different streets in Ninh Kieu district are presented in the description of the flooding situation and observed inundation depths in the report of Can Tho Water Supply and Drainage Construction Company regarding this event. The comparison shows a good agreement between the simulated flood extent and inundation depths and observations. Figure 9 shows the shape of the 24-h time series of all the modeled extreme water level (i.e., peaks greater than the threshold value of 2.15 m) time series at Can Tho (of all 36,000 annual time series) corresponding to RCP 4.5 and RCP 8.5. There are two dominant hydrograph shapes in the two RCPs considered (indicated by Pattern 1 and Pattern 2 in Figures 9A,B), but it is clearer in RCP 8.5.

Boundary Conditions for 1D/2D Model
Equations (1) was as used to calculate the 24-h long water level time series for a water level corresponding to 100-year return period ( Table 3) for Pattern 1 and Pattern 2 to create the 24-h boundary condition time series for the 1D/2D flood model to examine the response of different river water level hydrograph shapes on flooding.
The results (Figure 10) indicate that flood extent and inundation depths associated with Pattern 1 are greater than those associated with Pattern 2. The inundated area (for cells with maximum inundation depths ≥ 0.02 m) corresponding to 100-year return period water level for RCP 4.5 and RCP 8.5 with Pattern 1 are 1.96 and 2.30 km 2 , respectively. The corresponding flooded areas for Pattern 2 are lower at 1.41 and 1.81 km 2 . Therefore, in this study, the typical water level time series following Pattern 1 above was used to calculate flood hydrograph shape for each water level corresponding to each return period ( Table 3) for each scenario based on Equation (1) to create the 24-h boundary condition time series for the 1D/2D. This was a necessary simplification to reduce the complexity and computational burden of computing flood inundation. Figure 11 shows flood hazard maps corresponding to 5 and 100 year return period water levels (at Can Tho) for the present (model scenario #1), and 2050 under RCP 4.5 and RCP 8.5 (without land subsidence) (model scenarios #3 and #5, respectively). The remaining flood hazard maps which   correspond to 0.5, 1, 2, 10, 20, and 50 year return period water levels for these scenarios are shown in Appendix Figures A1-A3. Figure 12 and Appendix Figure A1 show that, under present conditions, 4.9% of the study area will be flooded even with the 0.5 year return period water level at Can Tho, especially some areas along the canals and lakes in the city that are connected to the Can Tho River and Cai Khe channel. Notably, there is a large area of Cai Khe ward (shown by the red circle of Appendix Figure  A1a), which is located close to the junction of Hau river, Can Tho river and Cai Khe channel (see Figure 1C) that is inundated under this condition.

Flood Hazard in Can Tho City
At the other higher end of the modeled return periods under present conditions, Figures 11B, 12 show that the inundated area with a 100 year return period water level at Can Tho, is more than double that with a 0.5 year return period water level at Can Tho increasing from 4.9% (0.5 year RP) to 12.7% (100 year RP).
For the future modeled scenarios (without land subsidence; model scenarios #3 and #5), the % inundated areas for the 100 year RP water level at Can Tho are 29.2 and 34.2%, under RCP 4.5 and RCP 8.5, respectively, representing more than a doubling of the inundated area by 2050, relative to the present.
Flood hazard maps corresponding to 5 and 100 year return period water levels (at Can Tho) for the future scenarios (2050under RCP 4.5 and RCP 8.5) that do take into account land subsidence (model scenarios #2 and #4) are shown in Figure 13. The remaining flood hazard maps which correspond to 0.5, 1, 2, 10, 20, and 50 year return period water levels for these scenarios are shown in Appendix Figures A4, A5. The difference is immediately visible with a large increase in the inundated area compared to when land subsidence is not taken into account. Figure 13 show severe flooding in a vast majority of the study area for all return periods of water level. The percentage area flooded by the 0.5 year RP water level under RCP 4.5 and RCP 8.5 are 73.4 and 75.3%, respectively (Figure 14), almost a 10fold increase relative the comparable projections without land subsidence and a 15-fold increase relative to present-day flooding due to the same RP water level. For the 100 year RP water level, the percentage area flooded under RCP 4.5 and RCP 8.5 projected to be 95.2 and 96.4%, respectively, representing a 3-fold increase relative to the comparable projections without land subsidence and a 8-fold increase relative to present-day flooding by the same RP water level. It is also noteworthy that the maximum inundation depths for the 100-year return period water level under RCP 4.5 and 8.5 are not much different at 3.17 and 3.20 m, respectively.
The above results highlight that while climate change will increase the flood hazard in the study area, land subsidence has a much greater effect than climate change driven variations in river flow on the flood hazard in the study area. Furthermore, clearly, the existing urban drainage network is not able to effectively drain flood waters even for present-day conditions, and this FIGURE 12 | The inundated area and percentage of the flooded area relative to total study area (also indicated by numerics along the solid lines in this figure) corresponding to each return period of water level (without land subsidence).
will be felt more severely in the coming decades. A significant reduction in groundwater extraction, which is the main cause of land subsidence in Can Tho (Erban et al., 2014) combined with a new and substantially efficient urban drainage network may be able to mitigate the projected flood hazard in the study area. It is recommended that the efficacy of these mitigation measures be investigated in detail in future modeling studies.

DISCUSSION
In this section, we discuss three aspects that are important in interpreting the findings of this paper: First, we explain the reason for using a short period of flow data in conjunction of a synthetic streamflow generator as opposed to a longer duration of observed data. Next, the rationale and limitations of using a probability distribution function (Gumbel distribution in this case) instead of empirical probabilities are discussed. There were several studies that discussed flood hazard/risk of Can Tho city. Of particular relevance is the study by Apel et al. (2016). So, finally we compare the current study with this past study.

Reliability of Using the Synthetic Streamflow Generator and Short Discharge Data
The most obvious reason to use synthetic hydrology stems from the fact that there is little or no data for the system (Lamontagne, 2015). This study used the synthetic streamflow generator to generate 1,000 annual synthetic river flow time series based on 7 years of discharge data at the upper boundary point, Kratie.
The river flow data is publicly available (e.g., the website of the Mekong River Commission) covering periods 1924-1970 and 2000-2018. Therefore, the rationale for using short data series to generate synthetic river flow data series in this study need to be justified.
The purpose of this study was to quantify climate change driven variations in the flood hazard between the present period and future time periods. Therefore, it is important to ensure that the selected baseline period (and baseline simulations) are in fact representative of the present-day period. This is important because, not only has the climate change signal emerged in several climate variables over the last 50 years or so (i.e., signal is clearly discernible from the inter-annual variability) (King et al., 2015), but also human activities (e.g., reservoirs) have led to noticeable changes in the natural regimes that may have existed earlier in the twentieth century (see Ranasinghe et al., 2019 for example in China). Both of these phenomena may change the probability distribution of climate variables over time (Chadwick et al., 2019).
The discharge time series at Kratie was analyzed, based on the 66 years of data (1924-1970 and 2000-2018). The analysis showed that the peak discharge at Kratie has indeed noticeably decreased over time, and particularly after 2000 (Figures 15,  16), likely due to irrigation expansion and upstream dam construction in recent years (MRC, 2010;Piman et al., 2013).
The use of the full discharge time series at Kratie to develop flood frequency curves is therefore inappropriate in the present study which aims to quantify climate change driven variations in the flood hazard, and further, risk relative to present-day conditions, in order to inform the development of climate-resilient flood risk reduction measures for the urban center of Can Tho city. The use of the full observed discharge data at Kratie, including pre-2000 flow with large flood peaks, can lead to an overestimation of flood hazard and risk. Therefore, for the purposes of this study, only the post-2000 discharge data were used to represent baseline conditions. Figure 17 shows several representations of synthetic flow time series that are generated based on the 7 years of observed discharge data used here, and their corresponding statistics and extreme values.
In all, 1,000 synthetic flow time series were created, which were then combined with 36 sea level time series to have 36,000 different water level time series at Can Tho. Using these 36,000 water level time series, does not add any information that was originally not present in the observed data. However, as the sealevel and river flow time series are independent of each other, these combinations of statistical realizations of streamflow with observed sea-level are assumed to improve the joint-probability manifestation in the resulting longer time series. It should be noted that the synthetic generator is not the only approach to achieve this. For example, a time-shifting one set of series against the other may give a similar outcome.

Limitation of Probabilistic Distribution Function
In this study, Gumbel distribution was used to model this distribution of the maximum water levels of the water level time series to select design water levels corresponding to each return period. However, there is a difference in probabilities between the empirical and distribution quantiles corresponding to the present and RCP 4.5 scenarios (Figures 6A,B). Using the log-log-linearity of the Gumbel distribution might introduce a bias to the very extreme values. In contrast, using the empirical distribution might result in less information loss. However, it will include all the random artifacts in the observed (and generated) data. Additionally, such an analysis needs a long series of data that are largely devoid of non-stationarities (or them being carefully removed). The flow database used here, which was limited in length due to the climate change impact focus of the study does not provide the necessary data quality or quantity to do such an analysis. Therefore, using the Gumbel distribution here is more reasonable in our view. Additionally, the difference in water levels corresponding to large return periods between the empirical and distribution quantiles is small. Furthermore, the Gumbel distribution was well-fitted in scenario RCP 8.5 with higher water levels ( Figure 6C).

Comparison With Previous Studies
Computation of present-day flood hazard and probabilistic flood maps using 2D models for Ninh Kieu district has also been done before by Apel et al. (2016). However, the approach adopted in the present study differs from that adopted in previous studies and has added value by improving the computation of flood hazard of Ninh Kieu district. Furthermore, this study takes a step forward from previous studies, being the first study to probabilistically compute future flood hazard in the study area under climate change. The main value additions of this study, compared to Apel et al.'s (2016) study, are discussed below.

Difference in Using Flood Probabilities to Develop Probabilistic Fluvial Flood Hazard Maps
One of the biggest differences between the present study and that reported by Apel et al. (2016) is the length of the river discharge time series used for flood frequency analysis. This difference is, in part, due to the different aims of the two studies: Apel et al.'s (2016) aim was to develop flood hazard maps for the present-day while the focus of the present study is to quantify climate change driven variations in the flood hazard.
Consistent with the aim of their study, Apel et al. (2016) used flood frequency curves at Kratie of Dung et al. (2015), which were constructed based on the longest possible time series of river discharge at Kratie, spanning 88 years . In contrast, as mentioned in Section Reliability of Using the Synthetic Streamflow Generator and Short Discharge Data, the purpose of the current study was to quantify climate change driven variations in the flood hazard between the present period and future time periods. Therefore, it is important to ensure that the selected baseline period (and baseline simulations), is in fact representative of the present-day period.
While Apel et al.'s (2016) introduced flood probabilities at Kratie into a combined large-scale inundation model (arguably a more realistic representation than the current 1D model) for the entire Mekong Delta (Dung et al., 2011) with a (more detailed) 2D representation of the Ninh Kieu district. Floods strongly vary over space (Nied et al., 2017;Vorogushyn et al., 2018). This spatial variability of flooding would influence the flood levels at Can Tho which is about 430 km downstream of Kratie. Moreover, the river water level at Can Tho and the resulting flood extent and inundation depth in the Ninh Kieu district are affected by the downstream sea level, especially high tides and storm surges (Huong and Pathirana, 2013). Thus, using flood probabilities at Kratie to develop probabilistic fluvial flood hazard maps for the Ninh Kieu district without considering the effects of spatial variation of flood probabilities and downstream sea level could lead to some uncertainties in the flood hazard computed at Can Tho. The present study overcomes these shortcomings by undertaking 2D flood modeling for Ninh Kieu district based on flood frequency analysis at Can Tho (as opposed to Kratie) and by taking into account both river discharges and downstream sea level in computing the river water levels near Can Tho.

The Difference in Flood Extent
Comparison of the results between the two studies shows substantial differences in the flood extent corresponding to different return periods for present-day. The inundated area corresponding to 2, 5, 10, 20, 50, and 100 year return period in Apel et al.'s (2016) study was 2.37, 3.33, 3.71, 4.30, 4.98, and 5.29 km 2 , respectively. In contrast, the inundated area for the presentday in this study was 0.42, 0.49, 0.54, 0.60, 0.74, and 0.85 km 2 , respectively. Apart from the two key methodological differences between the two studies highlighted above, there are also two other reasons that may have led to these differences in estimated present-day flood extents.
While the present study explicitly accounted for the effect of the urban drainage system in Ninh Kieu district on flooding, Apel et al. (2016) considered the entire district to be impervious. This has significant implications in terms of flood hazard estimations. The river water level in Can Tho varies following the downstream tidal fluctuation (semi-diurnal tide), as the urban center of Can Tho (Ninh Kieu district) is connected with the Hau River and Can Tho River via the open sewer channel and urban drainage system. Therefore, for, e.g., if during the flood phase of tide, the river water level rises above the lowest elevation of the top of the manholes in the city, although without necessarily being higher than the crest elevation of the river embankment, this will lead to flooding in the city center due to backwater flow through the urban drainage/sewer systems (Figure 18) (note: while the drainage outlets have flap-gates provided, they are largely dysfunctional due to solid waste blockages and lack of maintenance in Ninh Kieu district). This is consistent with the flood situation in Can Tho, which was described by Nguyen (2016) (http://www.cantholib.org.vn:84/Ebook.aspx?p= 27B9F975353796A6E64627B93B65654746C6B65637B91B857557). When the river water level drops during the ebbing phase of the tide, the inundation level is also reduced mostly as flood water is drained through the urban drainage/sewer systems. Hence, incorporating the effects of the flood drainage system, as done in the present study is crucial for correctly estimating flooding in this study area.
Both studies used the DEM presented by Huong and Pathirana (2013) for the study area as the input data of the 2D model. However, stemming from the above mentioned lack of consideration of the effects of the urban drainage/sewage systems, Apel et al. (2016) adjusted the elevation of the DEM data by subtracting 0.5 m from the original DEM in order to achieve an acceptable validation of their 2D model. Apel at al. (2016) FIGURE 16 | Gumbel distribution of discharge peaks at Kratie corresponding to three periods (1924-1950), (1951-1970), and (2000-2018). justify this decision referring to the two large fluvial flood events that occurred in 2011, with "extraordinary" peak water levels, but "the banks as given in the DEM were not overtopped, and thus no inundation would occur". However, revisiting the data of water levels at Can Tho station in 2011, used in Ngo et al. (2018) to validate the 1D simplified model for the entire Mekong Delta, the peak water levels of these two events occurred on the 28th of September and 27th of October with peak water levels of 2.04 and 2.15 m, respectively. Both these water levels are higher than the bank elevation extracted from the original DEM data (∼1.9-2.0 m) at the surveyed point in Apel et al. (2016). Thus, these two flood events would, in reality, have caused flooding in the Ninh Kieu district by both backflow through the urban drainage system and by direct overtopping of the river embankment. Urban floods like those of Can Tho are complex and highly non-linear processes. Therefore, compensating with a different variable (DEM datum) to address the issue of neglecting another (numerous connections between river and the city below the embankment level) does not result in reliable model performance. The lowering of the entire DEM is therefore the likely cause for the substantially larger present-day flood extents estimated by Apel et al. (2016), relative to those computed in the present study.

CONCLUSIONS
An efficient modeling approach that combines a simplified 1D hydrodynamic model with a detailed 1D/2D coupled model was developed and demonstrated at Can Tho city in the Mekong Delta. Key features of the modeling approach include (a) Model reduction-a substantially simplified 1D model for the entire Mekong Delta (area of 40,577 km 2 ) which can simulate 1 year of riverflow (with an hourly time step) in under 60 s, (b) strategic use of limited input data (7 years data of river flow) to generate 1,000 time series data of river flow at Kratie using the synthetic streamflow generator, (c) reduction of the number of 1D/2D coupled model runs required by performing flood frequency analysis while also considering flood hydrograph patterns, and (d) derivation of probabilistic flood hazard quantification for the present and for 2050 under RCP 4.5 and 8.5, with and without land subsidence. The detailed 1D/2D coupled model was successfully calibrated and validated against measured flood depths at two manholes, inundation depths and flood extent at many streets in Ninh Kieu district during the October 2016 flood event. Important improvements of this study over those of the past covering the same area are: (1) Consideration of the drainage network that provides hydraulic connections between the city and the river below the embankments; (2) Establishing the inundation hazard and risk as a function of the river levels near Can Tho, which is influenced both by upstream flow and downstream sealevel impact; (3) Develop a relatively numerically inexpensive approach to connect probabilistic urban flood hazard to up and downstream boundary conditions of the larger Mekong Delta area.
Flood hazard maps showing the maximum inundation depth during a flood event were developed for water level return periods ranging from 0.5 to 100 year. Analysis of the flood hazard maps indicates that even under present conditions, more than 12% of the study area will be inundated by the 100 year return period water level. With climate change, but without land subsidence, the 100 year return period flood extent is projected to more than double by 2050, with not much of difference between the two climate scenarios considered (RCP 4.5 and RCP 8.5). However, if the present rate of land subsidence will continue in the future, by 2050 and under both RCP 4.5 and RCP 8.5, the 0.5 and 100 year return period flood extents are projected to increase by around 15-and 8-fold, respectively, relative to the present-day flood extent.
These results indicate that reducing the rate of land subsidence, for example, by limiting ground water extraction, would substantially mitigate future flood hazards in the study area. Combining such a measure with a new and more efficient urban drainage network would further reduce the flood hazard. Future modeling studies are needed to quantitatively assess the hazard and risk reduction afforded by these adaptation measures, which could directly feed into risk informed adaptation measures and pathways. For Can Tho, the "donothing" management option does not appear to be an option given especially the 15-fold increase in flood extent projected by 2050 for even the twice per year (0.5 year return period) flood event.

DATA AVAILABILITY STATEMENT
The related data can be made available upon publication.
Requests to access these datasets should be directed to AP a.pathirana@un-ihe.org.

AUTHOR CONTRIBUTIONS
HN, AP, and RR prepared the initial version and the edited versions of the manuscript. All authors reviewed the manuscript and provided editorial input. All authors contributed to the article and approved the submitted version.