Application of data-based mechanistic modelling for flood forecasting at multiple locations in the Eden catchment in the National Flood Forecasting System ( England and Wales )

The Delft Flood Early Warning System provides a versatile framework for real-time flood forecasting. The UK Environment Agency has adopted the Delft framework to deliver its National Flood Forecasting System. The Delft system incorporates new flood forecasting models very easily using an “open shell” framework. This paper describes how we added the data-based mechanistic modelling approach to the model inventory and presents a case study for the Eden catchment (Cumbria, UK).


Introduction
Between 2004 and 2008 a component of the Flood Risk Management Research Consortium (FRMRC) investigated databased mechanistic (DBM) real-time flood forecasting and forecast uncertainty (Young et al., 2006).FRMRC phase 2 investigated uncertainty in inundation modelling using the Eden catchment as a test site.During this time, the UK Environment Agency launched the Risk-Based Probabilistic Fluvial Flood Forecasting study (Sene et al., 2009).The two projects exchanged ideas leading to the present authors developing a DBM module for the Delft Flood Early Warning System (Delft-FEWS).The module takes advantage of the Delft-FEWS open shell framework, which provides tried-and-tested code to drive models with real-time hydrology/meteorology data and permit data assimilation.The end user is also given access to a sophisticated interface for visualizing and interacting with forecasts and data.
This paper provides a description of the DBM real-time flood forecasting method (Appendix A), a section describing the design and implementation of the FEWS module, and a case study building a semi-distributed DBM representation of the Eden catchment with forecast locations at Carlisle, which was extensively flooded in 2005(Mayes et al., 2006;;Neal et al., 2012), together with intermediate forecast points at Great Corby, Greenholme, Linstock, Cummersdale and Harraby Green.

The Eden catchment case study site
The Eden catchment (2280 km 2 ) is located in Cumbria, North West England.The main channel is approximately 145 km long and discharges into the Irish Sea at the Solway Firth.Average discharge at Sheepmount (Carlisle) is 52 m 3 s −1 ; maximum recorded discharge occurred during the January 2005 flood event and was estimated at 1500 m 3 s −1 .The principle sub-catchments are the Upper and Lower Eden, Eamont, Irthing, Petteril and Caldew; these are shown in Fig. 1.The central Lake District to the west, the Pennine moors to the east, and parts of the Kielder Forest region to the northeast generate significant runoff within the catchment.The Upper Eden forms the largest sub-catchment at ∼ 600 km 2 , with a number of tributaries originating from Mallerstang to the west and the Pennines to the east; both groups containing peaks over 700 m.The Lower Eden catchment contains a number of Pennine peaks on its eastern side, including Cross Fell (893 m).The western section is mainly improved farm land.The Eamont sub-catchment (158 km 2 ) carries runoff from the central Lake District (including Helvellyn 950 m), where total rainfall is the highest in the catchment (approximately 1700 mm yr −1 ).Two lakes regulate flow along the Eamont: Ullswater (884 ha) and Haweswater (387 ha).The Irthing sub-catchment (335 km 2 ) is important from a flood management perspective as it is the only river carrying runoff from the northeast, close to the headwaters of the Tyne.This area was a major runoff generating region during the flooding in January 2005, contributing flows of 250 m 3 s −1 .The Caldew sub-catchment (244 km 2 ) is important for flooding at Carlisle as it drains from the high altitude regions of Skiddaw (931 m), which is a region of steep impermeable igneous rock resulting in a fast catchment response.The Petteril sub-catchment (160 km 2 ) is mainly improved lowland farms.The population within the catchment is approximately 240 000; the principal population centres are Carlisle, Penrith and Appleby.
Forecast lead times in the catchment are limited to a maximum of seven hours by steep relatively impermeable uplands, areas of glacial moraines, high rainfall, and steep channel geometry.Numerical weather predictions can extend this horizon but introduce further uncertainty.The distribution of rain-generating regions increases the forecasting challenge: Inputs from the central Lake District, Pennines, Skiddaw, and Kielder convolve to a complex signal observed at Carlisle.The UK Environment Agency closely monitors the catchment with a telemetered network of 31 level and 16 rainfall gauges generating a data field at 15 min intervals.Figure 1 shows the Eden catchment, sub-catchments and gauges used by the DBM model described in this paper.

The DBM flood forecasting model
The FEWS DBM case study uses a network of 6 model nodes broadly representing the Eden sub-catchments.The configuration of the nodes is shown in Fig. 2 together with the Environment Agency gauge sites supplying input (rain or level) and output (level) to each node.Appendix A gives a description of the DBM real-time flood forecasting method.The individual transfer functions, input non-linearity functions and Kalman filter hyperparameters of the Eden DBM model were identified, estimated and optimised using the RIV, RIVID and SDP functions from the Captain ™ toolbox (Taylor et al., 2007) and the Matlab ™ lsqnonlin function.The semi-distributed representation adds complexity but has the attractive properties of (1) accounting for spatial variability in the rainfall field, (2) providing forecasts at intermediate sites, and (3) introducing robustness should one or more data streams fail.The DBM approach models the relationship between an upstream level (or rainfall) as input and the downstream level.The resulting forecasting model is not, therefore, subject to a direct mass balance constraint, but this can be an advantage in times of flood when the rainfall monitoring network will not necessarily give an accurate input signal.It also means that a rating curve does not introduce additional uncertainty -particularly significant during extreme floods.At Carlisle, the peak water level in 2005 was 1.66 m above previous discharge measurements, this resulted in significant changes to the original rating curve estimates (Spencer et al., 2006).In addition, it is often forecasts of water level that are actually required for flood warning purposes.A number of previous studies have shown good results by forecasting level rather than discharge (see for example Romanowicz et al., 2006;Leedal et al., 2009;Alfieri et al., 2011).

Forecast uncertainty
The Kalman filter algorithm described in Appendix A operates as an estimator for the probability distribution of the model states.Online data assimilation updates the state, error covariance, and forecast estimates.Equation (A6) uses this information to estimate a variance for forecast error.The statistical assumptions of the Kalman filter and the form of the heteroskedastic variance model influence the accuracy of uncertainty estimates.Each is a simplification of the real system, and therefore calibration and ongoing testing should be used to monitor performance.
Uncertainty estimation is important in operational flood forecasting (see Beven, 2012, 289-311).Other research investigating probabilistic forecasting within FEWS includes Weerts et al. (2011) describing the quantile regression method.The FEWS data visualisation screen presents uncertainty to the end user by displaying the 50th percentile forecast together with the 5 to 95th percentile range.Large uncertainty clearly communicates a volatile catchment or poorly performing model.

Overview of FEWS
The FEWS open shell framework links a central database of hydrological and meteorological time series to an evergrowing inventory of forecasting models (Werner et al., 2013).An advanced set of utilities encapsulates data using extensible markup language (XML) and distributes these in a format specified by the model's adapter module.The end user is supplied with a graphical interface to view and interact with data.
The FEWS database adapter performs sophisticated merging of time series to ensure that, where multiple pieces of information are available, observational data takes precedence followed by the most recently generated forecast values (Deltares, 2010).This is important where upstream DBM nodes cascade forecasts to downstream nodes in order to achieve the maximum catchment lead time.In addition, FEWS substitutes the most logical value if a data point is missing or marked as corrupt.This improves operational robustness.
The FEWS framework and DBM module are purely for operational use.We used the open source R language for the DBM module as it has all necessary functionality and the liberal terms of the GNU General Public Licence allow straightforward, cost-effective distribution (R Development Core Team, 2008).Model identification, estimation and optimisation are carried out separately using the modeller's tools of choice.
The FEWS DBM module specifies the directory tree and naming convention shown in Fig. 3.Each DBM node is fully described by the components in Table 1.This information is stored as name-value pairs in each node's initialConditions file.The DBM model's MAIN SCRIPT calls functions in turn that transfer the appropriate values from initialConditions to the mainProcessingLoop.This function performs Eqs.(A4), (A6), and (A9).At present, the mainProcessingLoop function is identical in each node and could have been a shared function; however, for maximum flexibility, we chose to maintain a copy in each node so each can run a unique version of the Kalman filter algorithm in the future if required.

How the DBM module operates
The FEWS DBM nodes operate recursively: As sample period k becomes k − 1, new observational data (u k , y k ) are assimilated to produce the output forecast horizon, ŷk+1 to ŷk+δ , where δ is the maximum forecast lead time available at the node (see Appendix A).The algorithm also estimates the uncertainty of each point in the forecast horizon, Var( ŷk+1|k , . . ., ŷk+δ|k ).The nodes must be evaluated in the order determined by the model cascade.Details of the algorithm are provided in Appendix A.
Operational flood forecasting models must be able to stop and start execution between observations.The DBM model suits this mode of operation because it is a Markov process with a small number of model states.Storing the present model state is all that is required to form the next state and forecast horizon once new input data arrive.A small binary file can conveniently store the state data between samples.Also stored are the model parameters listed in Table 1; these are presently time invariant, but the storage-retrieval cycle adds little overhead and would allow time-varying parameters in the future.
At each time step FEWS extracts three data files from the central database: (1) observed input data made of present and δ past values -for multiple input sites this file will have as many columns as required; (2) output data containing present and δ future time steps, with future values marked as missing; and (3) auxiliary data (presently unused but available for future development).FEWS shifts the input data forward in time by δ time steps so the DBM module can always assume a δ of 0. The time shift appears complicated but is easily performed by FEWS, making the DBM processing straight forward.
A forecast is generated for any segment of the output series marked as missing; therefore, a forecast is always made for the final δ output values beyond the present sample time.The DBM module returns a file containing date, time, and a set of forecast percentiles (5, 50 and 95th by default).The FEWS adapter transfers this to the central database.Figure 4 shows the order of operation for a single time step.The process takes only seconds.The algorithm performs checks and logs progress at each stage of execution.The local machine can be powered down between samples if required.A Work directory stores input, output and log files.The states directory stores each node's internal state data between execution steps.Individual scripts are labelled with "F" for a function, "D" for a data file, "T" for a text file, or "M" for the main executable script.mechanistic interpretation: We consistently identified two dominant modes for rainfall to level models implying parallel fast and slow pathways (on the order of hours and days, respectively).Level to level nodes each returned a simple first-order response.The input non-linearity function for rainfall to level nodes roughly followed a logarithmic curve, indicating sensitivity to catchment antecedent conditions (increased runoff generation per unit of rainfall for a wet catchment).
We incorporated each node into the data assimilation algorithm shown in Appendix A and optimised the Kalman filter hyperparameters according to Smith et al. (2012a), calculating cost at the maximum forecast lead time.We manually specified adaptive gain hyperparameters to produce low-frequency, small amplitude forecast gain adjustments.Finally, we packaged the nodes into the full nested catchment structure.Example results for Sheepmount (Carlisle) are shown in Fig. 5.

Testing
The Environment Agency supplied testing data up to December 2011.This includes a significant flood in November 2009.We drove the DBM module with this data simulating the online data assimilation and forecasting process.The results for each location and lead time were stored.Figures 6 Figures comparing forecast and observation using a compressed time axis can appear overly convincing; instead we provide two alternative analyses of the Sheepmount data: (1) probability of detection (POD) and false alarm ratio (FAR) statistics (see Appendix B), and (2) Nash-Sutcliffe efficiency measures for the full data set and the sub-set above 3 m level.
Using the Environment Agency flood warning thresholds for Sheepmount, we counted 57 threshold-crossing events.Table 2 gives a summary of POD and FAR results for a 6h forecast at Sheepmount.The DBM module matched each threshold-crossing event.Data assimilation ensures the forecast tracks the observations reasonably well; occasionally the  forecasts stray above a threshold not crossed by observations, increasing the false alarm rate.As expected, using the upper percentiles of the uncertainty range in the FAR calculation increases the false alarm ratio.Table 3 gives Nash-Sutcliffe efficiency measures for the entire test data (35007 samples) and for the sub-set with a level over 3 m (334 samples).The global goodness-of-fit measure is excellent; however, it is considerably degraded for high-flow events.In theory, approximately 32 and 5 percent of observations should fall outside the 1 and 2 standard deviation ranges, respectively.The global data roughly conforms to this expectation; again, the high-flow events demonstrate higher variance.While the heteroskedastic forecast uncertainty estimates go some way to addressing this, there are clearly further challenges.However, as Figs.6 and 7 show, the DBM module provides valuable forecast information.It is worth noting the efficiency measure is reduced for highflow events by a general tendency to over-predict level, which is arguably preferable from a precautionary perspective.

Conclusions
We have shown the DBM method produces an effective and computationally efficient semi-distributed catchment model for real-time flood forecasting.The data assimilation process is based on the Kalman filter and is therefore probabilistic.The DBM module provides an estimate of forecast uncertainty.However, this study clearly illustrates the challenge of probabilistic flood forecasting by showing how the upper percentiles of the uncertainty range generate a high false alarm rate while often underestimating the true variance of highflow events.
We have also shown how the Delft-FEWS framework provides a flexible and extensible mechanism for incorporating new flood forecasting models.Developing the DBM module for FEWS provided a test case for knowledge transfer between the research and operational environments.This process was quite straightforward, which demonstrates the strength of the FEWS philosophy.
where B(z −1 ) and A(z −1 ) are polynomials of order m and n, respectively: Here z −1 is the discrete time backwards shift operator, i.e. z −i u k = u k−i ; δ is an integer value representing the time lag of the system, i.e. output (y) at time k is a response to the input stimulus (u) applied at time k − δ; ξ k is a noise input representing all the stochastic components of the system not accounted for by the model.The model orders (m and n) and the polynomial coefficients are identified and estimated from time series data (see Young, 2011).

A2 Mechanistic interpretation
The steady-state gain and time constant of a first-order transfer function, b 0 /(1 + a 1 z −1 ), are calculated respectively by b 0 /(1 + a 1 ) and − t/ ln(a 1 ), where t is the discrete sampling time (60 min in this example).A high-order transfer function response can be exactly matched by a combination of first-order transfer functions.The number of acceptable combinations grows rapidly as the order of the parent transfer function increases.For example, the second-order Eden transfer functions could be represented as two firstorder components in either parallel or feedback formation; however, the feedback form has no mechanistic justification.The first-order properties tell the modeller if he or she has identified a sensible representation of the process: In flood forecasting, steady-state gains should be positive, and highorder systems should form first-order components with realvalued time constants.The DBM approach calls for mechanistic interpretation -this guards against overfitting and poor extrapolation.

D. Leedal et al.: Application of data-based mechanistic modelling in the Eden catchment
how quickly the gain is allowed to adjust; ŷk is the estimate of level prior to applying the adaptive gain; and y k is the observed level.In the Eden example described in this paper, the q g terms were tuned manually taking values between 1 × 10 −2 and 1 × 10 −5 , depending on the characteristics of the node.It is useful to note that the adaptive gain method can be used on its own as a form of data assimilation for any deterministic model (e.g.Smith et al., 2012b).

Appendix B POD and FAR calculation
Probability of detection (POD) and false alarm ratio (FAR) (see Sene et al., 2009, p. 25) gauge a model's forecast skill.When a threshold crossing is forecast, that prediction can turn out to be either (a) correct, or (b) wrong.Alternatively, a forecast of a threshold not crossed can turn out to be (c) wrong, or (d) correct.POD and FAR are calculated by summing each outcome over a number of events, and then POD = a a+c and FAR = b b+a .

D.
Leedal et al.: Application of data-based mechanistic modelling in the Eden catchment Allen et al. (2010) describes the geology and hydrogeology of the region as part of the Eden Demonstration Test Catchment Project.

Fig. 2 .
Fig. 2. Configuration of the nodes making up the Eden FEWS DBM module.

Fig. 3 .
Fig. 3.The FEWS DBM module makes use of a directory and file naming convention.A catchment has a top level DBM <name> directory.Inside this, the functions of each node are contained within a separate Folder<ID> directory.The Config directory contains the main executable script and common processing functions.A Work directory stores input, output and log files.The states directory stores each node's internal state data between execution steps.Individual scripts are labelled with "F" for a function, "D" for a data file, "T" for a text file, or "M" for the main executable script.

Fig. 4 .
Fig. 4. Flowchart showing the operations performed during a single data assimilation and forecast cycle.The initialization step loads the k = 0 conditions.

Fig. 5 .
Fig. 5. Results for the January 2005 event at Sheepmount (Carlisle).Dots show observed levels spaced at 2-h intervals; the solid line is the f -step ahead forecast (where f = 2, 4, 6 and 7 h); the grey area is the ±2 standard deviation estimate of forecast uncertainty.

Fig. 6 .Fig. 7 .
Fig.6.Results the November 2009 event at Sheepmount (Carlisle).Dots show observed levels spaced at 2-h intervals; the solid line is the f -step ahead forecast (where f = 2, 4, 6 and 7 h); the grey area is the ±2 standard deviation estimate of forecast uncertainty.

Table 1 .
The data structure describing a FEWS DBM node.Each node also includes a hard-coded input non-linearity function.

Table 3 .
Summary of model performance for observations above 0 m (row 1) and 3 m (row 2).Column 3 (N-S) is the Nash-Sutcliffe efficiency measure.Columns 4 and 5 show the percentage of observations outside the estimated 1 and 2 standard deviation range.Results are for the 6-h forecast at Sheepmount.