A coupled hydrological and hydrodynamic model for flood simulation

A growing number of severe floods have been witnessed over the past two decades in the UK. In the autumn of 2000, tropical cyclones caused heavy rainfall and consequently severe flooding across the UK. The successive series of floods caused two deaths and the inundation of 7,406 properties with a total economic loss of £500 million (Met Office 2014). In January 2005, Carlisle experienced the worst flood since 1822, leading to 1,844 properties being flooded, three deaths and over £400 million of economic losses. Two years later in 2007, England and Wales were hit by severe flooding, which caused 14 deaths and left thousands of people suffering from prolonged misery. More than 55,000 homes and 6,000 businesses were flooded, resulting in the highest number of search-and-rescue missions in the UK since the Second World War (Marsh & Hannaford 2007). Nationwide flood events also occurred during the course of 2012 and through the winter into 2013, which caused at least nine deaths and an economic loss of about £1 billion. With the ongoing climate changes and more intensive human activities, Evans et al. (2004) predicted that river and coastal flood risks could increase up to 20-fold and the resulting economic losses could increase to between £1.5 billion and £21 billion by the 2080s. Completely eliminating flood risks is technically and economically impossible, and so integrated strategies to manage flood risk must be in place to avoid or mitigate the adverse impacts of floods on individuals and communities (Pitt 2008). Flood modelling has provided an indispensable tool to inform the development of flood risk management strategies.

Substantial research efforts have been made to improve the numerical accuracy and computational efficiency of 2D hydrodynamic flood models. However, the existing hydrodynamic models in many cases are still computationally prohibitive for large-scale applications, especially in urban environments where high-resolution representation of complicated topographic features is necessary (e.g. Hunter et al. 2007). Hydrological models have been widely used in predicting flooding and hydrological processes in rural catchments, which require less computational time at the price of representing less detailed physical processes. Coupling hydrodynamic models with hydrological models may overcome the shortcomings of either type of modelling approach. Although hydrodynamic and hydrological models have been separately used in many practical applications, coupling of these two types of models for urban flood simulation has rarely been reported in the literature. Therefore, the aim of this work is to develop a robust and efficient tool for practitioners to perform urban flood modelling by coupling a hydrodynamic/hydraulic model with hydrological models.

There are normally three main types of method for coupling hydraulic and hydrological models: external coupling, internal coupling and full coupling (Morita & Yen 2002). The simplest and most common type is external coupling which usually employs the pre-acquired hydrographs from hydrological models as the upstream and/or lateral boundary conditions for the hydraulic models, to provide a one-way but seamless transition (e.g. Anselmo et al. 1996; Correia et al. 1998; Lastra et al. 2008; Gül et al. 2010; Bravo et al. 2012). This coupling method has also been widely applied in flood routing through complicated river network systems (e.g. Whiteaker et al. 2006; Lian et al. 2007; Bonnifait et al. 2009; Mejia & Reed 2011; Paiva et al. 2011; Kim et al. 2012; Lerat et al. 2012).

The internal coupling method has also been reported in the literature (Thompson 2004), in which the governing equations of the hydraulic models and hydrological models are solved separately, with information at the shared boundaries updated and exchanged at each or several computational time steps (Morita & Yen 2002). For example, Thompson (2004) linked MIKE-SHE hydrological and MIKE 11 hydraulic models at several prescribed points along the river. Water levels calculated from MIKE 11 can be transferred to MIKE-SHE, and overland flows calculated by MIKE-SHE are fed back to MIKE 11 at these points throughout a simulation. Few studies report full coupling, due to the complication of reformulating and simultaneously solving governing equations in a single code base.

To relax the restriction of the computationally demanding hydrodynamic models for high-resolution urban flooding modelling, this study presents a new method for coupling hydraulic/hydrodynamic models with hydrological models for application in urban catchments. The proposed numerical method falls into the external coupling category. During a simulation, hydraulic and hydrological zones are first specified according to a design flood event. The runoff or net rainfall calculated by a lumped conceptual model, known as a water balance model (WBM) (Walker & Zhang 2002), in the hydrological zones is routed with a group of pre-acquired ‘unit hydrographs’ (UH) to the hydraulic cells at the shared boundaries for higher-resolution simulations utilising a full 2D hydrodynamic model (Liang 2010) in hydraulic zones. It should be noted that although drainage system network modelling is usually an important component in urban flood modelling, it is not included in the current modelling strategy as the flood event under consideration was mainly caused by the fluvial process and overtopping and damage of flood defences.

This section introduces the coupled hydrological and hydrodynamic model as proposed, in which a catchment is required to be first divided into hydraulic zones and hydrological zones according to a design flood event.

Before an actual simulation, the full 2D shallow flow model (Liang 2010) is used to predict a design flood event in the research domain, driven by design rainfall and design inflow hydrographs at the river boundary points. To reduce computational time, the design flood simulation is run on lower-resolution grids to rapidly obtain the preliminary inundation extent. The areas inundated by the design flood are preliminarily regarded as hydraulic zones, while the rest of the domain is specified as hydrological zones. The design flood event will be chosen to be more severe than the simulated event to ensure the validity of the domain partition. In the adopted one-way external coupling approach, flow from hydraulic zones to hydrological zones is not expected. Areas less likely to be inundated, with relatively higher elevations and steeper local slopes, may also be specified as hydrological zones, which can be determined by checking the digital elevation model (DEM) and land cover maps of the research domain.

The above shallow water equations are solved using a finite volume shock-capturing Godunov-type scheme, incorporated with an HLLC approximate Riemann solver for estimating the interface fluxes. Detailed implementation and validation of the numerical scheme can be found in Liang (2010).

The external coupling method is adopted in this work, which generally employs pre-acquired hydrographs from a hydrological model as the upstream and/or lateral boundary conditions to drive the hydraulic model, providing a one-way but seamless transition. Herein, the net rainfall that leads to surface runoff in hydrological zones (the left part in Figure 2) is first calculated using the WBM hydrological model. A UH method is then used to derive the corresponding hydrographs at the hydraulic cells at the border between the hydraulic and hydrological zones to drive the high-resolution hydrodynamic simulations in hydraulic zones (the right part in Figure 2). Whilst providing a reliable approach to derive the boundary hydrographs to drive the full 2D hydrodynamic model, the UH method is easy to implement in this context. This work introduces a novel means to apply this classic method in flood modelling. As mentioned previously, flow from the hydraulic zones to hydrological zones is not considered.

It should be noted that the WBM model is not only applicable to the hydrological zones for calculating the net rainfall and lumped runoff production but also to the hydraulic zones to account for the hydrological processes in each of the hydraulic cells (e.g. the infiltration as mentioned before). In the hydraulic cells, the net rainfall resulting from the WBM model will be directly fed into the full 2D shallow flow model as an external source term.

The UH represents the hypothetical response of a catchment at the outlet to a unit input of uniform net rainfall. In this work, the concept of a UH has been extended so that it is no longer limited to providing a hydrograph at the catchment outlet, but is applicable to the catchment zone border. In order to generate a UH at the zone border, the hypothetical net rainfall (10 mm for 15 minutes) is applied to the hydrological zones, and the integrated full 2D shallow flow and WBM model is used to simulate the rainfall induced surface flow and subsequently derive the hydrographs at each of the bordering cells. These hypothetical hydrographs can be regarded as a group of generalised UHs. During real simulations, these UHs will be scaled and superimposed according to real net rainfall pattern in the hydrological zones (as shown in Figure 3) to provide accumulative hydrographs along the zone border. These accumulative hydrographs will then provide external source terms in the full 2D shallow flow model in those bordering hydraulic cells to create flooding in the hydraulic zones.

In summary, to apply the new coupled model to support flood modelling, the following steps will be taken:

• (1)

  Collect and process spatial and non-spatial data for setting up models.

• (2)

  Run the full 2D shallow flow model on a coarse grid to produce the preliminary inundation extent for a design flood event.

• (3)

  Specify hydrological zone and hydraulic zones according to the inundation extent obtained in Step 2, as well as DEM and land cover maps.

• (4)

  Use the full 2D shallow flow model to generate a group of UHs at the zone border cells, and subsequently obtain the scaled hydrographs according to ‘real’ net rainfall.

• (5)

  Run the coupled model to simulate observed/design events, driven by the scaled hydrographs and rainfall input (if necessary).

• (6)

  Process and interpret simulation results.

In this study, the newly developed coupled flood modelling tool is applied to reproduce a major flood event which occurred in Morpeth, UK in 2008. This section provides the necessary description of the case study site, the event and the data required for the simulations.

Morpeth is an ancient market town located in the Wansbeck Catchment of northeast England, UK. The Wansbeck Catchment covers a total of 331 km², and the catchment area upstream of Morpeth is approximately 287.3 km² (CEH 2018). In Morpeth, the Wansbeck River is joined by three smaller tributaries, i.e. Cotting Burn, Church Burn and Postern Burn. The main reach of the Wansbeck River has an active floodplain ranging from 100 to 300 m in width, in which Morpeth is located (Environment Agency 2005). Due to its location, Morpeth is particularly vulnerable to flooding, with 1,407 properties in the town centre identified to be at high risk of flooding (Environment Agency 2009b).

In history, Morpeth has experienced many flood events of varying levels of severity, with the most recent floods occurring in 2008 and 2012. Between 4th and 6th September 2008, northeast England was hit by an extremely heavy storm and the Morpeth weather station recorded 152.3 mm of rainfall, equivalent to 235% of the average rainfall for September (Met Office 2008a). The catchment was already saturated due to the antecedent rainfall in July and August, which caused more surface runoff and rapid rise of river levels and discharges (JBA Consulting 2009). The peak discharge reached 357 m³/s at the Mitford flow station located about 2.2 km upstream of Morpeth. The high upstream discharge combined with local overland inflow led to the overtopping and damage of flood defences for the first time since they were built about 50 years ago. The resulting floods inundated more than 950 properties in the town centre. All road networks within the town centre were largely affected. Hundreds of residents were forced to evacuate, and many of them were not able to live in their home for months. Causing over £10 million of direct damage (Coolgeology 2018), the event was estimated to have a one in 137 year return period and was recorded as the worst flood in Morpeth (Environment Agency 2009a).

The data used in this work for setting up and validating models are summarised in Table 1, including non-spatial data (observed records of rainfall, river discharge, water depth) and spatial data (digital terrain model (DTM), land cover map, soil types and flood extents). Figure 4 shows the extracted simulation domain and the location of all rain gauges in the catchment. High-quality continuous rainfall records for the current event are mainly available from four rain gauges (i.e. Wallington, Harwood Burn, Font Reservoir and Newbiggin). An inverse distance weighting method, i.e. Shepard's method (Shepard 1968), is utilised to obtain the rainfall distribution to support flood simulations at the town centre. The resulting hourly rainfall is shown in Figure 5, in which the highest daily rainfall of 72.75 mm is observed to occur on 6th September, correlating well with the recorded peak flow discharge of 357 m³/s (which also occurred on 6th September). The discharge hydrograph recorded at the Mitford gauge station at a 15-minute interval is also illustrated in Figure 5. The Mitford station (denoted as ‘A’ in Figure 4) is located at the confluence point of Upper Wansbeck River and Font River with an upstream catchment area of approximately 282 km². The recorded discharge hydrograph provides perfect upstream inflow boundary conditions for the research domain.

A 2 × 2 m high-resolution LiDAR DTM is available for the case study site. The DTM has been validated to have a vertical RMSE of about 20 cm (Parkin 2010). Surveyed river cross-section data are used to correct the DTM to give better representation of river geometry. Extra fieldwork has been carried out to measure the location and height of the flood defences to further correct the DTM. With the building outlines and locations available from OS MasterMap, the buildings are represented on the DTM by raising the ground level by 10 m. Land cover data are also obtained from OS MasterMap to define permeable and impermeable zones in the computational domain. Based on the USDA soil texture, the general Green-Ampt infiltration parameters can be decided (Rawls et al. 1982; Rawls & Brakensiek 1985), as shown in Table 2. Due to the almost saturated soil prior to the flood event, the initial moisture contents are all set to be ‘saturated’.

Finally for model validation, flood depth and extent data are made available from a post-event investigation carried out by Parkin (2010). The Morpeth Flood Action Group has reviewed the inundation extent maps and provided further clarifications (further details can be found in Parkin (2010)).

In this section, the coupled model is applied to reproduce the 2008 Morpeth flood event, and the simulation results are compared with observations and also the high-resolution numerical predictions obtained from the full 2D shallow flow model. As explained previously, the major cause of the 2008 Morpeth event was fluvial flooding from defence damage and overtopping; the effect of drainage systems is not considered in the simulation results presented here.

The simulation domain is selected to be a 19.16 km² sub-catchment where Morpeth is located. The location of the sub-catchment within the entire Wansbeck Catchment is illustrated in Figure 4 (enclosed by the red line). On the map, the Mitford flow station is located at point ‘A’, through which all of the upstream flows converge. The sub-catchments downstream of point ‘B’ have no effect on Morpeth, which provides a natural downstream boundary for the selected computational domain.

The extracted sub-catchment is divided into hydraulic and hydrological zones to facilitate the setup of the coupled model. The design rainfall over the catchment and design hydrograph at the Mitford flow station, with the return period of 200 years, are first obtained using the Depth-Duration-Frequency curves and the Revitalised Flood Hydrograph model for the Wansbeck Catchment, available in the Flood Estimation Handbook (Kjeldsen 2007). Subsequently, a design flood event is then generated by the full 2D shallow flow model based on the design rainfall hyetograph and design hydrograph. The resulting flood extent can be then employed as a reference to divide the extracted sub-catchment into hydraulic and hydrological zones.

When classifying the hydraulic and hydrological zones, only rough inundation extent is required and so simulations can be performed on a coarse grid at 5 m resolution to reduce computational cost. The resulting inundation extent is displayed in Figure 6. The fluvial inundation area can be preliminarily regarded as a hydraulic zone, and local adjustments may be made according to the land cover. For example, in Figure 6, the inundation area around point ‘A’ may be caused by local rainfall, but it is still included in the hydraulic zone because it is mostly covered by buildings and is close to the fluvial inundation area. The final hydraulic zone is outlined in Figure 6 and the remaining areas are regarded as the hydrological zone. Since the return period of the September 2008 flood event is evaluated as 137 years, the catchment division scheme from the one in 200 year design flood event is considered suitable for the assigned simulations.

The extracted sub-catchment/simulation domain is discretised at 2 m resolution to represent the complex urban topography. The numbers of different types of computational cells in the domain are listed in Table 3. In the simulation using the full 2D hydrodynamic model, a total of 4,789,781 cells are used to solve the shallow water equations to predict rainfall-induced flooding process. In the coupled model simulation, only 269,480 cells in the hydraulic zone are handled by the full 2D model but the hydrological zone covered by 4,520,301 cells is treated as a lumped area using hydrological models.

For all of the simulations, a steady flow in the river channel, driven by a steady low inflow of 10 m³/s at the inlet, is obtained as the initial condition for the full 2D shallow flow model. The Manning coefficient, n = 0.02 m^–1/3 s, and the acceleration due to gravity, g = 9.81 m/s², are used.

Before reproducing the 2008 Morpeth flood event using the coupled model, the full 2D shallow flow model is employed to create the UHs at the hydraulic and hydrological zone border cells under a uniform net rainfall of 10 mm of 15-minute duration in the model domain. Figure 7 displays a UH in one of the border cells as produced. The UH is then scaled and superimposed according to the real net rainfall pattern to generate an accumulative hydrograph from the hydrological zone at that particular cell, as illustrated in Figure 8. The accumulative hydrographs in the border cells will then be used to drive the simulation in the hydraulic zone.

The 2008 Morpeth flood is reproduced for 50 hours starting at 4 pm on 5th September on a desktop PC with an Intel^(R) Core^(TM) i5 CPU. Six critical points located at different parts of the town centre are selected to record the time histories of water depths to compare with measurements. Details of these critical points are documented in Table 4 and Figure 9.

In Figures 10–13, the inundation maps at 11 am, 1 pm, 3 pm and 5 pm on 6th September obtained from the coupled model are compared with those obtained using the full 2D shallow flow model, as well as the post-event investigation. The full 2D shallow flow model has produced similar inundation extents in comparison with the coupled model in the hydraulic zone at 11 am and 1 pm. However, the inundation extent at 3 pm from the coupled model is larger at the left bank in Central Morpeth and at the right bank in High Stanners. At 5 pm, this difference becomes marginal again. This may be explained by the coupling method used in the coupled model. The adopted UH approach may route the flow from the hydrological zone to the hydraulic zone faster than using the full 2D shallow flow model, and therefore the flow from the hydrological zone predicted by the coupled model may join the fluvial flood in the hydraulic zone earlier to produce a larger inundation extent at around 3 pm. However, this effect is less pronounced when the rainfall reduces, and so the two inundation extents become similar again at 5 pm, after the intense rainfall period terminates.

In comparison with the post-event investigated inundation maps, the two simulated inundation extents, respectively from the coupled model and the full 2D shallow flow model, are observed to be larger at 11 am and 1 pm. A closer match has been produced for flood extents at 3 pm. At 5 pm the simulated inundation extents from both models achieve the closest overall match with the post investigation survey. The cause of this discrepancy may be due to limitations in the measured data. The Environment Agency debris data and the crowd-sourced information may not represent the full flood propagation process, and this means that the maximum inundation extent (at 5 pm) is likely to be much more reliable than the records at other times.

Water depths at the six critical points are recorded by both models and plotted in Figure 14. At all of the critical points the water depths generally rise until reaching a peak, and then fall. After this, the water depths increase again to a second but lower peak before a final retreat of the flood water. The ‘twin peak’ shape of the water depth series reflects the similar shape of the inflow hydrograph, as shown in Figure 5. The water depths from the coupled model match reasonably well with those obtained by the full 2D shallow flow model in all of the critical points, particularly at P1, P4 and P6. However, deviations are seen in P2, P3 and P5, where the coupled model generally produces slightly larger depths than the full 2D shallow flow model. As a matter of fact, points P2, P3 and P5 are located relatively nearer to the zone border than the other critical points, and are therefore more influenced by the lateral inflows from the hydrological zone. The faster routing method of UH may explain the larger depth as well as the earlier arrival time predicted by the coupled model. In practice, the use of slightly larger hydraulic zones is therefore recommended and may effectively minimise the effect.

When compared with the measured water depths in the post-event investigation, obvious discrepancy can be found. Results from both models give shallower depth than the measured data at P1 and P2 but deeper depth at P6. The full 2D shallow flow model performs better at P3 and P5, while the coupled model behaves better at P4. Negating elements of the draining process along with limitations of the crowd-sourced observation data may explain these differences. As a whole, the results from the two models are considered to be consistent and the coupled model is preferable for wider applications due to its high computational efficiency while maintaining satisfactory simulation accuracy.

In terms of model performance, the coupled model costs 175.2 hours of CPU time to complete the simulation, while the full 2D model takes 2,137.4 hours. For this particular simulation, the coupled model is therefore 12.2 times computationally more efficient, i.e. the computational efficiency has been dramatically improved by the coupled model, against the full 2D shallow flow model for the same simulation, without significantly compromising the simulation accuracy. This improved simulation efficiency is consistent with the use of far fewer computational cells for hydrodynamic calculation in the coupled model. As shown in Table 3, the full 2D simulation involves about 17 times more hydrodynamic computational cells than the coupled model. The model speedup is hence related to the size of hydrological zones and larger hydrological zones will clearly lead to more efficient simulations.

This paper presents a new coupled hydrodynamic and hydrological model for urban flood simulation including flow from upstream steep catchment areas. The model has the following features:

• The computational domain/catchment is divided into hydraulic and hydrological zones according to a design flood event to ensure flows only occur from hydrological zones to hydraulic zones and hence facilitate an effective one-way coupling method.

• A lumped conceptual hydrological model is used to improve computational efficiency as well as maintain reasonable representation of hydrological processes.

• Overland flows into the hydraulic zone are estimated and distributed using a series of UHs generated by pre-running a hydrodynamic model.

• Inside the hydraulic zone, hydrological processes caused by precipitation and soil infiltration are approximated using a lumped conceptual hydrological model and taken into account in the full 2D hydrodynamic model as external source terms.

The new coupled model has been successfully applied to reproduce a major flood event occurring in Morpeth, northeast England in September 2008. The results are consistent and compare reasonably well with the numerical predictions obtained using a full 2D hydrodynamic model and also the post-event investigation data. For the simulations conducted in this work, the coupled model is over 10 times computationally more efficient than the full 2D model when conducting the same simulation with similar prediction accuracy. The current coupled model therefore provides a more efficient simulation tool for large-scale high-resolution flood simulation.

The study is partly supported by the National Key R&D Program of China (2017YFC1502703, 2017YFC1502706). Particular acknowledgement also goes to Dr Geoff Parkin at the Newcastle University, UK who has provided access to observed data and shared useful information about the case study site.