Integration of hillslope hydrology and 2D hydraulic modelling for natural flood management

There has been a growing number of investigations attempting to model the effectiveness of NFM measures for flood risk reduction (Nicholson et al. 2012; Ghimire et al. 2015; Metcalfe et al. 2017, 2018), although relatively few have focussed on modelling strategy (Hankin et al. 2017a). Such measures comprise distributed changes in the landscape aiming to emulate natural processes to store, slow or infiltrate water to help reduce flooding and have been collectively referred to as working with natural processes (Pitt 2008; Burgess-Gamble et al. 2017), natural flood management (UK and used hereinafter), Natural Water Retention Measures (Europe) and Nature-Based Approaches (international, e.g. World Wildlife Fund 2016). Such measures include in-stream leaky barriers (Metcalfe et al. 2017), engineered log jams (Addy & Wilkinson 2016), riparian or widespread tree-planting (Carrick et al. 2018), runoff attenuation features (RAFs) in the landscape (Nicholson et al. 2012; Metcalfe et al. 2018), peatland restoration to increase roughness (Holden et al. 2008; Shuttleworth et al. 2019), and river restoration including re-meandering and floodplain reconnection to hold back and temporarily store more water which has been popular across Europe (e.g. Office International de l'Eau 2013). They present a challenge to hydrologists and hydraulic modellers to first identify the catchment processes that might be altered, and then to justify changes or ‘shifts’ to effective parameter values controlling these processes in their models.

NFM has often been seen as a low-cost alternative to traditional flood risk management techniques and is being deployed by communities at risk unable to secure funding for measures such as hard flood defences. Many such examples are included within 65 UK case studies recently compiled in Burgess-Gamble et al. (2017), with online mapping (naturalprocesses.jbahosting.com/) and a dedicated European website (http://nwrm.eu/). The recent UK Evidence Directory (Burgess-Gamble et al. 2017) and review papers (Dadson et al. 2017; Lane 2017) summarise the current state of knowledge, pointing to evidence at the small (10 km²) to medium (100 km²) scales and for small (<10% annual exceedance probability – AEP) to medium floods (1% < AEP < 10%), but indicate little direct evidence for effectiveness at larger scales (>100 km²) and major flood events (<1% AEP). European research (Office International de l'Eau 2013) and other international advice on nature-based methods (World Wildlife Fund 2016) provide much advice on incorporating risk reduction measures based on green engineering.

Strategies for distributed modelling of NFM (e.g. Hankin et al. 2017a) and testing the resilience against realistic spatial extremes (Hankin et al. 2017b) and synchronisation effects (Pattison et al. 2014) rely on modellers first being able to represent distributed hydrological responses accurately, and then represent the changes incurred to those processes as a result of small-scale distributed nature-based interventions realistically. Since the effective parameterisation of these processes and changes to processes are very uncertain, a framework for uncertainty analysis such as the generalised likelihood uncertainty estimation methodology (Beven & Binley 1992; Beven & Binley 2014) becomes a necessity, and methods for detecting change in ensemble modelling results with and without nature-based measures important (Hankin et al. 2017a).

Physics-based models can offer a useful lens for upscaling small-scale distributed processes such as increases in friction, storage and infiltration to study their combined effect at the catchment scale, yet also reveal a blurry image of our knowledge of exactly how to represent often variable and difficult-to-measure processes. The crux for modellers has therefore been how to represent very distributed and evidenced small-scale changes to hydraulic and hydrological processes, using efficient tools to capture the main interactions, whilst allowing for ensemble simulations and uncertainty estimation. Often a single type of model is used to represent the whole system, whereas using a combination of models (which are becoming more versatile to integrate) can improve the overall accuracy of process representation. The ‘evidenced’ changes to effective distributed parameters require further research, currently the focus of three ongoing UK Natural Environmental Research Council research projects (https://nerc.ukri.org/research/funded/programmes/nfm/).

This paper demonstrates a new integrated modelling framework, isolating one large runoff attenuation feature in a relatively small catchment, to test the effectiveness of a new modelling strategy to assess changes in catchment response whilst also quantifying uncertainty. The paper aims to assert the need for such a framework with the following three objectives:

• Bring together two fast technologies that improve integrated modelling of the distributed hydrology and hydrodynamics required to simulate small-scale NFM interventions and scale the effects up to the catchment scale, whilst allowing for ensemble simulations to explore uncertainties. Hillslope processes are modelled using the R version of Dynamic Topmodel (Metcalfe et al. 2016) with channel and floodplain hydrodynamics modelled using JFlow (Lamb et al. 2009).

• Undertake cascade modelling of an ensemble of ‘acceptable’ model structures (Beven 2006; Blazkova & Beven 2009) and define the modelling uncertainties in our estimates of the effectiveness of NFM. Many investigations are using modelling to scale up the effects of distributed catchment measures, to understand the change in response at the outlet, but few are attempting to quantify uncertainties in their answers.

• Use more of the 2D modelling outputs to drive a better understanding of processes having a strong influence on flooding and flood risk reduction – here we use 2D depth and velocity data to estimate shear stresses and generate erodibility maps.

Swindale is a small, 15 km² upland valley near Shap in Cumbria, UK, with water draining from upland peatland cascading down to a valley with runoff characteristics dominated by impermeable glacial till. It has a high annual average rainfall (standard period annual average rainfall (SAAR) = 2,445 mm) and low baseflow index (BFIHOST = 0.29). Flows are gauged upstream of the United Utilities (UU) abstraction point (NY5146113169) that feeds water to Haweswater in the adjacent valley. The UU flow gauge (station 761113) has telemetry at 15 min 1997–2015; 2016 to present, with a gap whilst the intake and fish pass were upgraded. This unfortunately results in no data for the extreme event of Storm Desmond in winter 2015, so a large storm November 2009 is studied here with an estimated annual exceedance probability of between 50% and 20% of the scale where there is more evidence for NFM interventions being effective (Burgess-Gamble et al. 2017).

The Swindale catchment has been the subject of recent investigations by UU and the Royal Society for the Protection of Birds (RSPB) into land management and peatland restoration strategies, river restoration and riparian and upland tree-planting. The effects of these distributed NFM measures were previously modelled (RSPB 2017), identifying that the most effective single measure was a large RAF created by restricting flow out of a large storage area at Dodd Bottom (marked as blue in Figure 1). For the purposes of modelling at a fine temporal resolution, the closest 15-minute Tipping Bucket Rainfall gauge in the vicinity was Deeming Moss (NGR NY5545206553, SAAR = 1,982 mm), which was proportionately scaled using the annual average rainfall based on catchment descriptors, although there will be unknown local orographic and rain shadow effects. A digital terrain model (DTM) with detailed resolution (0.5 m) was available from the RSPB (collected by Glasgow University), in the vicinity of Dodd Bottom and the restored river meandering. This was resampled to 2 m resolution using bilinear interpolation, and feathered into coarser scale topography (2 m composite DTM based on LiDAR and where unavailable, Synthetic Aperture Radar) to increase the computational efficiency of the model calculations.

We demonstrate a cascade modelling approach that couples Dynamic Topmodel (Beven & Freer 2001; Metcalfe et al. 2016) with JFlow (Lamb et al. 2009). It takes advantage of the computational efficiency and parsimonious parameterisation of hillslope hydrology in Dynamic Topmodel, coupled with the fast 2D hydrodynamic modelling offered by the graphics-card implementation of the Shallow Water Equations (SWE) in JFlow (for benchmarking study, see Hunter et al. 2008 and the updated speed testing for the SWE solver in Lamb et al. 2009). This highly efficient combination allows us to undertake ensemble simulations to explore uncertainties, whilst capturing essential surface/subsurface interactions in the wider catchment, and the complex channel–floodplain hydrodynamic interactions needed to understand the influence of in-channel NFM measures. Other rainfall runoff models could be used in a similar approach, but here Dynamic Topmodel has been used based on it having a low number of parameters, whilst being able to represent the key processes in a set of distributed hydrological response units (HRUs) which dynamically interact with one another. Many rainfall runoff models such as PDMs, hydrological predictions for the environment (HYPE) and SWAT similarly use efficient HRUs to generate responses reflective of local topography or soils, but few use a flow distribution matrix to pass flows between the landscape units to more properly reflect celerities in the saturated zone responses. In the R version of Dynamic Topmodel (Metcalfe et al. 2016), saturation excess overland flow can be generated in an upslope HRU and be passed downslope to an HRU with a deficit and re-infiltrate. Here, we require a general NFM modelling framework that can represent such processes since these are processes that we seek to influence through, for instance, tree-planting or reduced grazing intensity. In the discussion we also demonstrate a similar coupling between Dynamic Topmodel and HEC-RAS2D, giving similar results and highlighting the increasing versatility of 2D hydrodynamic packages.

The model cascade is a natural extension of the modelling in the Brompton catchment (Metcalfe et al. 2017), in which the distributed runoff from Dynamic Topmodel was fed into a network solution of the 1D St Venant equations to explore leaky barriers, that were designed largely to store water in the channel. Here, we require that the general modelling framework is capable of representing the effects of NFM interventions that promote additional flood flow attenuation by encouraging water onto the floodplain where it was previously poorly connected. The method also builds upon the direct rainfall with losses 2D modelling approach (e.g. Hankin et al. 2017a) that use Revitalised Flood Hydrograph losses model (Kjeldsen et al. 2005) or approaches that simulate distributed infiltration using an effective hydraulic conductivity. None of these approaches are capable of modelling surface–subsurface interactions, such as saturation excess runoff flowing overland until it reaches an unsaturated area, whereupon it re-infiltrates (run-on). Such interactions can, however, be represented in Dynamic Topmodel, so here we combine this more realistic hillslope process representation with the 2D hydrodynamic routing to capture the complexity of flows over steep and diverse terrain and distributed interactions between channel and floodplain, where many NFM measures are currently being experimented with (for examples, see Office International de l'Eau 2013; Clilverd et al. 2013). The cascade of models, therefore, provides a useful compromise between detail in process representations and computational efficiency and a framework for exploring NFM and uncertainties.

Before Dynamic Topmodel is used to simulate a catchment, some preprocessing of the Digital Terrain Model (Figure 1) is required in order to divide the catchment into areas that respond in a hydrologically similar way. These are the HRUs (shown in Figure 2) in Dynamic Topmodel and are typically based on accumulated upslope area, slope or distance to a channel. A flow connectivity matrix provides the relationship between the direct runoff between these HRUs and into the channel, represented initially through a time delay histogram in order to route the flows to the outlet, and permitting an assessment of performance at the UU flow gauge. Driven by rainfall and evapotranspiration for the winter storm modelled (November 2009), the different parameters controlling the runoff response in the HRUs were sampled from uniform random distributions between typical ranges expected for this upland catchment. Swindale has peat in the uplands and slowly permeable, shallow soils lower down dominated by glacial till deposits synonymous with slowly permeable soils. More details of the parameters in the Dynamic Topmodel model are provided in Beven & Freer (2001), but the most sensitive parameters include the saturated surface transmissivity, the exponential transmissivity shape parameter, m, the route zone storage coefficient and the vertical time delay parameter, Td.

Solutions meeting a relatively high Nash–Sutcliffe efficiency (NSE) (>0.85) were classed as behavioural. The more stringent limits of acceptability approach (Beven 2006; Blazkova & Beven 2009; Liu et al. 2009; Hollaway et al. 2018) was not used at this stage, but the intention is to use this for new applications to larger catchments. Insisting on a relatively high performance (NSE > 0.85) resulted in 1,482 behavioural simulations out of the ten thousand Monte-Carlo simulations that were undertaken, and these were taken forward and fed into the 2 m resolution JFlow model. This was achieved by taking the simulated distributed runoff (subsurface and surface flows) and feeding the time series into along-reach internal flow boundaries in the 2 m resolution JFlow mesh, represented in Figure 2. Here, the colour shading represents the different HRUs where each produce similar surface and subsurface runoff responses based on the model of Metcalfe et al. (2016). In this investigation, eight in-stream reaches were used, represented as HRUs that receive flows via the flow distribution matrix in Dynamic Topmodel. These time series are distributed over internal flow boundaries as through-time source terms in the SWE in JFlow.

In the discussion, we also compare the same approach to model integration using a combination of Dynamic Topmodel and HEC-RAS 2D (using a diffusion wave), using a much coarser resolution numeric grid (10 m). We discuss how this scheme is able to incorporate detailed subgrid topographic data (2 m) despite the coarse grid (10 m), and explore additional outputs from the 2D hydrodynamic models, such as generating spatial distributions of shear stresses and erodibility maps. These outputs can help understand sediment processes that are also important to understanding risk reduction and the long-term performance of RAF-based NFM measures, which may, if designed incorrectly fill up with sediment.

Figure 3 shows the initial spatiotemporal results of a single coupled simulation, with and without an in-channel obstruction downstream of Dodd Bottom (Figure 1) to throttle flows using a reduced outflow width of 4 m so as to store more peak flows temporarily in the areas known as Dodd Bottom labelled in the lower left. We focus effort here on this single large RAF and model it in isolation from other measures to test the framework. The re-meandering of the watercourse in the valley bottom undertaken to improve biodiversity is included in both scenarios and can be seen in darker blue (deeper) water in Figure 3 where it is evident that the floodplain in this area is inundated for the November 2009 event. Critically, there is no evidence of unrealistic stepping in the predicted depths near the boundaries between the river reaches receiving different time series from Dynamic Topmodel, apart from close to the intervention, where the red shaded depth grid shows the additional storage due to constricting the width of the outlet to the RAF to 4 m. Use of Dynamic Topmodel alone would not have been able to accurately represent and show this level of important detail in the floodplain. The inset hydrograph shows a slight underprediction compared with the observed 48.3 m³ s⁻¹ peak flow at the UU river gauging station, although the ensemble results are analysed further below, showing a range of under- and overpredictions of the peak.

A sample of 120 of the 1,482 ensemble baseline behavioural simulations is shown in Figure 4 highlighting a range of peak flows in comparison with the gauged flow. The overall performance is strong (NSE ∼ 0.8–0.9), although there is a delay in the early part of the rising limb for all the model runs and a slight difference in response near the peak with observations suggesting a broader peak.

Rather than plot the time series for the intervention scenario, it is more informative to look at the difference between baseline and intervention (following Metcalfe et al. 2018). Figure 5 shows the ‘difference hydrograph’ between the baseline and the RAF intervention scenario for the sample of 120 simulations, for which positive values represent a reduction on the rising limb of the hydrograph, and negative values indicate successfully attenuated flows returning back to the watercourse. The maximum difference in the peak is 1.8 m³ s⁻¹ ± 1.4 m³ s⁻¹, although the largest changes tend to occur in advance of the peak in the hydrograph that occurs around hour 20 at the UU flow gauge. This confirms what was suspected from observations of the basic geometry of the watercourse (see Figure 1), in that Dodd Bottom is on a faster-rising part of the catchment, and the additional attenuation tends to result in taking volumes from the rising limb of the hydrograph as considered at the UU flow gauge.

The reduction in the peak flow afforded by the additional temporary storage is not so pronounced by the time the flows reach the UU flow gauge (Figure 5). This indicates that the restriction applied (4 m) needs to be narrower, and that restricting flows in this area from what is a faster rising part of the catchment is not an effective strategy and to improve it further, flows from the main channel may also need to be diverted into Dodd Bottom and multiple interventions may be needed. However, this is a sensitive site, and building bunds, structures and diversions across unique moraine formations is not likely to be acceptable, with the added constraint that the volume retained if impounded would come under the UK Reservoirs Act and thereby require regulation.

In addition to the spaghetti time series plots of the ensemble predictions in Figure 5, it is important to view the ensemble results spatially and develop methods for constraining model outputs which can also be helpful for expressing additional limits on the acceptability of model structures. An example of this can be found in Aronica et al. (1998) where the approximate timing of part of the flood wave exiting one area of the Imera flood basin was used as a fuzzy acceptability criteria, which if not simulated, was used to effectively reject models that might otherwise have generated an acceptable behaviour.

To understand the differences between Monte-Carlo simulations, it is once again more informative to look at the differences in predicted depths in this catchment, since the high performing simulations (NSE > 0.85) have similar extents. Figure 6 shows the performance-weighted depth grid in the central pane for the NFM scenario (6e—on scale of 0.0–3.0 m), with the eight outside panes showing the difference between this most likely depth grid and eight other NFM Monte-Carlo simulations (panes 6a–6i excluding 6e on a scale of −0.1 to +0.1 m). Some of these show positive differences (6a, 6 h, 6i), negative differences (6c, 6d, 6 g), and some not much difference (6b, 6f), but the range of predicted depths is expected based on the range of differences in Figure 5. This example illustrates how there can be differences in spatial behaviour, including potentially flood pathways within our acceptable simulations, which may coincide with knowledge of flooding that can be used to constrain the model behaviour further and reject certain modes.

The key advantages of the new approach demonstrated here include efficient, distributed representation of hillslope processes using Dynamic Topmodel, combined with more accurate 2D hydraulics of flows over a complex terrain. The previous coupling of Dynamic Topmodel with 1D hydraulic routing has proven successful in a similar scale UK catchment (Metcalfe et al. 2017), including representation of the hydraulics of in-stream leaky barriers. Such measures are also being designed to spill flows onto the floodplain to enhance flood attenuation, and so a 2D model becomes more desirable. The use of the fast 2D SWE model, JFlow, allows the impacts of in-channel modifications such as ‘leaky barriers’ and floodplain reconnection, and failures of these features, to be more adequately explored at larger scales. This includes assessing their positive and negative effects on sediment transport, in ways that take into account the large uncertainties inherent in the modelling process (including the parameters representing change within the models).

To show that the coupling Dynamic Topmodel with JFlow is a generic approach, we present a preliminary simulation of Dynamic Topmodel with an alternative 2D hydraulic model, namely, the licence-free HEC-RAS2D (https://www.hec.usace.army.mil/software/hec-ras/downloads.aspx). The HEC-RAS2D model has a new type of numerical scheme that permits efficiencies, whereby a relatively coarse numerical grid can be used (e.g. 10 m), for which the subgrid topography (2 m) is taken into account. This is achieved through precomputing hydraulic tables for each coarse numerical grid cell, comprising the conveyance across each face as a function of subgrid geometry, and then also a table of storage volume in the cell as a function of water elevation. This permits the modelling of channels that pass through the numerical cells if they are represented in the subgrid topography. Figure 7 shows the results of using HEC-RAS2D with a 10 m numerical grid based on a 2 m resolution subgrid topography compared with JFlow. The solution shown is driven by 8 time series distributed along 8 inflow reaches that are identified as HRUs in the Dynamic Topmodel modelling scheme. The red-dashed HEC-RAS 2D output is compared with the gauged data (black dash) and one of the JFlow ensemble runs. The slight differences between JFlow and HEC-RAS2d are likely to stem from the different resolutions: JFlow uses a 2 m grid, HEC-RAS 2d uses a 10 m grid with hydraulic tables based on the 2 m sub-grid topography; and the equations: HEC-RAS 2D uses an implicit solution to the Diffusion Wave, whilst JFlow solves the full SWE using an explicit scheme.

The approach yields very similar behaviour and performance at the gauge, and this increases confidence in the outputs from both models that either could be used in this model cascade approach.

To further illustrate the advantages of using 2d hydraulic models to better represent channel/floodplain flows resulting from the hillslope outputs Dynamic Topmodel, we very briefly show the ease to which the hydraulic modelling could be extended to understand the spatial distribution of shear stresses and sediment erodibility, critical to understanding the performance of NFM measures longer term. A problem with RAF measures is the potential for reduction of flow velocities and subsequent sediment deposition within the temporary storage area, or localised increases in velocity at the outlet to the storage area, creating scour and unwanted bed erosion, which can both lead to a loss of performance in terms of flood risk reduction.

The study demonstrates the feasibility of utilising the best aspects of two models, using Dynamic Topmodel for hillslope hydrology and JFlow or HEC-RAS2D for 2D channel–floodplain hydraulics, in order to simulate the detail of the processes that need to be represented with uncertainty in an NFM scheme evaluation. The credibility of the simulations of subsurface flows beneath hillslopes is increased with the use of Dynamic Topmodel for this component, while the credibility of the 2D floodplain hydraulics increased with the use of JFlow (or alternative 2D hydraulic models such as HEC-RAS2D). We have demonstrated that the approach is capable of a detailed simulation of interventions applied to floodplains (e.g., storage impoundments or channel–floodplain reconnection).

The approach has successfully advanced the coupled distributed hydrology and 1D hydraulics approach taken by Metcalfe et al. (2017) to 2D, permitting investigation of the upscaled impact of very small-scale (2 m²) NFM measures at the catchment scale (15 km²), whilst accounting for dynamic interactions in the drainage network. We have made use of the increasing versatility of different 2D hydrodynamic software packages to incorporate lateral inflows and hydraulic structures directly into a 2D mesh and allowing flexible mesh around key floodplain features such as embankments. Whilst there are numerous examples of broadscale coupling of hydrology and hydraulics models (Nguyen et al. 2015; Cea & Rodriguez 2016), these recent developments make integration of distributed hillslope-hydrological processes with distributed 2D hydrodynamic modelling considerably more feasible for a range of catchment investigations.

We have demonstrated that the approach can be used for identifying the quantum of change due to an isolated runoff attenuation feature, and also the uncertainty associated with this, given as: 1.8 m³ s⁻¹ ± 1.4 m³ s⁻¹. This represents a 4% ± 2% change at the flow gauge, which can help decision makers understand whether the intervention is worth investing in or re-engineering as a risk reduction strategy. In reality, further scenarios can be investigated, some diverting more of the flow from the main channel into Dodd Bottom, or through adding multiple features.

The approach will next be scaled up to the 210 km² Kent catchment, with 5 flow gauges and 72 km of watercourse as part of the NERC-funded Q-NFM project (https://www.lancaster.ac.uk/lec/sites/qnfm), to investigate the wider applicability of the approach. This larger catchment poses a computational problem, whereby it will no longer be feasible to simulate all of the behavioural Dynamic Topmodel scenarios using 2D hydraulic modelling at 2 m resolution, and a further sampling strategy will need to be devised. Since many of the hydraulic outputs are similar, we intend to record statistics of each of the inflow time series from Dynamic Topmodel for each reach (such as peak, volume, time of peak) and use these to distinguish between scenarios likely to result in different spatial behaviours. These can then be further constrained with local knowledge of particular flooding pathways during recent flooding.

This project has been supported by NERC grant NE/R004722/1 and has been a collaborative effort between JBA Consulting and Lancaster University in the UK. The Environment Agency is thanked for the use of rainfall and terrain data under licence CL77737MG, while Lee Schofield of the RSPB is thanked for the use of the additional terrain data and streamflow. Thanks also to members of the Q-NFM team at Lancaster, Trevor Page, Ann Kretzschmar and Rob Lamb for their helpful discussions. The ensemble Dynamic Topmodel simulations were undertaken by Peter Metcalfe, who sadly is no longer with us to see the results.