Reproducing an extreme flood with uncertain post-event information

Studies for the prevention and mitigation of floods require information on discharge and extent of inundation, commonly unavailable or uncertain, especially during extreme events. This study was initiated by the devastating flood in Tegucigalpa, the capital of Honduras, when Hurricane Mitch struck the city. In this study we hypothesized that it is possible to estimate, in a trustworthy way considering large data uncertainties, this extreme 1998 flood discharge and the extent of the inundations that followed from a combination of models and post-event measured data. Postevent data collected in 2000 and 2001 were used to estimate discharge peaks, times of peak, and high-water marks. These data were used in combination with rain data from two gauges to drive and constrain a combination of well-known modelling tools: TOPMODEL, Muskingum–Cunge–Todini routing, and the LISFLOOD-FP hydraulic model. Simulations were performed within the generalized likelihood uncertainty estimation (GLUE) uncertainty-analysis framework. The model combination predicted peak discharge, times of peaks, and more than 90 % of the observed highwater marks within the uncertainty bounds of the evaluation data. This allowed an inundation likelihood map to be produced. Observed high-water marks could not be reproduced at a few locations on the floodplain. Identifications of these locations are useful to improve model set-up, model structure, or post-event data-estimation methods. Rainfall data were of central importance in simulating the times of peak and results would be improved by a better spatial assessment of rainfall, e.g. from radar data or a denser rain-gauge network. Our study demonstrated that it was possible, considering the uncertainty in the post-event data, to reasonably reproduce the extreme Mitch flood in Tegucigalpa in spite of no hydrometric gauging during the event. The method proposed here can be part of a Bayesian framework in which more events can be added into the analysis as they become available.


Introduction
Losses caused by natural hazards have a significant impact on the world economy, and floods account for around half of all disasters globally (UN/ISDR, 2016).Prevention and mitigation of floods require information on discharge and extent of inundation.Such information is commonly unavailable or uncertain, especially during extreme events when gauging equipment becomes insufficient or is lacking.Data scarcity is further aggravated in developing countries with weak infrastructure.
Nearly 11 000 people were killed in Central America during Hurricane Mitch because of extreme flooding, an estimated 2.7 million lost their homes, and flood damages were estimated to more than 6 billion USD (McCown et al., 1999).This study was initiated by the flood in Tegucigalpa, the capital city of Honduras, on 30-31 October 1998 when Mitch struck the city.The estimated 500-year return period rainfall produced by Mitch (JICA, 2002) caused significant damage to Tegucigalpa, where 1000 casualties were reported and Published by Copernicus Publications on behalf of the European Geosciences Union.
approximately 40 % of its capital stock was damaged (Angel et al., 2004;JICA, 2002).In addition to these calamities, many of Honduras' hydrological archives were swept away from their premises at SANAA (Servicio Autónomo Nacional de Acueductos y Alcantarillados) which was located close to the main channel of the upper Choluteca River.
Simulations of water-level dynamics caused by disastrous events are needed for preparedness, to produce floodinundation maps useful for urban planning, and to prioritize investments (Pappenberger et al., 2006;Schanze, 2006).Such simulations are also relevant to better comprehend the hydraulic mechanism of large flood events in order to improve model structure (Beven et al., 2011;Jarrett, 1990).However, given that simulations of extreme floods are generally associated with limited data availability and large uncertainties, the question arises as to whether it is possible to achieve simulations that can be useful for contingency planning and prevention.
When hydrometric measurements of discharge and water levels during an event are lacking or highly inaccurate, such information may be inferred from post-event surveys.These can be done through eye-witness accounts and field campaigns (Brandimarte and Di Baldassarre, 2012;Ciervo et al., 2015;Gaume and Borga, 2008;Horritt et al., 2010;JICA, 2002;Smith et al., 2002), sometimes in combination with additional methods such as searches into historical documentation and paleo-flood techniques (Mård Karlsson et al., 2009;Smith et al., 2012;Valyrakis et al., 2015).Such surveys have been useful to estimate hydrometric data of the floods.Pictures and movies can be used to identify locations, flow type, depth, flow velocity, and discharge at the time they were taken (e.g.Ciervo et al., 2015;Le Boursicaud et al., 2016).Post-event information of channel topography and maximum water level can be used to estimate maximum peak discharge (Dalrymple and Benson, 1968;Matthai, 1968).
Post-event-estimated maximum peak discharge can be used to produce probabilistic regional envelope curves (Castellarin, 2007;Gaume et al., 2009) and discharge series for flood-frequency analysis (Coeur and Lang, 2008).These provide design-flood estimates used for inundation mapping (e.g.Brandimarte and Di Baldassarre, 2012).However, an assessment of flood development in time is required for early-warning systems (Schanze, 2006).The development of a flood in time can be obtained through a strategically planned post-event survey of peak discharge and the associated time of the peak (e.g.Delrieu et al., 2005).Detailed hydrographs can also be obtained from rainfall time series in conjunction with post-event hydrometric data, by the use of a rainfall-runoff model (RRM).A RRM in turn can be coupled with a hydraulic model to estimate the water-level development along a floodplain (Bonnifait et al., 2009;JICA, 2002;Montanari et al., 2009;Pappenberger et al., 2005a).Results from hydraulic models can be validated against postevent-estimated peak discharge, time of the peak, maximum water level, and flood-extent data (e.g.Bonnifait et al., 2009;Brandimarte and Di Baldassarre, 2012;Horritt et al., 2010).
Post-event data have been used with deterministic calibration within hydraulic models (e.g.Horritt et al., 2010;JICA, 2002), and for coupling RRMs with hydraulic models (e.g.Ciervo et al., 2015).Using post-event data, Bonnifait et al. (2009) present a multi-variable assessment to find a group of best parameter sets for the TOPMODEL RRM and a 1-D hydraulic model.Borga et al. (2008) and Pappenberger et al. (2006) suggest that post-event data should be used within an uncertainty-analysis framework given their large uncertainties.Di Baldassarre et al. (2010) discussed the advantages of distributed uncertainty mapping, as first proposed by Romanowicz and Beven (1998), in comparison with deterministic mapping.Uncertainty-analysis techniques have been used to account for uncertainty in hydraulic models (Aronica et al., 1998;Bozzi et al., 2015;Brandimarte and Di Baldassarre, 2012;Pappenberger et al., 2005aPappenberger et al., , 2007) ) and for the coupling of a RRM with a hydraulic model (Montanari et al., 2009;Pappenberger et al., 2005a) using eventmeasured data.
Uncertainty-analysis techniques account for possible errors involved in the modelling process, e.g.errors in model parameters and input data, due to lack of knowledge of their true values, spatio-temporal variability, or inaccurate estimation, and errors related to limited knowledge of the behaviour of the real system, i.e. epistemic uncertainty (Beven, 2016(Beven, , 2009)).Thus, in uncertainty-analysis techniques, uncertainties can be associated with several sources that interact among them, in which each interaction is associated with a likelihood dependent on how well it fits the observations.The formal Bayesian approach is a widely used method for uncertainty analysis, with different set-ups available (e.g.Smith and Roberts, 1993).Bayesian techniques have been commonly applied in hydraulic and hydrological modelling (e.g.Hall et al., 2011;Renard et al., 2008) and can be used within a global sensitivity analysis (see summaries by Iooss andLemaître, 2015, andSarrazin et al., 2016) to assess the effect of each source of uncertainty on the output (e.g.Abily et al., 2016).An informal Bayesian approach is the generalized likelihood uncertainty estimation (GLUE) framework (Beven and Binley, 1992), which differs in the way likelihood is defined and in that it does not require prior knowledge on the correlations or distributions of the parameter errors, yet with GLUE it is possible to get posterior information in the parameter combinations.In this study we hypothesize that it is possible to reasonably estimate, considering the large uncertainties in the observations, the extreme 1998 flood discharge in Tegucigalpa and the extent of the inundations that followed, from a combination of models and post-event data.We are aware of works that use the combination of hydraulic models and RRMs to assess flood dynamics or others that use post-event data to calibrate either RRMs or hydraulic models, both deterministic and through uncertainty analyses.We are not aware of any previous study combining a RRM, hydraulic modelling, and post-event data within an uncertainty-analysis framework to prove that reasonable estimation of an extreme flood is possible when hydrometric data are lacking.The methodology suggested in this paper integrates TOPMODEL (Beven and Kirkby, 1979;Kirkby, 1997), Muskingum-Cunge-Todini (MCT) (Todini, 2007) routing, and the LISFLOOD-FP (Neal et al., 2012a) hydraulic modelling tool in a GLUE framework.
2 Study area and data

Area description
The study area is the floodplain at Tegucigalpa, approximately 13 km of river length downstream from the upper part of the Choluteca River catchment.The area draining to the floodplain is around 811 km 2 and is composed of five sub-catchments: Grande River (448 km 2 ), Guacerique River (243 km 2 ), Chiquito River (71 km 2 ), Salada Creek (25 km 2 ), and Las Lomas Creek (12 km 2 ) (Fig. 1).Rainfall in the region is affected by high hurricane recurrence (Alvarado and Alfaro, 2003;Strobl, 2009) and convective activity.These two features in combination with the mountainous nature of the terrain (Amador et al., 2006) might lead to a high spatial variation of rainfall.Westerberg et al. (2010) found that daily precipitation has a high spatial variability and that bias in the estimations is likely due to insufficient gauge stations to measure in space and at different elevations.The land use and geology are relatively uniform in all sub-catchments.The land use is mainly composed of sparse coniferous for-est at higher-elevation lands, and fallow, pastures, and urbanized areas in the low land (CIAT, 2007).The geology at the surface is mainly composed of tuff and limestone to a minor degree; the superficial aquifer is classified as poor to moderately productive (ING, 1996).The average basin slope estimated in the Grande River, Guacerique River, Chiquito River, Salada Creek, and Las Lomas Creek sub-catchments is 19.5, 18, 25, 17.5, and 11 % respectively.Two reservoirs operated by SANAA are established upstream of the Tegucigalpa floodplain: Concepción reservoir, located at the Grande River sub-catchment, and Los Laureles reservoir, located at the Guacerique River sub-catchment (Fig. 1).

Topography
An airborne light-detection and ranging (lidar) survey in Tegucigalpa was conducted in 2000 by the University of Texas in cooperation with the US Geological Survey (USGS) during their survey in Honduras in response to Hurricane Mitch (Mastin, 2002).They generated a 1.5 m cell-resolution digital-terrain model (DTM) with an estimated vertical accuracy of 0.14 m (Fig. 2).These lidar data were used by Haile and Rientjes (2005) to investigate the effect of a digital elevation model (DEM) resolution on simulated flood extension using the SOBEK modelling tool.In 2001, JICA (2002) also conducted a topographic field survey as part of a flood/landslide-mitigation master plan and a total of 99 cross sections along the rivers in the floodplain at Tegucigalpa, surveyed at intervals of approximately 100 m, were used in this study (Fig. 2).In addition, orthographic pictures were taken at Tegucigalpa by JICA (2002).
The topography of the Tegucigalpa floodplain upstream sub-catchments was available from the 90 m spatial resolution Shuttle Radar Topography Mission (SRTM) data described by Reuter et al. (2007) (Fig. 1).

Precipitation
Upstream of the Tegucigalpa floodplain, two stations measured hourly rainfall during the Mitch event (Figs. 1 and 3).One of the stations is operated by Servicio Meteorológico Nacional (SMN, national weather service) and the other by the Universidad Nacional Autónoma de Honduras (UNAH).

Discharge
Discharge at three locations was estimated post-event by Smith et al. (2002) using the standard USGS techniques by Benson and Dalrymple (1967).The peaks at the Chiquito River and the Grande River (points 1 and 2 in Figs. 1 and 2 and Table 1) were estimated using the width-contraction analysis that uses the continuity and energy equations between a cross section approaching the contraction section under a bridge (Matthai, 1968).The peak at Choluteca (point 3) was estimated using the slope-area analysis, in which discharge is computed on the basis of the uniform-flow equation involving channel geometry, high-water marks, and roughness coefficients (Dalrymple and Benson, 1968).The measurements of discharge using the width-contraction analysis and the slope-area analysis can be associated with a 25 % error for unfavourable field-data conditions (Benson and Dalrymple, 1967), but up to 100 % overestimation might be associated with the slope-area analysis for slopes greater than 0.2 % (Jarrett, 1987).
A deterministic reproduction of the flood produced by Hurricane Mitch was done by JICA (2002) by setting a rainfall-runoff analysis using a linear reservoir model driven with hourly rainfall data from the SMN station.The produced hydrograph was used as input for the 1-D Mike 11 modelling tool (DHI, 2000) for unsteady flow conditions.In Hydrol.Earth Syst.Sci., 21, 3597-3618, 2017 www.hydrol-earth-syst-sci.net/21/3597/2017/ addition to the flood extent (Fig. 2), JICA (2002) reported the maximum peak discharge at the points 4, 5, and 7 in Figs. 1  and 2 and Table 1.Controlled flow release through the spillway at the Concepción reservoir was conducted and recorded by SANAA during the Mitch event (Fig. 3).The outflow over Los Laureles dam was not recorded.However, SANAA reported that its gate was overtopped at 22:30 LT on 30 October, reaching a maximum of approximately 1200 m 3 s −1 (JICA, 2002;Smith et al., 2002).Peak times in Table 1 except at point 5 were obtained by interviewing witnesses.The time of the peak at point 5 was estimated by propagating the peak reported at Los Laureles reservoir.

Maximum water levels
High-water marks during the Mitch flood were surveyed post-event by JICA (2002); the data were obtained by interviewing residents who experienced the event.The survey was carried out at the same locations where the topographic cross sections were made (Fig. 2).

Consistency in the post-event measured data
An inspection of the consistency of the data was done prior to the analysis.The inspection was done by plotting the maximum water-level profile to detect possible outliers.The consistency in timing and magnitude along the river network for the post-event maximum peak discharge was also checked.The flood-wave peak and time of the peak were expected to be larger and later downstream from the river confluences respectively.

Uncertainty and evaluation function
To quantify the propagation of uncertainty, the GLUE method was used.The assumptions of more formal statistical approaches can not be justified in data-scarce cases with high epistemic uncertainties.Within the GLUE methodology, pa- rameter sets were generated using a Monte Carlo technique, assuming a uniform prior distribution of the parameters.
Behavioural parameter sets, those that perform well in predicting the observations, were selected using a likelihood measure that reflected the performance of individual simulations with respect to one or several evaluation variables (o i ).Likelihoods were inferred by using the degree of belief (d i ) of a trapezoidal fuzzy membership function (Fig. 4), whose shape was chosen to account for uncertainties in the postevent estimated values, which are not considered crisp estimations.Thus the degree of belief for a difference smaller than a between the simulated and post-event estimated values is equal to 1, and it declines linearly to 0 for differences larger than b.

Modelling framework
The dynamic of the water level along the river channel and floodplain was reproduced with the sub-grid channel formulation of the LISFLOOD-FP hydrodynamic model (Neal Evaluation with post-event data: -maximum water level -maximum peak discharge and -time of the peak: Fuzzy measure within GLUE

Selection of class hydrographs:
K-means flat algorithm

Likelihood of inundated area during an extreme flood event
To consider rapidly varied flow (e.g. from a dam)

Representative hydrographs for the upstream boundary condition
Topographic information is a basis to set up TOPMODEL, which was one reason to select it in our mountainous catchment.Additionally, the version used here (Fuentes-Andino et al., 2017) has been shown to improve model prediction by considering the uncertainty associated with the spatially averaged estimation of rainfall.The mass-conservative version of the Muskingun-Cunge routing, the MCT, was incorporated to consider the sudden release of water from the Concepción reservoir, and it was chosen since a more complex routing could not be applied given the lack of data in the upstream area of the floodplain.The effect of Los Laureles dam on simulating the hydrograph of the Guacerique River sub-catchment was assumed to be negligible since the dam was overtopped long before the most intensive period of the storm.
The TOPMODEL and MCT combination assumes slopedependent variable velocity at a hillslope, constant velocity at a normal channel, and a variable velocity (according to the diffusive wave model of the MCT) at the main channel (whose length was estimated to have a minimum drainage area equal to 65 km 2 ).For each sub-catchment, the main channel was sub-divided into reaches of approximately 2.5 km to execute the MCT routing.For the MCT routing at Grande River, the inflow for the most upstream reach was set equal to the outflow hydrograph from the reservoir, and for other sub-catchments, to be equal to the hydrograph draining to that reach using TOPMODEL.For the subsequent reaches, this inflow was estimated as the sum of the outflow from the MCT routing at the immediate upstream reach and the hydrograph produced by TOPMODEL in the area draining to that reach (excluding the area draining to the upstream reaches).The modelling time step was equal to 5 min, smaller than the estimated travel time of the flood wave along the reach, as required by the MCT routing.
For the TOPMODEL, a network width function for each reach was created using topography from the SRTM raster.Only two rain-gauge stations were available, which made it difficult to infer the spatial distribution of rainfall.However, rainfall registered at the two stations was similar; thus, rainfall was assumed to be spatially uniform and estimated as the average of the two time series.Given the large magnitude of the event, it was expected to be associated with little spatial variation.
Uncertainty in rainfall input was taken into account by a multiplier (R); in addition, uncertainty of six model parameters was considered: the rate of decline of transmissivity (m), horizontal transmissivity (T o ), time constant (t d ), land-use coefficient (l u ), flood-wave celerity (v c ), and maximum soil infiltration rate (i max ) (Appendix A).The MCT method required information of the river slope and the geometry of the cross sections (Appendix B).The former was approximated from SRTM data, while the latter was inferred here as a function of discharge using the Manning equation for a wide parabolic channel as in Tewolde and Smithers (2007), with the channel roughness coefficient (n cu ) assumed uniform along all the reaches to make the modelling system simple and, in view of the lack of data, to constrain localized values.
All parameters were sampled from uniform distributions with ranges considered large but possible in the literature (Table 2) and each generated parameter set was used to simulate the Chiquito, Grande, and Guacerique river subcatchments (outlets at points 1, 2, and 5 in Figs. 1 and 2 and Table 1).A stopping criterion as in Pappenberger et al. (2005b) was used to decide the number of simulations required.For every 500 behavioural simulations added, a cumulative distribution function (CDF) of the predicted peak discharge and one of the time of the peak were estimated (evaluation variables; see Sect.3.4.1).These estimated CDFs were compared with the previous one and the number of runs was considered sufficient when the addition of behavioural simulations did not change the CDF significantly (i.e.P < 0.05) using the Kuiper (1960) statistical test (Appendix C).This statistical test was considered suitable since it is sensitive to changes in the tail and to the median values of the distribution; therefore, it makes sure that the distributions did not change along the whole range of values.
Las Lomas Creek and Salada Creek (points 8 and 9 in Figs. 1 and 2) did not have data to constrain the simulations and, by proximity, the behavioural parameters found at both Grande and Chiquito were used to simulate them.This is expected to not greatly affect the system as the contributing areas for Las Lomas Creek and Salada Creek are relatively small in comparison to the three sub-catchments where postevent data were available (Fig. 1).In addition, these two areas were smaller than the threshold drainage area for applying MCT; therefore, only parameters from TOPMODEL were transferred to those sub-catchments.

Output evaluation
To decide on behavioural hydrographs for the Chiquito, Guacerique, and Grande river sub-catchments, the maximum peak and time of peak post-event observations, together with their associated uncertainty, were used (refer to points 1, 2, and 5 in Table 1).The assumed uncertainty range was b = ( s m −1/3 ) 0.001-0.08 a = ±50 % of the peak flow for the peak magnitude and b = a = ±2.5 h for the time of peak (Fig. 4).For the evaluation of the hydrographs, a was set equal to b; thus, every hydrograph within the uncertainty bounds was considered behavioural and to have an equal degree of belief.The uncertainty in peak discharge at points 1 and 2 was chosen considering, and was assumed larger than, the value suggested in Benson and Dalrymple (1967).The discharge at point 5, although it was estimated by running a RRM by JICA ( 2002), was considered reliable for calibration since its magnitude was similar to the maximum peak outflow measured at Los Laureles dam, located in the same river and with nearly equally contributing upstream areas as in point 5 (Fig. 1).It was expected that 50 % of the flow uncertainty, as well as for points 1 and 2, was also reasonable at point 5.
All times of peak came from the same source, i.e. witness accounts, and there was no additional information on their uncertainties; thus, we allowed up to 2.5 h uncertainty considering that the survey was carried out 2 years after the event and because of the expected difficulties in witnesses identifying the exact times when the peak occurred.
To reduce computational costs and avoid redundancy, 100 representative hydrographs (class hydrographs) were obtained for each sub-catchment by clustering the full behavioural ensemble.Clustering was done using the K-means flat algorithm also called Lloyd's algorithm, originally developed by Lloyd (2006), described in Madhulatha (2012), and with a tool available for use at Mathworks (2011).Following the K-means algorithm, the number of groups (K) to cluster an ensemble of data (here the behavioural hydrographs) were defined (here equal to 100).Then, a number of K hydrographs were randomly chosen from the ensemble to represent the cluster centroids.Each of the hydrographs in the ensemble was assigned to one of the centroid hydrographs according to the smallest distance, here taken as the sum of the absolute differences between hydrographs.Subsequently the centroid for each of the clusters was replaced with the hydrograph whose sum of distances from all hydrographs within the cluster is minimized; and then each hydrograph in the ensemble was assigned to the new centroid found.The procedure of moving centroids and assigning hydrographs to new centroids is repeated until there is no change in the clusters.To consider the extreme cases, the hydrograph from each cluster with the largest sum of the distances to all other centroid hydrographs was chosen.

Flood-wave propagation
The LISFLOOD-FP was used to propagate the flood waves along the channels and across the floodplain.Here the subgrid channel formulation following Neal et al. (2012a) was used, where the floodplain and the channel have a 2-D square grid representation and flow is conveyed using the local inertia formulation (de Almeida et al., 2012).Thus, the continuity equation (Eq. 1) and a simplified version of the momentum equation (where the convective-acceleration term was assumed negligible) (Eq.2) were used to keep the continuity of mass and momentum respectively in each cell and between cells.
where Q x and Q y , and S x and S y are the volumetric flow rates and the slopes respectively in the x and y directions, h the water depth, t time, A the cross-sectional area of flow, g gravity, n the Manning coefficient, and R the hydraulic radius, taken as the cell cross-sectional area divided by the wetted perimeter.
Equations ( 1) and (2) are solved using an explicit forward difference scheme on a staggered grid (Bates et al., 2010) which requires fewer numerical operations (about an order of magnitude) than a full 2-D dynamic model (Neal et al., 2012b).The former numerical procedure was computationally more efficient than the latter and therefore more suitable for uncertainty analysis.In addition, the model-grid representation made it possible to obtain the discharge and water-Hydrol.Earth Syst.Sci., 21,2017 www.hydrol-earth-syst-sci.net/21/3597/2017/ level time series output at any grid along the channel or floodplain.
The basic input data for the LISFLOOD-FP are topography, hydrographs at the upstream boundary conditions, a downstream boundary condition, and Manning roughness coefficients.To use as topographic input to the model, the lidar data were aggregated to 21 m cell resolution, as a tradeoff between high resolution and the speed of simulations.The surveyed cross sections and orthographic pictures from JICA (2002) were used to define channel depth and width respectively.Test simulations of this event were performed within the HEC-RAS 1-D hydraulic model (Brunner, 2001) considering the topography of the bridges, and preliminary results showed that bridges had a negligible effect on the overall flood profile.Thus, the geometry of bridges in the LISFLOOD-FP implementation was neglected by assuming a limited and localized impact on flood levels as in e.g.Castellarin et al. (2009), especially since the calibration data are associated with large uncertainties so that the localized effect of structures is not possible to detect (Fewtrell et al., 2011).
Uncertainty of the input hydrographs at each of the upstream boundary conditions was considered by sampling from the 100 class hydrographs.By assuming normal flow, the overall downstream valley slope, b c , was used as the downstream boundary condition.This assumption was considered in view of the lack of hydrograph information at the downstream boundary; however, water-level predictions at the most downstream cross sections can be associated with larger uncertainties due to this assumption (Pappenberger et al., 2006).Besides b c the channel-roughness coefficient, assumed uniform along all the channel length, n c , and the floodplain-roughness coefficient uniform along all the floodplain, n f , were also considered to be uncertain parameters.Ideally roughness coefficients would be allowed to vary spatially to reflect changes in channel and floodplain characteristics (e.g. one value per reach or per each side of the floodplain), but this would have led to an increased number of parameters in the hydraulic model when there was not enough information at each reach to constrain the local roughness.
Using a one-at-a-time (OAT) design for sensitivity analysis, the effect that uncertainty in the channel depth and channel width (through a multiplying factor) had on the outputs was explored, which led to the incorporation of the channelwidth multiplier (w f ) in the uncertainty analysis.For the hydraulic simulations, a total of 130 000 parameter sets were sampled from a uniform distribution with ranges considered large but possible in the literature (Table 3) in the same way as for the RRM.

Output evaluation
Different degrees of belief (d i ) (Fig. 4) were obtained by comparing the simulations with the following evaluation data.
-One degree of belief value, d 1 , as performance in predicting the maximum peak discharge value of point 3 (Figs. 1 and 2 and Table 1).
-Two degrees of belief values, d 2−3 , as performance in predicting the time of the maximum peak discharge of points 3 and 6 (Figs. 1 and 2 and Table 1).
-Ninety-nine degrees of belief values, d 4−102 , as performance in predicting maximum water levels along the main river and two tributaries (Fig. 2).
The fuzzy set values of a and b for evaluating the simulated peak discharge were set to 20 and 50 % of observed values respectively.Thus, for differences between observed and predicted peak discharge within 20 % of the observation, the degree of belief was assumed to be equal to 1, and decreased to 0 for differences larger than 50 %.These values were chosen taking into account those values suggested by Benson and Dalrymple (1967) and Jarrett (1987).The fuzzy set values of a and b for evaluating the time of the peak were set equal to 0.5 and 2.5 h respectively; thus, the degree of belief for differences between observed and predicted times of the peak smaller than 0.5 h was assumed to be equal to 1, and it decreased to 0 to allow for up to 2.5 h of difference, an error considered possible in the observations.And finally, the fuzzy set values of a and b for evaluating the water levels were set equal to 0.5 and 1.8 m respectively.Thus, 0.5 m was chosen to account for error in topography representation (Neal et al., 2009), and 1.8 m was chosen considering the magnitude of the observed water level and that 2 years after the event witnesses' memories might have been associated with large uncertainties.A parameter set was considered behavioural if the degree of belief was larger than 0 for each of the 102 evaluation points.For every parameter set, a global score (GS) was calculated based on a weighted average of the degrees of belief obtained for each evaluation criterion.

GS
where w i are the weights associated with the degrees of belief corresponding to the observations.The weight associated with the peak discharge and the two times of the peak data (d 1−3 ) were set equal to 0.1 each; thus, 0.7 was the weight corresponding to the sum of the degrees of belief associated with all the observed maximum water levels (d 4−102 ).A larger aggregated weight was given to predict the observed water marks in comparison to the peak discharge and times of the peaks to reflect the larger number of observed water marks (99) and because the focus was on predicting flood extent.The weights could be changed according to the purpose of the study, which might also result in different ensembles being behavioural for different purposes (Pappenberger et al., 2007).Subsequently, likelihood values were obtained by scaling the global scores by a constant C, so they will sum to  unity over all behavioural sets (Beven, 2009).Finally, the behavioural parameter sets were used to generate a fuzzy likelihood water-level profile and map of the maximum flood extension during the Mitch event as in Di Baldassarre et al. ( 2010).

Consistency in the post-event measured data
From prior inspection of the data, it was found that information about the maximum peak discharge and time of the peak was consistent (i.e. in comparison to locations at the upstream reaches): discharge values and times of the peaks were larger and later at downstream locations after the confluences.A plot of the high-water marks showed sudden jumps at some observation points without any obvious physical explanation, but this is perhaps to be expected given the origin of those observations (witness accounts from memory).Thus, we did not eliminate any of the observations, but instead allowed an uncertainty range associated with all observation points.

Representative hydrographs for the upstream boundary condition
Behavioural hydrographs to use as the upstream boundary conditions of the hydraulic model were obtained for the subcatchments of the Grande, Guacerique, and Chiquito rivers and, by using behavioural sets at the Grande and Chiquito river sub-catchments, at the Salada Creek and Las Lomas Creek sub-catchments (Fig. 6).The cumulative distribution function (CDF) of the predicted peak discharge and of the time of the peak of 2000, 8000, and 9000 behavioural simulation for sub-catchments of the Chiquito, Guacerique, and Grande rivers respectively did not change significantly by adding 500 behavioural simulations more.Thus a total of 3000, 9000, and 10 000 behavioural simulations, obtained from a total of 61 205, 60 237, and 60 833 samples respectively, were considered enough to infer 100 class hydrographs for the Chiquito, Guacerique, and Grande river subcatchments respectively.When comparing the prior and pos-   terior distributions of the rainfall-runoff model parameters, five out of eight parameters were sensitive: the rainfall multiplier (R), rate of depletion (m), time constant (t d ), main channel roughness coefficient (n cu ), and maximum soil infiltration rate (i max ) (Fig. 7).

Flood-wave propagation
There were no simulations for which all degrees of belief were larger than 0. Criteria d 1−3 were fulfilled by 47 894 out of 130 000 total simulations, but some observed water marks (criteria d 4−102 ) were constantly and largely under-or overpredicted.To allow for special cases, i.e. larger error in the observations or in the hydraulic simulations, the constraints  were relaxed by allowing 10 % of observed water marks (10 out of 99 observations) to be outside the fuzzy bounds, i.e. the degree of belief was allowed to be equal to 0. By relaxing the constraints a total of 6357 parameter sets were found; the degrees of belief for those parameters varied between 0.001 and 1, 0. Change in the posterior distributions of the parameters showed that the channel roughness coefficient and floodplain roughness coefficient were more sensitive than the channel width factor and the slope for the downstream boundary condition (Fig. 8).Changes in the posterior distribution of the peak and time of the peak showed that the model was unsurprisingly more sensitive to input hydrographs from the larger sub-catchments than from small sub-catchments (Fig. 9).Flood-wave propagation of different input-hydrograph combinations led to prediction of two markedly different times of the peak at the floodplain, resulting in under-(over-) prediction when the earliest (latest) peak of input hydrograph combinations prevailed (Fig. 10).
There were three observed high-water marks in the Chiquito River reach that were constantly under-predicted and outside the uncertainty bounds of the observations (Fig. 11).The propagation from the water-level uncertainty to the flood extent was more evident in urban areas, where the flood extent varies more with changes in the water level due to the presence of structural features such as buildings (Fig. 12).From behavioural simulations, the 90 % confidence interval for prediction of the discharge at the floodplain outlet was 2708 to 4619 m 3 s −1 , encompassing the 3880 m 3 s −1 value estimated in JICA (2002) (reference point 7 in Fig. 1 and Table 1).For reference point 4, at the Chiquito River, the 90 % confidence interval was 247 to 482 m 3 s −1 , also encompassing the 436 m 3 s −1 value estimated in JICA (2002).

Discussion
A field campaign after a large flood event is an opportunity to collect information useful for flood forecasting and subsequent contingency planning in places where hydrometric measurements are lacking because of non-existing or broken gauges.
Our study demonstrated that it was possible, in a datascarce situation, to reproduce an extreme flood event that was within the bounds of the uncertainty in the evaluation data.Our results support those of Bonnifait et al. (2009) and Ciervo et al. (2015) about the possibility of reproducing an extreme flood event by a suitable combination of RRM and hydraulic modelling tools with only event-based rainfall data and post-event hydrometric data.Here we additionally incorporated the GLUE methodology to account for expert knowledge of uncertainties in model parameters, rainfall input, and evaluation data.Thus, the combination of a RRM with a hydraulic modelling tool within an uncertainty framework as in Montanari et al. (2009) and Pappenberger et al. (2005a) proved to be useful also in the case with only post-eventestimated hydrometric data.
After considering the uncertainties and their interaction it was possible to identify behavioural parameter sets that were used to obtain a realistic probabilistic reproduction of the flood-water level (Fig. 11) and flood extension (Fig. 12).In comparison to the deterministic estimates made by JICA (2002) using different modelling tools, in this work it was possible to obtain predictive ranges of the water level that encompassed most of the observations.The flood extent here, associated with a likelihood at each flooded cell, generally extended beyond the extent of the JICA (2002) mapping.
The combination of TOPMODEL and MCT allowed us to estimate behavioural hydrographs for the Chiquito, Guacerique, and Grande sub-catchments.The simulations could be constrained (Fig. 6) in spite of the wide uncertainties in the data and the simplified assumption of the MCT routing for ungauged basins applied here.The rainfall multiplier (R), rate of depletion (m), time constant (t d ), main channel roughness coefficient (n cu ), and maximum soil infiltration rate (i max ) were more important in selecting the resulting hydrographs (Fig. 7), whereas horizontal transmissiv-ity (T o ), land-use coefficient (l u ), and flood-wave celerity (v c ) were less sensitive.
The rainfall multipliers were sensitive and the means of their posterior distributions varied across sub-catchments (0.93, 1.5, and 1.3 for Chiquito, Guacerique, and Grande respectively) (Fig. 7), suggesting that the spatial average rainfall estimated from the two available gauges was overestimated at Chiquito and underestimated at the Guacerique and Grande sub-catchments.The Guacerique and Grande subcatchments are larger and have a higher topographic elevation than the Chiquito sub-catchment.Underestimation of rainfall for these sub-catchments might be the result of lack of stations to represent the rainfall spatial pattern, highly variable in the area (Westerberg et al., 2010).Thus, a simplistic account of a space-and time-averaged rainfall multiplier as in Fuentes-Andino et al. (2017) was also useful here to account for bias estimation of the spatially averaged rainfall.The posterior distribution of the rainfall multiplier at the Chiquito and Guacerique sub-catchments clearly aggregated to different mean values.The sensitivity to the multiplier was different in the case of the Grande sub-catchment, which also showed a different posterior marginal distribution shape for the rate of depletion (m) and the time constant (t d ) (Fig. 7).
Different shapes of posterior marginal parameter distributions at the Grande River sub-catchment relative to the Guacerique and Chiquito river sub-catchments could be caused by parameter adjustment to fit the observations or by different hydrological processes going on in the different sub-catchments.The sudden release of water from the dam could also be a reason for these differences.The posterior marginal parameter distributions for the Grande River sub-catchment suggest that it has a shallower effective soil depth (low m) and a faster channel response in the MCT routing (low n cu ) than the other two sub-catchments.Hydrographs from a total of five sub-catchments (Fig. 6) from the TOPMODEL and MCT combination were used as upstream boundary conditions for the hydraulic simulations.
Even if more detailed post-event observations of flood extent might do better than water levels in constraining the LISFLOOD-FP (Fewtrell et al., 2011;Horritt and Bates, 2002), the modelling tool predicted the observed high-water marks, peaks, and times of peaks well.Behavioural simulations for which the degree of belief for the peak discharge, time of the peak, and at least 90 % of predicted high-water marks (89 out of 99 observations) were above 0 were identified.
The channel and floodplain roughness coefficients were the most important parameters for the hydraulic model (Fig. 8).As roughness coefficients directly affect the estimation of discharge and water level, the impact of their uncertainty has been shown previously in other studies (Dimitriadis et al., 2016;Pappenberger et al., 2005b;Warmink and Booij, 2015;Wohl, 1998).Here, uncertainty is expected to be particularly large as these coefficients interacted with uncertain post-event estimated discharge and high-water marks and also because they were assumed to be spatially aggregated due to data limitations.For example, a more localized calibration of such coefficients could have helped to tackle the problem of localized channel erosion during flood events common in the area (Guerrero et al., 2012).Given the assumed spatial representation of the roughness coefficients and the uncertainty they are associated with, they interacted with all other sources of uncertainty in a complex way that is difficult to separate.Such complex interactions are contained implicitly in the resulting ensemble of behavioural simulations (Beven, 2016).
The effect of the input hydrographs from the Grande River and Guacerique river sub-catchments on the resulting outputs is evident in Fig. 9. Thus, as in Dimitriadis et al. (2016), here the roughness coefficients and input flow were the most important sources of uncertainties.Two peaks in the input rainfall (Fig. 3) led to two main large peaks in the hydrographs as input boundary conditions (Figs. 6 and 9).The propagation of input hydrographs along the floodplain led to under-or over-prediction of the times of peak (Fig. 10).This suggests that the spatial pattern of rainfall was not well represented by the gauge average, as also suggested by the posterior distribution of rainfall multipliers in the RRM.Since rainfall data played an important role in predicting the times of peak, investment to improve the rainfall measurement system, e.g.radar estimates or a denser rain-gauge network, should be prioritized in the study area, especially because these data are easier to collect relative to discharge in a high-magnitude event.
Some observed high-water marks were constantly largely underpredicted in the estimates by JICA (2002) and outside the prediction bounds produced here, even when allowing for significant uncertainty in the evaluation data (Fig. 11).Inspection at the points that were constantly underpredicted showed that no man-made structure could have been the reason for such disagreement.Thus the problem of predicting at those locations could be caused by the inability of the hydraulic modelling tool to simulate the system under extreme conditions where effects such as sharp river bends might have an important local effect on the flow.However, a previous experiment using the 1-D HEC-RAS model on the same river also agreed with the results obtained here, and no localized effect in the under-predicted places was obtained.Another reason for the disagreement could be large errors in the postevent data.
In general, minor errors between prediction and observations in this work could be caused by a weak spatial representation of topography and roughness coefficient, i.e. special topographic details in a highly populated area with manmade structures that could not be captured by the DEM.However, those local features might not affect the general flood extent (Haile and Rientjes, 2005).
The peak discharge at point 3 (Figs. 1 and 2) was underpredicted by most of the simulations (Fig. 10).However, the high-water mark was over-estimated at that location Hydrol.Earth Syst.Sci., 21, 3597-3618, 2017 www.hydrol-earth-syst-sci.net/21/3597/2017/ (Fig. 11).The reasons for this could be an over-estimation of the post-event peak discharge, or an under-estimation of the observed high-water mark, or the simplistic representation of the downstream boundary condition assumed.
A general under-prediction of the water level in the Chiquito River reach could be due to the low (perhaps under-estimated) post-event-estimated peak discharge; as in the comparison with the Grande and Guacerique subcatchments, most of the hydrograph simulations for the former were rejected because the simulated peaks were larger than the evaluations (even considering the uncertainty) (Fig. 6).This could also be the reason for a lower rate of behavioural sets for the Chiquito River sub-catchment when comparing with the other two.
A detailed inspection of model structure, model set-up, and data at specific points where the modelling tools did not perform well even after considering possible uncertainties in the parameters, input, and evaluation data, could reveal areas for improvement.
This study was set up to demonstrate the use of post-event data and a combination of suitable RRM and hydraulic modelling tools with uncertainty analysis to reproduce an extreme flood in a data-scarce area.The behavioural ensemble found here depends on the uncertainties coming from the model structure (Dimitriadis et al., 2016), quality of the data (Pappenberger et al., 2006), topographic resolution (Haile and Rientjes, 2005), and spatial aggregation of the parameters (Beven, 1995).Considering the dependency with those sources of uncertainties and their interaction, the post-event data proved to be useful in reproducing the Hurricane Mitch flood event.High-water marks obtained from personal memories of an event are a good source of information.To decrease uncertainty of such information, institutions in charge of disaster prevention should be prepared to carry out such surveys soon after flood events when memory is fresh.In fact, soon after extreme events it is also possible to collect that information by surveying the marks left by the flood (e.g.Neal et al., 2009).Post-event-estimated peak discharge, though it is known to be associated with large uncertainties (Benson and Dalrymple, 1967;Jarrett, 1987), was a valuable source of information in this work.A higher spatial availability of flood peak discharge and time of the peak estimates would greatly benefit this methodology as it will allow a better quality control of individual estimates, to leave some of the estimates out for validation, and to estimate more localized patterns of roughness coefficients.
The use of this methodology can be done within a Bayesian framework in which the posterior distribution of the parameters is updated when more events become available.Data from more events could further reduce the predictive uncertainties and help us to learn from the flow behaviour at some localized areas where the errors were large.Post-event estimates in the future could likely also come from social-media information which is becoming gradually more available (Fraternali et al., 2012;Triglav-Čekada and Radovan, 2013).
The flood-hazard map presented here can be used by the committee in charge of disaster contingency and management in the city of Tegucigalpa (CODEM-DC) as a complement to the 5-, 10-, 25-, and 50-year return period hazard map produced in JICA ( 2002) and the the 50-year hazard map produced in Mastin (2002) and Mastin and Olsen (2002) for spatial planning and to prioritize investment.If real-time discharge measurements are available to calculate the initial saturation of a catchment, behavioural parameter sets updated from a range of events can be used for forecasting the flood extent as shown by Romanowicz and Beven (2003) and Montanari et al. (2009).In the absence of such measurements, a guess of the initial discharge may also work since it will not significantly affect the prediction for the intense period of the event.Furthermore, for that period, our methodology can give a better performance since calibration is done against discharge, time, and water level at the peak.It is also tempting to consider this methodology for forecasting fed both by an improved rain-gauge network and water-level information coming from social media.

Conclusions
In this study we tested the possibility of reproducing an extreme flood disaster in a data-scarce area, the devastating flood in Tegucigalpa triggered by Hurricane Mitch in 1998.It was possible to realistically reproduce this large ungauged flood event by using post-event hydrometric data in combination with rainfall data and various modelling tools, demonstrating the value of post-event field campaigns to constrain the uncertainties in estimates of hydrometric data, model parameters, and output.A methodology has been proposed where post-event-estimated data are used to drive and constrain a combination of rainfall-runoff and hydraulic modelling tools to reproduce floods within a GLUE uncertaintyanalysis framework.Results of the flood extent proposed here were comparable to the deterministic mapping produced by JICA (2002) using different modelling tools.However, here more information was embedded as likelihoods of inundation associated with each cell in the floodplain.
Combining the TOPMODEL with the MCT routing to reproduce hydrographs in catchments with rapidly varied flow, e.g. release from a dam, resulted in hydrographs that were within the uncertain bounds of the observations.The predictive capability of the TOPMODEL and MCT combination warrants further exploration with more detailed and less uncertain event data.The rate and bias in the rejection of the hydrographs due to over-estimation indicated underestimation of post-event-estimated discharge at one location.The propagation of estimated hydrographs through the 2-D hydraulic LISFLOOD-FP resulted in successful predictions of observed high-water marks, discharge peaks, and times of peaks within the uncertainty bounds for most of the evaluation variables.A few critical locations in the floodplain were identified where the model set-up could not reproduce the maximum water level.Locations of disagreement between simulations and evaluations, after considering all important sources of uncertainties, can provide information useful for improving model structure or post-event data-estimation methods.Results showed the importance that rainfall data have in simulating the times of peaks; thus, results would be improved by a better spatial assessment of rainfall.Improvements of this methodology can be done by using it within a Bayesian framework of updating the parameters' posterior distribution when more events become available.The methodology proposed here can be useful for planning, prioritizing investments, and flood forecasting.
Data availability.The lidar data on Tegucigalpa described in Mastin (2002) are available on request to the US Geological Survey (USGS).Post-event-estimated peak discharge obtained by the USGS is found in Smith et al. (2002).Post-event-estimated peak discharge obtained by JICA ( 2002) is found in http://libopac.jica.go.jp/images/report/P0000054205.html.The topography field survey and maximum water level survey described in JICA ( 2002) are found in http://libopac.jica.go.jp/images/report/P0000054208.htmldata provided by the Servicio Meteorológico Nacional (SMN, national weather service) and by the Universidad Nacional Autónoma de Honduras (UNAH), found in the Supplement.Data of the flow release at the Concepción reservoir provided by the Servicio Autónomo Nacional de Acueductos y Alcantarillados (SANAA) can be found in the Supplement.
where P τ = 4.75 √ Q τ is the wetted perimeter estimated for stable river channels, S the reach slope, and n the Manning roughness coefficient.
where a coefficient equal to 1.4 was chosen as the average between a parabolic channel and wide rectangular channel (1.2 and 1.6 respectively).
where W τ is the top flow width, assumed to be approximately equal to the wetted perimeter (P τ ).The specialization factor for correction of the Courant and Reynolds number, β τ , following Todini (2007), is Thus, the corrected Courant number, C * τ , is estimated as and the corrected Reynolds number, D * τ , as which yields the following MCT parameters:

Figure 1 .
Figure 1.Study area and data location; topography data from the Shuttle Radar Topography Mission (SRTM).

Figure 2 .
Figure 2. Geometry set-up for hydraulic simulation at the Tegucigalpa floodplain.Lidar data from Mastin (2002).

Figure 4 .
Figure 4. Fuzzy membership function for evaluation of model performance: a and b depend on the uncertainty associated with the evaluation (o i ).

Figure 5 .
Figure5.Scheme of the modelling framework used to reproduce an extreme flood event using post-event-estimated data to drive and constrain a combination of modelling tools within an uncertaintyanalysis framework.

Figure 6 .
Figure 6.Precipitation (bars) and 100 class hydrographs chosen from the behavioural ones (black plots) for five sub-catchments upstream of the floodplain.Predictive range of the 100% probability limits for all hydrograph simulations (grey shaded area) and rectangles representing the fuzzy set to allow for uncertainty for peak discharge and time of the peak for the sub-catchments of the Chiquito, Guacerique, and Grande rivers.

Figure 7 .
Figure 7. Prior (grey) and posterior (black outlined) relative frequency distribution for the most sensitive rainfall-runoff parameters: rainfall multiplier (R), rate of depletion (m), time factor (t d ), the main channel roughness coefficient (n cu ) and the maximum soil infiltration rate ( imax ) for the Chiquito, Guacerique, and Grande catchments (first, second, and third rows respectively).

Figure 8 .
Figure 8. Prior and posterior relative frequency distribution (grey and black outlined bars respectively) of the LISFLOOD-FP parameters (width factor, slope for the downstream boundary condition, channel roughness coefficient, and floodplain roughness coefficient: w f , b c , n c , and n f respectively).

Figure 9 .
Figure 9. Prior and posterior relative frequency distribution (grey and black outlined bars respectively) of simulated maximum peak and time of the peak of input hydrographs for boundary conditions.

Figure 10 .
Figure10.Performance of the model in predicting high-water marks, average (d 4−102 ), against predicted maximum peak discharge and two times of peak at the Choluteca River (reference points 3 and 6 at Figs.1 and 2and Table1) for non-behavioural simulations (grey dots) and behavioural ones (black dots).Observed values and their limits of acceptability are plotted in continuous and dashed vertical lines respectively.

Figure 11 .
Figure 11.Likelihood of high-water marks during the Mitch event, considering uncertainty in model parameters, model input, and evaluation data to drive and constrain a combination of rainfall-runoff and hydraulic modelling tools.

Figure 12 .
Figure 12.Likelihood of inundated area during the Mitch event on 30-31 October 1998, considering uncertainty in model parameters, model input, and evaluation data to drive and constrain a combination of rainfall-runoff and hydraulic modelling tools.The deterministic flood extent was obtained by digitalization of the flood extent in JICA (2002).

Table 1 .
Post-event estimated peak discharge and time of peaks.

Table 2 .
Sampling parameter ranges to run the rainfall-runoff model.

Table 3 .
Sampling range of parameters to run the hydraulic model.
Ôt+ t = C 1 I t+ t + C 2 I t + C 3 O t .(B12)All the estimations for the time τ = t + t are computed twice to eliminate the influence of the first guess Ôt+ t in Eq. (B1).