Coupling a global hydrodynamic algorithm and a regional hydrological model for large-scale flood inundation simulations

In recent years, a series of destructive flood events affected Europe and raised public and scientific interest in the risk of floods. As flood risk is defined as the product of flood hazard, vulnerability and exposure (Kron 2005), the simulation of the extent and depths of floods is a key component for flood risk analysis.

Flood inundation models provide a valuable tool for flood hazard mapping and flood risk analysis. They vary with regard to their different dimensionality and one dimensional (1D) models, 2D models and a coupling of 1D/2D models are widely used in literature (Hunter et al. 2007; Tayefi et al. 2007; Chatterjee et al. 2008; Liu et al. 2015; Vozinaki et al. 2016). Generally, the high dimensional models (e.g. 2D models) are limited to small-scale applications (e.g. Pappenberger et al. 2008) or specific river reaches (e.g. Veijalainen et al. 2010). In addition, the long duration of floods and the need to represent precipitation and evaporation processes make it difficult to apply these models for large-scale studies (Wilson et al. 2007). However, larger scale is of more interest for integrated flood management and for climate impact studies than smaller scale studies. For example, the European Directive on Assessment and Management of Flood Risks requires that flood hazard and risk assessments are carried out at the basin scale (European Commision 2007), which has encouraged modeling of the extent of floods in large river basins in Europe.

Despite the difficulties, various modeling approaches have been explored for large river basins over the last 10 years. Due to the low computational cost, the fully 1D river hydrodynamic models with simple representation of floodplain/polders and the quasi-2D models are favored in large river basins (Castellarin et al. 2011; Paiva et al. 2011; de Paiva et al. 2013; Wicks et al. 2013). Another widely used approach is coupling a 1D hydrodynamic model and a 2D raster-based model, for example, LISFLOOD-FP (Wilson et al. 2007; Biancamaria et al. 2009; Trigg et al. 2009). da Paz et al. (2011) developed the LISFLOOD-FP model further by considering the precipitation, evaporation and infiltration processes over the floodplain. However, to avoid excessive computing time, the raster resolution is relatively coarse in large-scale applications (ca. 0.3–2 km).

In order to obtain a fine resolution flood hazard map for the whole of Europe, Alfieri et al. (2014) downscaled the estimated 100-year floods to the river network at a resolution of 100 m and ran small-scale floodplain hydraulic simulations for every 5 km river reach using a 2D hydraulic model (LISFLOOD-ACC). The final flood hazard map is a mosaic of output maps for all river reaches. However, such a reach-based approach cannot consider the attenuation of a flood wave passing along a long and ramified reach.

Recently, Falter et al. (2016) presented a continuous modeling chain for flood risk assessment in the German Elbe basin at a resolution of 100 m. This modeling chain includes a rainfall-runoff model, a 1D river network, 2D hinterland inundation and damage estimation models. The authors found that the hydrodynamic simulations generate significant uncertainties due to low data quality and dyke breach effects. Therefore, it remains challenging to precisely quantify the flood extent for large rivers even though the computation capability has been improved significantly.

Besides the regional hydraulic models, global river routing models were improved to include components of floodplain storage in recent years. These models were typically designed to convey runoff generated by global hydrological or land-surface models to the sea (e.g. Oki & Sud 1998), and they usually use relatively coarse resolution grids (0.25–0.5°). In order to represent detailed river hydrodynamics and inundated areas within such coarse grids, recent advanced global river routing models take into account sub-grid topographic information derived from high-resolution digital elevation models (DEMs) (Yamazaki et al. 2011; Decharme et al. 2012). Since both regional hydraulic models and the global river routing models were applied for large-scale floodplain simulations, Neal et al. (2012) introduced the sub-grid scale approach into the 2D model LISFLOOD-FP and aimed to compare it with the conventional model structures (1D model, 1D/2D model and 2D model) for an 800 km reach of the Niger River. The results show that the sub-grid scale 2D model increased simulation accuracy in terms of water level, wave propagation speed and inundation extent due to its ability to represent any size of river channels below the grid resolution.

However, as with the input data for the global river routing models, the outputs of global hydrological or land-surface models are often biased due to errors in inputs, in particular rainfall, and uncertainties in the parameterization of dominant hydrological processes (Haddeland et al. 2011). As a result, flood inundation maps based on direct global model outputs cause large uncertainty for regional decision-making and a downscaling procedure of the global model results needs to be applied prior to the local-scale inundation estimation (Winsemius et al. 2013).

Within the second phase of the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP2) (https://www.isimip.org/), an ensemble of regional hydrological models were successfully applied to 12 large river basins worldwide and they generally show better results than the global hydrological or land-surface models in terms of discharge amount and seasonal dynamics (Hattermann et al. 2017). Since the input data are mostly the same for both regional and global models in this project, one important reason for the poorer performance of global models is that they were not calibrated specifically for each river basin, in contrast to regional hydrological models. Hence, one can expect a more reliable inundation estimation based on the runoff simulated by the regional models than by the global ones for large river basins.

This study aims to propose a new coupling of a regional hydrological model and a global river routing model for flood risk assessments in large river basins. Such coupling is especially important because it can be used for long-term ensemble simulations (e.g. ensembles of climate change scenarios) at large scales, while the computational cost of regional hydraulic models can be very expensive for this purpose. The coupled models were tested in a meso-scale catchment as a proof-of-concept study. It takes the advantages of the sub-grid feature of the global routing method and the outputs from regional hydrological models. In addition, this study intends to fill the gap between the widely used regional 1D/2D hydraulic models and the global river routing methods for large-scale studies.

In this study, we coupled the process-based ecohydrological model SWIM (Soil and Water Integrated Model) (Krysanova et al. 1998) and the global river routing model CaMa-Flood (Catchment-based Macro-scale Floodplain model) (Yamazaki et al. 2011). Both models were applied for large-scale river basins with satisfactory results in previous studies. For example, the SWIM model showed a good performance for several large rivers worldwide in the historical period (1971–2000) within the framework of the ISIMIP project (Huang et al. 2016). The CaMa-Flood model generated promising simulated discharges across the globe and showed a good agreement between the simulated and the satellite-observed flooded area in the Amazon basin (Yamazaki et al. 2012).

The SWIM model was developed on the basis of the models SWAT (Arnold et al. 1993) and MATSALU (Krysanova et al. 1989). It simulates both discharge and runoff with a daily time step.

SWIM simulates all processes by disaggregating the basins to sub-basins and hydrotopes, where the hydrotopes are sets of elementary units in a sub-basin with homogeneous soil and land-use types. It is assumed that a hydrotope behaves uniformly regarding hydrological processes.

Water flows are calculated for every hydrotope and then lateral fluxes of water to the river network are simulated taking into account retention. After reaching the river system, water is routed along the river network to the outlet of the simulated basin using the Muskingum flow routing method (Maidment 1993).

The simulated hydrological system consists of three control volumes: the soil surface, the root zone of soil and the shallow aquifer. The water balance for the soil surface and soil column includes precipitation, surface runoff, evapotranspiration, subsurface runoff and percolation. The water balance for the shallow aquifer includes groundwater recharge, capillary rise to the soil profile, lateral flow, and percolation to the deep aquifer.

Surface runoff is estimated as a non-linear function of precipitation and a retention coefficient, which depends on soil water content, land use and soil type (modification of the Soil Conservation Service (SCS) curve number method (Arnold et al. 1990)). Lateral subsurface flow (or interflow) is calculated simultaneously with percolation. It appears when the storage in any soil layer exceeds field capacity after percolation. Potential evapotranspiration is estimated using the method of Priestley-Taylor (Priestley & Taylor 1972). Actual evaporation from soil and actual transpiration by plants are calculated separately.

There are 12 parameters used for calibration: two parameters of the Muskingum routing method, three parameters of the degree-day snowmelt process, two parameters as the temperature and precipitation laps rates, two parameters of the potential evapotranspiration method, two parameters controlling the subsurface and groundwater flow contribution to the streamflow and one correction parameter for the saturated conductivity of the soil layers.

The CaMa-Flood model is a distributed global river routing model, which routes runoff generated by land-surface models (Yamazaki et al. 2011). It calculates river and floodplain water storage, river discharge, water depth and inundated area at grid bases. Each grid has a river channel reservoir and a floodplain reservoir, which corresponds to the unit catchment of the river channel for that grid rather than the grid itself (see Figures 1 and 2(b) in Yamazaki et al. 2011). The river channel reservoir has three parameters, channel length, channel width and bank height, while the floodplain reservoir has two parameters, unit catchment area and a floodplain elevation profile. The floodplain elevation profile describes the floodplain water depth as a function of floodplain area (see Figure 2(c) in Yamazaki et al. 2011). Note that this simplified function ignored all topographic depressions including permanent lakes and wetlands. The output of this routing model is water storage, water level, flood areas, and routed discharge within each grid cell at daily time steps.

River discharge for each grid is calculated by a diffusive wave equation, which is a simplification of the fully one-dimensional St Venant momentum equation without the acceleration and advection terms. The default calculation interval is 5 min but it can be reduced automatically for preserving water balances and the numerical stability.

There are only two calibration parameters for the CaMa-Flood model: the Manning's values for river bed and floodplains, respectively.

The model coupling aims to replace the SWIM Muskingum routing method with the CaMa-Flood diffusive wave approach. Since the SWIM model simulates sub-basin-based runoff while the CaMa-Flood model routes the grid-based runoff, the first step is to reduce the 2D structure in the CaMa-Flood model to 1D. The original code was modified so that the CaMa-Flood runs after the sequence of sub-basin numbers. In addition, it needs to read the sub-basin-based input data, such as daily runoff generated by SWIM, river characteristics (river length, width and depth), sub-basin areas, flow directions that link sub-basins to the downstream sub-basins and floodplain elevation profiles. The core calculations for water transport in the CaMa-Flood model were not changed and the output variables remain the same. This modified version of CaMa-Flood model can be named the Sub-CaMa-Flood (sub-basin-based CaMa-Flood) model.

Here h is the flow indices for flow exceedance probabilities lower than 0.02 and H is the index of flow exceedance probability of 0.02.

A value of 1 for the NSE and FAI denotes an absolute match of predicted and measured values, while the value 0 for all percent biases (BIAS, ΔFHV and ΔPeak) means no difference between the measured and simulated values.

As a proof-of-concept study, we selected a meso-scale river basin, the Mulde River, which is one of the major tributaries of the Elbe river (Figure 1). The Mulde River basin is mainly located in Germany (ca. 98% of the total area) and the rest is located in the Czech Republic. The total drainage area at the gauge Bad Dueben is about 6,171 km². The Mulde basin is an intensive agricultural region with about 52% of the total area used as arable land. About 23% of the drainage area is forested, and 7% is grassland.

The long-term average discharge at Bad Dueben is about 64 m³/s in the period 1960–2004. The maximum daily discharge of 1,570 m³/s occurred on 14th August 2002 during the extreme 2002 Elbe flood. The flood peak discharges exceeded a return period of 100 years in most Mulde stretches (Schröter et al. 2015), causing many dyke breaches and much damage in the cities along the main river.

Four spatial maps were required to derive the sub-basin and hydrotope structure for the SWIM model: the DEM, the soil map, the land-use map, and the sub-basin map in a grid format with 100 m resolution. The DEM was provided by the NASA Shuttle Radar Topographic Mission (SRTM). The soil map of the study area was merged from the general soil map of the Federal Republic of Germany ‘BÜK 1000,’ produced by the Federal Institute for Geosciences and Natural Resources (BGR), and soil map of the Czech Republic. The land-use map was obtained from the CORINE 2000 land-cover data set of the European Environment Agency. On the basis of the DEM and the stream network, a sub-basin map of 476 sub-basins was delineated (Figure 1) with an average drainage area of 12.6 km² and an average river stretch length of 3.7 km. The purpose of delineating such fine sub-basins was to avoid long river stretches and steep slopes between sub-basins. However, we noted that some large sub-basins cannot be avoided in the delineation process (see Figure 1). The sub-basin area and the flow direction between sub-basins was also generated with the sub-basin map.

The climate data (temperature, precipitation, solar radiation and air humidity) were interpolated to the centroids of every sub-basin by the inverse distance method using data from 39 climate and precipitation stations. The discharge data at Bad Dueben, which was used to calibrate and validate the models, was obtained from the Global Runoff Data Centre in Koblenz, Germany.

The sub-basin-based river characteristics are the essential input data for the Sub-CaMa-Flood model as they determine the water storage capacity in rivers and influence the flooding water volume in floodplains. The river length was calculated based on the river network, which was generated with the sub-basin map, while the river depth and width were estimated by empirical functions of drainage areas (Allen et al. 1994). The latter two parameters were corrected by comparing the estimated and measured values at the river outlet and assuming the percent error holds for all river stretches within the basin. The river elevation is determined as the average elevation of the whole river stretch in each sub-basin. The simplified function of floodplain elevation profiles were generated by plotting the cumulative distribution function (CDF) of the elevation within a sub-basin and connecting each 10th percentile of the CDF (see more explanation in Yamazaki et al. 2011).

The flood inundation map in 2002, which is used to compare the simulation results, is a satellite-based remote sensing product provided by the National Aeronautics and Space Research Centre of the Federal Republic of Germany (DLR).

The SWIM model was firstly calibrated for the discharge at the gauge Bad Dueben in the period 1991–1997 and validated for the period 1998–2003. The selection of the validation period aims to highlight the model performance of the 2002 flood event. Note that the validation period includes the extreme flood year 2002 and the extremely dry year 2003, so the model validation can be challenging when the model is calibrated for normal hydrological years. The four evaluation criteria (NSE, BIAS, ΔFHV and ΔPeak) were calculated after each calibration run and served as the model outputs. The parameter estimation routine PEST (Doherty & Skahill 2006) was applied to minimize the difference between the criteria results and their ideal values (1 for NSE and 0 for all biases).

Table 1 lists the NSE, the BIAS, ΔFHV and ΔPeak for both the calibration and validation periods and Figure 2 shows the comparison of the observed and simulated daily and average monthly discharges using the SWIM model. The results show that SWIM can well reproduce the river discharge during the calibration period with the NSE of 0.73 and total bias of −1%. The high flows and flood peaks were also well simulated with biases of −6% and −5%, respectively. For the validation period, the low NSE value (0.53) indicates that the daily discharge was poorly reproduced compared to that in the calibration period, though the simulated average monthly discharge has a good agreement with the observed one. The total water volume and peaks were overestimated by 8% and 14%, respectively. The poor model performance in terms of NSE in the validation period can be explained by the overestimation of the extreme flood peaks and the negative discharges in the falling limb of the 2002 flood hydrograph due to the numerical problem of the Muskingum routing method.

In order to directly compare the performance between the Sub-CaMa-Flood routing model and the Muskingum routing method, we took the simulated runoff from the SWIM model presented in the previous section as the input to the Sub-CaMa-Flood model. The Manning coefficients of 0.026 and 0.11 for river beds and floodplains, respectively, gave the best results for the discharge at Bad Dueben.

Figure 3 shows the same comparison as Figure 2, but with the simulated discharge by the Sub-CaMa-Flood model. It is difficult to distinguish the differences visually between the daily/monthly discharge routed by the Muskingum method and the Sub-CaMa-Flood model in the calibration period. The Sub-CaMa-Flood model output shows slightly poorer results in term of NSE (0.70) than the Muskingum method (Table 1). The total bias remains −1% as the routing methods do not influence the total water yield. The errors in high flows and peaks are also satisfactory (−4 and −8%, respectively).

Interestingly, the Sub-CaMa-Flood model outperforms the Muskingum approach in the validation period, with the NSE of 0.77 and a peak error of 7%. The extreme flood hydrograph in 2002 was better reproduced with a reasonable peak level and a smooth falling limb. No negative discharge was generated. Such results show that the Sub-CaMa-Flood model performs more stably and reasonably for extreme floods than the Muskingum routing method in this study.

The most important feature of the Sub-CaMa-Flood model, compared to many routing methods that only route runoff to rivers, is its ability to estimate flood inundation areas. The water level and storage for both rivers and floodplains were calculated and visualized using the DEM map of the basin.

Figure 4 compares the observed and simulated flood inundation area in the Mulde basin. The red color shows the flooded area agreed by both observations and simulated results, while the brown color indicates the overestimation by the model. For the main river stretches (see the purple box in Figure 4), the simulated flood inundation area generally matches the observed flooding area with the FAI of 0.6. However, the flooded areas are overestimated along the tributaries and headwaters, showing that such large-scale routing method may be inappropriate for small-scale flood simulations or for river characteristics that are poorly estimated. The overestimation leads to a low FAI of 0.3 for the whole Mulde River basin. Since the coupled models are supposed to be applied in large river basins in the future, the present results for the main river indicate that they can be a useful tool for preliminary flood risk assessment along large rivers.

The negative discharges are a well known problem of using the Muskingum routing method (Hjelmfelt 1985). However, in this study, the negative discharge can be partly attributed to the small size of sub-basins where the water accumulation time is less than one day. In our previous study (e.g. Huang et al. 2013), the Mulde River basin was delineated into 44 sub-basins and no negative discharge was generated for this event by the SWIM model. In contrast, the Sub-CaMa-Flood model applied finer time intervals and proved to be stable for simulating such a large flood event. We should also note that the input runoff data is daily (although empirical approaches are applied in SWIM to get realistic fast runoff components), so that the Sub-CaMa-Flood model could not precisely simulate the sub-daily variation of discharge in this study.

It is always challenging to simulate large-scale flood extent at fine resolutions. Falter et al. (2016) applied the SWIM model and a 1D/2D raster-based inertial model to simulate the 2002 Elbe flood extent. The results showed an underestimation of flood extent along the main Mulde River; the underestimation can be attributed to dyke breaching which was not considered in the simulation. Alfieri et al. (2014) compared the simulated 100-year flood hazard map with the regional flood hazard map for the state of Saxony in Germany, which includes the Mulde River basin. They showed generally good agreement between the pan-European and regional maps using small-scale 2D floodplain hydraulic simulations for every 5 km along the river network. These 2D models calculating at a fine resolution are expensive in terms of computational resources. For example, Falter et al. (2016) reported that the 2D hydrodynamic simulation for the German Elbe River and its large tributaries using high performance clusters took around 19 hours for a simulation period of three months. In contrast, the Sub-CaMa-Flood routing model, which can be characterized as 1D with sub-grid representation of floodplains, took ca. 3 minutes to simulate 476 river stretches over 14 years. Taking the promising results along the main river and the fast computational time into account, the Sub-CaMa-Flood routing model can become a useful tool for preliminary flood risk assessments in large river basins, especially under an ensemble of climate scenarios.

We should also mention the weakness/limitations of the coupled models before analyzing the outputs. Firstly, as a global routing algorithm, the CaMa-Flood model itself is a simplified hydraulic model, which uses the diffusive wave equation instead of the fully 1D St Venant momentum and continuity equation. The river channel and the floodplain cross-sections are assumed to be of a simplified shape, while the floodplain elevation profile was derived from the cumulative distribution of elevations. Such simplifications mean that some topographic depressions including permanent lakes and wetlands are also considered as rivers or floodplains (Figure 4). Secondly, the estimation of river channel characteristics from the empirical functions produces large bias compared to the measurements at the outlet of the basin. It implies that the river width and depth in other sub-basins may still differ from the actual values even after a global correction for all river stretches. Since the river characteristics are of particular importance in hydraulic modeling to determine the bank-full water storage, the overestimation of flood inundation area along the tributaries in Figure 4 can be partly due to the wrong river channel information. Ideally the terrestrial survey cross-section data should be applied for each river stretch, but such data is not usually available for large-scale river basins. Thirdly, the current model version can only simulate natural flooding without considering the effects of dykes, reservoirs and other flood protection measures. It should be possible to simulate reservoir operations, water diversion and effects of polders in the future, since the SWIM model has already had a reservoir module and a wetland module in combination with the Muskingum routine (Krysanova et al. 2015). The dyke effects can also be considered by calculating flow in river channels with adjusted river depth if water volume is below the dyke-full water storage. Another approach is to consider the design floods of the dykes as an upper threshold below which inundation does not happen. However, detailed dyke information is then required along rivers. It is challenging to derive accurate flood risk maps even if dyke breaching is specifically simulated (Falter et al. 2016). Hence, we can only consider our results as potential flood risk maps.

The present study coupled a regional hydrological model (SWIM) and a large-scale hydraulic model (CaMa-Flood) to simulate river discharges and flood inundation areas in a meso-scale river basin. The Sub-CaMa-Flood routing method proved to be a good alternative to the Muskingum method for the SWIM model as it showed good performance for normal hydrological conditions as well as extreme flood events. The inundated area along the main river was also satisfactorily simulated, though an overestimation of the flooding area was found in headwaters and tributaries. The presented results show good potential to apply this coupling of models for preliminary flood risk assessments in large river basins to identify vulnerable areas, especially using large ensembles of climate change scenarios, which are required in state-of-art climate impact studies. Once the potential flooded areas are identified, smaller scale studies can be performed to produce more precise flood risk maps given the future extreme flood events. Hence, our next work is to apply the coupled models in large river basins in Europe, for example, the Elbe and the Danube, in both historical and climate scenario periods.

This work shows an example of coupling a regional hydrological model with the CaMa-Flood routing model. Since the Sub-CaMa-Flood model can be applied independently of SWIM, it is in principle possible to couple it with other regional hydrological models, which simulate sub-basin-based daily runoff. It offers the possibility not only to compare the river discharges generated by an ensemble of hydrological models, as has been done under the ISIMIP project, but also to analyze the uncertainty of estimated flood-prone areas driven by different hydrological model outputs.

The authors would like to thank Dr Yamazaki for sharing the original code of the CaMa-Flood model.