DECIPHeR v1: Dynamic fluxEs and ConnectIvity for Predictions of HydRology

This paper presents DECIPHeR (Dynamic fluxEs and ConnectIvity for Predictions of HydRology), a new model framework that simulates and predicts hydrologic flows from spatial scales of small headwater catchments to entire continents. DECIPHeR can be adapted to specific hydrologic settings and to different levels of data availability. It is a flexible model framework which includes the capability to (1) change its representation of spatial variability and hydrologic connectivity by implementing hydrological response units in any configuration and (2) test different hypotheses of catchment behaviour by altering the model equations and parameters in different parts of the landscape. It has an automated build function that allows rapid set-up across large model domains and is open-source to help researchers and/or practitioners use the model. DECIPHeR is applied across Great Britain to demonstrate the model framework. It is evaluated against daily flow time series from 1366 gauges for four evaluation metrics to provide a benchmark of model performance. Results show that the model performs well across a range of catchment characteristics but particularly in wetter catchments in the west and north of Great Britain. Future model developments will focus on adding modules to DECIPHeR to improve the representation of groundwater dynamics and human influences.

1 Introduction

Water resources require careful management to ensure adequate potable and industrial supply, to support the economic and recreational value of water, and to minimize the impacts of hydrological extremes such as droughts and floods on the economy, river ecosystems and human life. Robust simulations and predictions of river flows are increasingly needed across multiple temporal and spatial scales to support such management strategies (Wagener et al., 2010) that may range from the assessment of local field-scale flood mitigation measures to emerging water challenges at regional to continental scales (Archfield et al., 2015). Such approaches are particularly important, indeed mandated, given national and international policies on water management, such as the European Union's Water Framework Directive (EC, 2000) and Floods Directive (EC, 2007). Specifically national and international information on water resources and low and high flows is needed to underpin robust environmental management and policy decisions. This requires the effective integration of field observations and numerical modelling tools to provide tailored outputs at gauged and ungauged locations across a wide range of scales relevant to policymakers and societal needs.

To address this need, a fundamental challenge for hydrologic sciences is to develop hydrological models that represent the complex drivers of catchment behaviour, such as space- and time-varying climate, land cover, human influence, etc. (Blöschl and Sivapalan, 1995). The hydrologic community has made substantial investments to develop and apply hydrological models over the past 50 years to produce simulations and predictions of surface and groundwater flows, evaporation, and soil moisture storage across multiple scales. These include gridded approaches (e.g. PCRaster Global Water Balance, PCR-GLOBWB – Wada et al., 2014; Variable Infiltration Capacity, VIC – Hamman et al., 2018; Liang et al., 1994; Grid-to-Grid model– Bell et al., 2007; multiscale hydrologic model – Samaniego et al., 2010; DK model – Henriksen et al., 2003), semi-distributed approaches that aggregate the landscape into hydrologic response units or sub-catchments (e.g. HYPE – Lindström et al., 2010; SWAT – Arnold et al., 1998; TOPNET – Clark et al., 2008b) and many conceptual models applied at the catchment scale (Beven and Kirkby, 1979; Burnash, 1995; Coron et al., 2017; Leavesley et al., 1996; Lindström et al., 1997; Zhao, 1984). The current generation of hydrological models can represent a range of natural and anthropogenic processes and various levels of spatial complexity. Furthermore, there are significant ongoing efforts to represent spatial heterogeneity at finer scales over national–global scales (Bierkens et al., 2015; Wood et al., 2011) and build multi-model frameworks to test competing hypotheses of catchment behaviour, such as FUSE (Clark et al., 2008a) and SUPERFLEX (Fenicia et al., 2011; Kavetski and Fenicia, 2011).

However, whilst these models have provided a wealth of useful insights and relevant outputs, they tend to have a fixed representation of spatial variability (i.e. a single spatial resolution or a single spatial structure such as a raster-based one), lack spatial connectivity between hillslope-to-hillslope and hillslope and riverine components, be computationally expensive, and/or employ a single model structure across the model domain or nested catchment scale. This impacts our ability to apply models to a wide range of scales, places and water challenges, as different model representations of hydrological processes (i.e. model structure, parameterizations, hydrologic connectivity or spatial variability) are needed to capture heterogeneous hydrological responses and changing landscape connectivity, particularly for local conditions. Consequently, there is a pressing need to develop new spatially flexible modelling tools that can be applied to a range of space and timescales and that are based on general hydrological principles applicable to a broad spectrum of different catchment types. The need for such approaches is well-documented in the literature (Clark et al., 2011, 2015; Mendoza et al., 2015), with calls for flexible hydrological modelling systems that can (1) incorporate different model structures and parameterizations in different parts of the landscape to represent a variety of processes; (2) change their spatial complexity, variability and/or hydrologic connectivity for hillslope elements and the river (Beven and Freer, 2001; Mendoza et al., 2015); and (3) be applied across a wide range of spatial and temporal scales and across places (Blöschl et al., 2013). However, few such models exist.

In line with these requirements, we have created a new model framework, DECIPHeR (Dynamic fluxEs and ConnectIvity for Predictions of HydRology), to simulate and predict hydrologic flows and connectivity from spatial scales of small headwater catchments to entire continents. The flexible modelling framework allows users to test different spatial resolutions, spatial configurations (i.e. gridded, semi-distributed or lumped), levels of hydrologic connectivity (i.e. representations of the lateral fluxes of water across model elements) and process representation (i.e. model structure and parameters). DECIPHeR has an automated build function that allows rapid set-up across required model domains with limited user input. The underlying code has been optimized to run large ensembles and enable model uncertainty to be fully explored. This is particularly important given inherent uncertainties in hydro-climatic datasets (Coxon et al., 2015; McMillan et al., 2012) and their impact on model calibration, regionalization and evaluation (Freer et al., 2004; Kavetski et al., 2006; Kuczera et al., 2010; McMillan et al., 2010, 2011; Westerberg et al., 2016). We have specifically made the model code readable, reusable and open-source to allow the broader community to learn from, verify and advance the work described here (Buytaert et al., 2008; Hutton et al., 2016).

In this paper, we (1) describe the key capabilities and concepts that underpin DECIPHeR, (2) provide a detailed discussion of the model code and components, (3) demonstrate its application at the national scale to 1366 catchments in Great Britain (GB), and (4) discuss potential future model developments.

2 The DECIPHeR modelling framework

2.1 Key concepts

The DECIPHeR modelling framework is based on the key concepts enshrined in Dynamic TOPMODEL originally introduced by Beven and Freer (2001). Since its original development, Dynamic TOPMODEL has been applied in a wide range of studies (Freer et al., 2004; Liu et al., 2009, p. 200; Metcalfe et al., 2017; Page et al., 2007; Younger et al., 2008) and integrated into other modelling frameworks (e.g. HydroBlocks – Chaney et al., 2016). The core ideas of Dynamic TOPMODEL were threefold (Beven and Freer, 2001): (1) to allow more flexibility in the definition of similarity in function for different points in the landscape, (2) to implement a non-linear routing of subsurface flow that simulates dynamically variable upslope subsurface contributing area and (3) to remain computationally efficient so that uncertainty in hydrological simulations can be estimated.

To realize this, Dynamic TOPMODEL uses hydrological response units (HRUs) to group raster-based information into non-contiguous spatial elements in the landscape that share similar characteristics (see Fig. 1). Each HRU maintains hydrological connectivity in the landscape via weightings that determine the proportions of lateral subsurface flux from each HRU to all connected HRUs and flows to river cells. This solution offers key advantages in capability to traditional grid-based or lumped approaches employed by many hydrological models. Firstly, the user can split up the catchment using, for example, different landscape attributes (e.g. geology and land use) and/or spatially varying inputs (e.g. rainfall, evaporation, etc.) to define spatial similarity. This capability allows the user to modify the spatial complexity, resolution and/or hydrologic connectivity of hillslope elements and river network reaches in any configuration. Secondly, each HRU is treated as a separate functional unit in the model which can have different process conceptualizations and parameterizations. This means that more process complexity can be incorporated where needed to better suit local conditions (e.g. to account for “point-source” human influences or more complex hydrological processes such as surface–groundwater exchanges). Finally, by grouping together similar parts of the landscape, HRUs minimize run times of the model compared to grid-based or fully distributed formulations while still allowing model simulations to be mapped back into space.

Figure 1 Digital terrain analysis and simplified examples of using classification layers to discretize a hypothetical catchment into hydrological response units, from (a) the gauge network; (b) landscape layer with a chalk outcrop for HRU 2; (c) the gauge network, ungauged flow point and landscape layer; and (d) same as (c), with individual river reach lengths specified.

While these key concepts that underpin Dynamic TOPMODEL address many of the challenges outlined in the introduction, for the most part the model has only ever been applied to a single catchment or very simple nested catchments in headwater basins (Peters et al., 2003). Consequently, we have completely restructured and rewritten the model code and added several new features to improve the flexibility and automation of the original Dynamic TOPMODEL code so that the model can be applied from single small headwater catchments to regional, national and continental scales. These changes include the following:

1. Both the legacy and new model code have been updated to a version that is compliant with Fortran 2003, with new array and memory handling to allow significantly larger and more complex gauging networks to be processed.

2. The model build process is now fully automated to allow national- continental-scale data to be easily and quickly processed and to build and apply models in complex multi-catchment regions.

3. New model code and functions have been written to (a) enable greater flexibility in the complexity and spatial characteristics of river network and routing properties (a newly developed river network scheme allows flow simulations to be produced for any gauged or ungauged point on a river network, and allows river reaches of any length for individual hillslope–river flux contributions), (b) ensure that multiple points on the river network can be initialized via local storage and fluxes in each HRU successfully, and (c) seamlessly facilitate digital terrain analysis (DTA) classification layers and result in a rainfall–runoff model configuration that allows each individual HRU to have a different model structure, parameters and climatic inputs.

4. A new analytical solution of the subsurface flow equations has been implemented, resulting in increased computational speed and numerical stability.

5. The model can be easily adopted and adapted because it is open-source, version controlled and includes a detailed user manual.

HRUs are defined prior to rainfall–runoff modelling, and DECIPHeR consists of two key steps where (1) digital terrain analyses are performed to define the gauge network, set up the river network and routing, discretize the catchment into HRUs, and characterize the spatial variability and hydrologic connectivity in the landscape and (2) HRUs are run in the rainfall–runoff model to provide flow time series. These two steps are described in the following sections. More detailed descriptions of the input and output files, code workflows, and codes can be found in the user manual.

2.2 Digital terrain analysis (DTA)

The DTA in DECIPHeR constructs the spatial topology of the model components to define hillslope and riverine elements. The DTA defines the spatial extent of every HRU based upon multiple attributes, quantifies the connectivity between these HRUs in the landscape, determines the river network and all downstream routing properties, and determines the extent to which and where simulated output variables (i.e. discharge) should be produced (including gauged or ungauged locations; see Fig. 1).

2.2.1 Data prerequisites

The minimum data requirement to run the DTA is a digital elevation model (DEM) and XY locations where flow time series are needed on the river network. The DEM must contain no sinks or flat areas to ensure that the river network and catchments can be properly delineated, as is common in digital terrain analyses. This means that any real inland sinks (such as lakes) will be filled. Accounting for these features in the modelling framework will be a focus for future model development.

Additional data can also be incorporated depending on data availability and modelling objectives. A river network can be supplied if the user wishes to specify headwater cells from a predefined river network, and reference catchment areas and masks can be used to identify the best station location on the river network. Depending on user requirements, topographic, land use, geology, soils, and anthropogenic and climate attributes can be supplied to define the spatial topology and thus differences in model inputs, structure and parameterization.

2.2.2 River network, catchment identification and river routing

DECIPHeR generates streamflow estimates at any point on the river network specified by the user. A river network is generated in DECIPHeR which matches the DEM flow direction and always connects to the boundary of the DEM or the sea. The river network is created from a list of headwater cells, which the functions can use and/or produce in three different ways depending on user requirements and/or data availability:

1. A list of predefined headwater (i.e. starting) river locations is read into the DTA algorithms from a file.

2. Headwater cells are found from a predefined river network.

3. Where no predefined river network or headwater locations are available, headwater cells are found from a river network which is derived from cells that meet thresholds of the accumulated area and/or topographic index.

Each headwater location is then routed downstream in a single flow direction via the steepest slope until reaching a sea outlet, other river or edge of the DEM to construct a contiguous river network for the whole area of interest. Gauge locations are then generated on the river network from the point locations specified by the user. If a reference catchment mask or area is available, catchment masks are produced for candidate river cells found in a given radius, and the catchment mask with the best fit to the reference mask or area is chosen as the gauge location. Otherwise the closest river cell is chosen as the gauge location.

Catchment masks are created from the final gauge list, with both individual masks for all the points specified on the river network and a combined catchment mask with the nested catchment masks created for use in the creation of the hydrological response units. From the river network and gauge locations, the river network connectivity is derived, with each river section labelled with a unique river ID. A suite of routing tables is also produced so that each ID knows its downstream connections and to allow multiple routing schemes to be configured (see Sect. 2.3.4 for a description of the current routing scheme implemented in the modelling framework). These codes also provide the option of setting a river reach length where output time series can also be specified at different reach lengths between gauges (see Fig. 1, HRU Set-up D).

2.2.3 Topographic analysis

Topography, slope, accumulated area and the topographic index are important properties of the landscape to aid the definition of hydrologic similarity and more dominant flow pathways. In DECIPHeR, they provide the basis for river routing and river network configuration, and they also can be used to help determine the initial separation of landscape elements for defining hydrological similarity using percentiles of accumulated area, elevation and slope (in addition to alternative catchment attributes such as urban extent, geology, land use, soils, etc.).

The topographic index is calculated using the M8 multiple-flow directional algorithm of Quinn et al. (1995). The DTA calculates slope, accumulated area and the topographic index for the whole domain. It uses the river mask to define the cells where accumulated area cannot accumulate downstream and the catchment mask to ensure that the accumulated area does not accumulate across nested catchment boundaries.

2.2.4 Hydrological response units

The most critical aspect of running DECIPHeR is to define HRUs according to user requirements. The HRU configuration determines the spatial connectivity and complexity of model conceptualization as well as the spatial variability of inputs and conceptual structure and parameters to be implemented in each part of the landscape. Any number of different spatial discretizations can be derived and subsequently applied in the DECIPHeR framework, allowing the user to experiment with different model structures and parameterizations and modify representations of spatial variability and hydrologic connectivity.

In the DTA, hydrologically similar points in the landscape are grouped together so that each HRU is a unique combination of four different classification layers. These specify (1) the initial separation of landscape elements from topographic information (e.g. slope, accumulated area and/or elevation), (2) inputs, (3) process conceptualizations, and (4) parameters implemented for each HRU store in the model (see Fig. 2). These classification layers can be derived from climatic inputs, such as spatially varying rainfall and potential evapotranspiration, and landscape attributes such as geology, land use, anthropogenic impacts, soil data, slope and accumulated area. The simplest set-up will consist of one HRU per catchment, while the most complex set-up can consist of one HRU for every grid cell (i.e. fully distributed).

Figure 2 DECIPHeR represents spatial heterogeneity in the landscape through hydrological response units (HRUs). Each HRU can have a different model structure, parameters or inputs.

To maintain hydrological connectivity in the landscape, the proportions of flow between the cells comprising each HRU are calculated based on accumulated area and slope. The flow fractions are then aggregated into a flow distribution matrix that summarizes the proportions (weightings) of lateral subsurface flow from each HRU to (1) itself, (2) another HRU or (3) a river reach. For n hydrological response units, the weights (W) are defined as

$$\begin{array}{}\text{(1)}& {\displaystyle }W=\left(\begin{array}{ccc}{w}_{\mathrm{1},\mathrm{1}}& \mathrm{\cdots }& w_{\mathrm{1},n}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ w_{n,\mathrm{1}}& \mathrm{\cdots }& w_{n,n}\end{array}\right),\end{array}$$

where each row defines how the HRU's output is distributed to other HRUs, any river reaches or itself and each column represents the total input to each HRU at every time step as the weighted sum of all the upstream outputs, each row and column sum being 1 to ensure mass balance. The weights are detailed in an HRU flux file (which is fixed for a simulation) as a flow distribution matrix along with tabulated HRU attributes to provide information on which inputs, parameter and model structure type each HRU is using.

2.3 Rainfall–runoff modelling

2.3.1 Data pre-requisites

To run the rainfall–runoff modelling component of DECIPHeR, time-series forcing data of rainfall and potential evapotranspiration are required. Discharge data can also be provided for gauged locations and are used to initialize the model.

Besides forcing data, the model also needs (1) the HRU flux file and routing files produced by the DTA, (2) a parameter file specifying parameter bounds for Monte Carlo sampling of parameters, and (3) project and setting files specifying the number of parameter sets to run, which HRU and input file to use, etc.

2.3.2 Initialization

Initialization is an important step for any rainfall–runoff model to ensure that subsurface flows, storage and the river discharge have all stabilized. This can be particularly problematic when modelling regionally over a large area, as not all HRUs will initialize at the same rate (depending on size and slope characteristics).

A simple homogenous initialization is currently implemented in DECIPHeR, where the storage deficits for all HRUs are determined from an initial discharge. This is calculated as a mean area weighted discharge of the starting flows at time step 1 for all output points on the river network. If a gauge does not have an initial flow, then the initial flow is either calculated from the mean of the data or set to a value of 1 mm d⁻¹ (as a representative starting flow for most catchments) if no flow data are available. The initial discharge is assumed to be solely due to the subsurface drainage into the river, so it is used as the starting value for Q_SAT (subsurface flow) and to determine the associated storage and unsaturated zone fluxes. The model is then run for an initialization period to allow its model stores and fluxes to fully stabilize with the catchment climatic information. Initialization periods depend in part on the parameterization of the model simulation run as well as the size and characteristics of the catchment being considered.

2.3.3 Parameters

DECIPHeR can be run either using default parameter values or through Monte Carlo sampling of parameters between set parameter bounds to produce ensembles of river flows. In the DTA, the user can set different parameter bounds for each HRU or sub-catchment, thus specifying areas of the landscape where different parameter bounds may be needed. Alternatively, a single set of parameter bounds can be applied across the model domain.

For the model structure provided in the standard build and described below, there are seven parameters that can be sampled or set to default parameters. These parameters describe the transmissivity of the subsurface, the water holding capacity and permeability of soils, and the channel routing velocity (see Table 1). More parameters can easily be added by the user if required for different model structures by changing the model source code.

Table 1 Overview of DECIPHeR's stores, fluxes and parameters.

2.3.4 Model structure

The description below details the model structure that is provided in the open-source code (see Fig. 3 and Table 1). While the code is built to be modular and extensible so that a user can easily implement multiple model structures if so wished, the aim of this paper and the initial focus of the code development was on applying the model across large scales and beginning with a release that has relatively simple representations of the core processes. Thus, we provide a single model structure in the open-source code that serves as a model benchmark to be built upon in future iterations.

Figure 3 Simplified conceptual diagram of the model structure currently implemented in DECIPHeR. All scientific notations are described in Table 1.

The model structure consists of three stores defining the soil profile (S_RZ, S_UZ and S_D in Fig. 3), which are implemented as lumped stores for each HRU. The first store is the root zone storage (S_RZ). Precipitation (P) is added to this store, and then evapotranspiration (ET) is calculated and removed directly from the root zone. The maximum specific storage of S_RZ is determined by the parameter SR_max. Actual evapotranspiration from each HRU depends on the potential evapotranspiration (PET) rate supplied by the user and the root zone storage using a simple common formulation where evapotranspiration is removed at the full potential rate from saturated areas (i.e. if the root zone storage is full) and at a rate proportional to the root zone storage in unsaturated areas:

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{ET}=\text{PET}\cdot \left(S_{\mathrm{RZ}}/{\text{SR}}_{max}\right).\end{array}$$

Once the root zone reaches maximum capacity (i.e. deficit of zero and conceptually analogous to field capacity), any excess rainfall input is added to the unsaturated zone (S_uz), where it is routed to the saturated zone (S_D). If the saturated zone is also full (as determined by S_max), Q_EXUS is added to the saturation excess storage (S_EX) and routed directly overland as saturated excess overland flow (Q_OF). The unsaturated zone links the S_RZ and saturated zones according to a linear function that includes a gravity drainage time delay parameter (T_d) for vertical routing through the unsaturated zone. The drainage flux (Q_uz) from the unsaturated zone to the saturated zone is at a rate proportional to the ratio of unsaturated zone storage (S_uz) to the storage deficit (S_D):

$$\begin{array}{}\text{(3)}& {\displaystyle }{Q}_{\mathrm{UZ}}=S_{\mathrm{UZ}}/\left(S_{\mathrm{D}}\cdot T_{\mathrm{d}}\right).\end{array}$$

Changes to storage deficits for each HRU are dependent on recharge from S_UZ(Q_UZ), fluxes from upslope HRUs (Q_IN) and downslope flow out of each HRU (Q_SAT), with subsurface flows for each HRU being distributed according to the DTA flow distribution matrix described in Sect. 2.2.4:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{D}}}{\mathrm{d}t}}=Q_{\mathrm{SAT}}-Q_{\mathrm{IN}}-Q_{\mathrm{UZ}},\end{array}$$

where S_D is the current deficit in the saturated zone, Q_SAT is outflow from this HRU, Q_IN is inflow into the HRU representing subsurface flow from other HRUs and Q_UZ is inflow into the HRU representing drainage from the unsaturated zone of this HRU. This equation is solved sequentially for each HRU and provides values for the deficit S_D and outflow Q_SAT at time step t for each HRU. In DECIPHeR, this equation is solved analytically (see Appendix for derivation of this solution), assuming a transmissivity profile that declines exponentially with depth and is truncated at depth S_max such that no flow is generated when the deficit is greater than S_max (Beven and Freer, 2001). The analytical solution provides better computational speed and increased numerical stability compared to the iterative four-point numerical scheme described by Beven and Freer (2001).

The exponential transmissivity profile takes the following shape (Beven and Freer, 2001; Eq. 6):

$$\begin{array}{}\text{(5)}& {\displaystyle }{Q}_{\mathrm{SAT}}=T_{\mathrm{0}}\tan \mathit{\beta }\exp (-fz)=Q_{\mathrm{0}}\exp \left(-\stackrel{\mathrm{‾}}{S}\left./\text{SZM}\right)\right..\end{array}$$

The truncated exponential transmissivity profile takes the shape (rewritten from Beven and Freer, 2001; Eq. 9):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{Q}_{\mathrm{SAT}}=\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{Q}_{\mathrm{0}}\cos \mathit{\beta }[\exp \left(-\cos \mathit{\beta }\stackrel{\mathrm{‾}}{S}/\text{SZM}\right)-\exp & \\ \left(-\cos \mathit{\beta }{S}_{max}/\text{SZM}\right)]& \stackrel{\mathrm{‾}}{S}\le S_{max}\\ \mathrm{0}& \stackrel{\mathrm{‾}}{S}>S_{max}\end{array}\right.,\end{array}$$

where β is the mean slope of the HRU and $\stackrel{\mathrm{‾}}{S}$ is the average deficit across the HRU. The parameter, SZM, sets the rate of the exponential decline in saturated zone hydraulic transmissivity with depth thereby controlling the shape of the recession curve in time. The parameter, S_max, sets the saturated zone deficit threshold at which downslope flow between HRUs no longer occurs. If the storage deficit is less than zero (i.e. the soil is at or above its saturation capacity), then excess storage (Q_EXS) is added to saturation excess overland flow (Q_OF). Q₀ is the maximum rate of Q_SAT from an HRU when the HRU is at saturation and is calculated from

$$\begin{array}{}\text{(7)}& {\displaystyle }{Q}_{\mathrm{0}}={\displaystyle \frac{T_{\mathrm{0}}}{e^{\mathit{\lambda }}}},\end{array}$$

where the parameter T₀ determines the lateral saturated hydraulic transmissivity at the point when the soil is saturated and λ is the average topographic index across the HRU.

Channel flow routing in DECIPHeR is modelled using a set of time delay histograms that are derived from the digital terrain analyses for the points where output is required. A fixed channel wave velocity (CHV) is applied throughout the network to account for delay and attenuation in the simulated flows (Q_SIM). DECIPHeR is a mass conserving model, and therefore the model water balance always closes (subject to small rounding errors).

2.4 Model implementation

The DECIPHeR model code is available on GitHub (https://github.com/uob-hydrology/DECIPHeR, last access: 28 May 2019) and is accompanied by a user manual which provides a detailed description of the file formats, how to run the codes and a code workflow. All the model code is written in Fortran for its speed, efficiency and ability to process large scale spatial datasets. Two additional bash scripts are provided as an example of calling the digital terrain analysis codes.

3 Great Britain national model implementation and evaluation

While the modelling framework has a wide range of functionality, in this paper we wanted to demonstrate the ability of the model to be applied across a large domain to generate ensembles of flows at thousands of gauging stations and evaluate its current capability across large scales to guide future model developments. Consequently, we applied DECIPHeR to 1366 gauges in Great Britain (GB), and in this section we describe the model set-up, input data, evaluation criteria and model results.

3.1 Great Britain hydrology

Catchments in Great Britain (GB) cover wide hydrologic and climatic diversity. Hydro-climatic characteristics were derived from rainfall, potential evapotranspiration and flow data described in Sect. 3.3.1. Figure 4 shows the mean annual rainfall, mean annual potential evapotranspiration, runoff coefficient, and slope of the flow duration curve between the 30 and 70 flow percentiles for the 1366 catchments in this study. Rainfall is highest in the west and north of GB and lowest in the east and south, ranging from 540 to 3400 mm yr⁻¹ (Fig. 4a), while potential evapotranspiration is highest in the east and south and lowest in the west and north, ranging from 370 to 545 mm yr⁻¹ (Fig. 4b). This regional divide of rainfall and potential evapotranspiration is reflected in the runoff coefficients (Fig. 4c), where generally runoff coefficients are lowest in the east and south and highest in the north and west. The slope of the flow duration curve (Fig. 4d) is a more mixed picture across GB, with lower values (i.e. a less variable flow regime) found in north-eastern Scotland, the Midlands and patches of the south-east and higher values (i.e. a more variable flow regime) in the west, with the highest values for ephemeral and/or small streams in the south-east.

Figure 4 Hydro-climatic characteristics of 1366 GB catchments. (a) Annual rainfall (mm yr⁻¹). (b) Annual potential evapotranspiration (mm yr⁻¹). (c) Runoff coefficient (–). (d) Slope of the flow duration curve between the 30th and 70th percentiles (–). Minimum and maximum values on colour bars have been chosen to show clear differences between catchments.

River flows vary seasonally, with the highest totals generally occurring during the winter months when rainfall totals are highest and evapotranspiration totals are lowest and the lowest totals during the summer months (April–September) resulting from lower precipitation totals and higher evapotranspiration losses due to seasonal variations in energy inputs. Snowmelt has little impact on river flows in GB except for some catchments in the Scottish Highlands where snowmelt contributions can impact the flows. River flow patterns are also heavily influenced by groundwater contributions from various regional aquifer systems. In catchments overlying the Chalk outcrop in the south-east of GB, flow is groundwater-dominated, with a predominantly seasonal hydrograph that responds less quickly to rainfall events. Land use and human influences also significantly impact river flows, with flows most heavily modified in the south-east and Midland regions of England due to high population densities.

Figure 5 Inputs and outputs of digital terrain analyses for GB: (a) 50 m hydrologically consistent digital elevation model, (b) DECIPHeR river network, (c) nested catchment mask, (d) topographic index and (e) 5 km input grid.

3.2 Digital terrain analyses for GB

To implement DECIPHeR across GB, the UK NEXTMap 50 m gridded digital elevation model was used as the basis of the Digital Terrain Analysis (Intermap Technologies, 2009). The first step was to ensure that the DEM contained no sinks or flat areas before being run through the DTA codes. Many freely available packages and codes exist to fill sinks in DEMs, but for use with large national datasets, a two-stage process is often necessary to ensure that there are no flat areas in the DEM and that important features, such as steep-sided valleys, are not filled due to pinch points in the DEM. For this study, we first applied an optimized pit-removal routine (Soille, 2004; code available on GitHub at https://github.com/crwr/OptimizedPitRemoval, last access 28 May 2019). This tool uses a combination of cut and fill to remove all undesired pits while minimizing the net change in landscape elevation. We then applied a sink fill routine to ensure that no flat areas remained in the DEM.

The inputs and outputs for the GB DTA are summarized in Fig. 5. To build the river network, we first extracted headwater cells from the Ordnance Survey MasterMap Water Network Layer, a dense national river vector dataset for GB. These headwater cells were then routed downstream via the steepest slope to generate the river network used by the model. This ensures that the DEM and the calculated stream network are consistent for flow accumulations based on surface slope. Locations of 1366 National River Flow Archive (NRFA) gauges were used to define the gauging network and specify points on the river network where output was required. We used NRFA catchment areas and masks as a reference guide to evaluate the best point for the gauge locations from potential river cell candidates within a local search area. Slope, accumulated area and the topographic index were then calculated for every grid cell and routing files produced.

Finally, we chose three classifiers to demonstrate the modelling framework while ensuring that the number of HRUs was still computationally feasible for modelling across a large domain:

1. The catchment boundaries for each gauge were used to ensure minimal fluxes across catchment boundaries.

2. A 5 km grid for the rainfall and potential evapotranspiration inputs was used to represent the spatial variability in climatic inputs across GB.

3. Three equal classes of slope and accumulated area were implemented, resulting in HRUs that cascade downslope to the valley bottom.

3.3 Rainfall–runoff modelling

3.3.1 Input and evaluation datasets

Daily data of precipitation, potential evapotranspiration and discharge for a 55-year period from 1 January 1961 to 31 December 2015 were used to run and assess the model. This period was chosen as an appropriate test for the model covering a range of climatic conditions and to demonstrate the model's ability to simulate long time periods within uncertainty analysis frameworks. The year 1961 was used as a warm-up period for the model; therefore no model evaluation was quantified in this period.

A national gridded rainfall and potential evapotranspiration product was used as input into the model. Daily rainfall data were obtained from the Centre for Ecology & Hydrology Gridded Estimates of Areal Rainfall dataset (CEH-GEAR; Keller et al., 2015; Tanguy et al., 2016).  This dataset consists of 1 km² gridded estimates of daily rainfall from 1961 to 2015 for Great Britain and Northern Ireland derived from the Met Office UK rain gauge network. The observed precipitations from the rain gauge network are quality controlled, and then natural neighbour interpolation is used to generate the daily rainfall grids. Daily potential evapotranspiration data were obtained from the Centre for Ecology & Hydrology (CEH) Climate hydrology and ecology research support system potential evapotranspiration dataset for Great Britain (CHESS-PE; Robinson et al., 2016). This dataset consists of 1 km² gridded estimates of daily potential evapotranspiration for Great Britain from 1961 to 2015 calculated using the Penman–Monteith equation and data from the CHESS meteorology dataset. Both datasets were aggregated to a 5 km grid as forcing for the national model run.

The model was evaluated against daily streamflow data for the 1366 gauges obtained from the National River Flow Archive (https://nrfa.ceh.ac.uk/, last access: 28 May 2019). These data are collected by measuring authorities including the Environment Agency (EA), Natural Resources Wales (NRW) and Scottish Environmental Protection Agency (SEPA) and then quality controlled before being uploaded to the NRFA site.

3.3.2 Model structure and parameters

To initially evaluate the model, DECIPHeR was run within a Monte Carlo simulation framework whereby 10 000 parameter sets were randomly sampled from a uniform prior distribution. This number of parameter sets was chosen to provide a reasonable sampling of the parameter space for demonstration purposes; however, for a full evaluation of the parameter space, more parameter sets would be needed.

These parameters were applied uniformly across the HRUs and used within a single model structure (as described in Sect. 2.3.4). Given the wide range of hydro-climatic conditions across GB, sampling of the feasible parameter space was ensured by using wide sampling ranges based on previous studies that have used Dynamic TOPMODEL (Beven and Freer, 2001; Freer et al., 2004; Page et al., 2007; Table 2).

Table 2 Parameter ranges.

3.3.3 Model evaluation

Daily time series of discharge for the 10 000 model simulations from each gauge were evaluated against daily observed flow for all 1366 gauges. This is a challenging test for the model, as these catchments cover a large range of hydrologic behaviour across GB and are impacted by a variety of climatic, geological and anthropogenic processes as outlined in Sect. 3.1. However, evaluating the model over such a large number of gauges acts as a benchmark of model performance and a means of identifying future areas for model development.

To benchmark model performance, we wanted to evaluate the model's ability to capture a range of hydrologic behaviour, including maintaining overall water balance, capturing flow variability, reproducing low and high flows, and observing the timing of flows. Consequently, multiple metrics, including hydrological signatures, standard hydrological model performance metrics and statistics of the flow time series were used to provide insights into model performance. Based on previous studies evaluating national-scale models (McMillan et al., 2016) and considering a diagnostic approach to model evaluation (Coxon et al., 2014; Gupta et al., 2008; Yilmaz et al., 2008), four metrics were chosen which are summarized in Table 3 alongside their equations: (i) the Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970), (ii) slope of the flow duration curve (Yadav et al., 2007), (iii) bias in runoff ratio (Yilmaz et al., 2008) and (iv) low flow volume (Yilmaz et al., 2008).

Table 3 Evaluation metrics used in the study.

These metrics are also used to determine a behavioural ensemble of parameter sets. The focus of this model application is to demonstrate that the model can be run in a Monte Carlo framework. Consequently, while many different approaches could be used to determine a behavioural ensemble of parameter sets (see, for example, Beven, 2006; Coxon et al., 2014; Krueger et al., 2010; Westerberg et al., 2011), in this study we adopt a simple approach to produce ensembles of flows. The four metrics described above are combined, and the behavioural ensemble was then taken as the top 1 % of the model simulations according to this combined score. To calculate the combined score, each metric was ranked in turn, these ranks were summed and all simulations were sorted by the total combined rank. Weaker and stricter performance thresholds in NSE and bias metrics were also defined to further explore the performance of the ensembles against a common set of criteria (see Table 3). These were chosen based on previous studies, and although subjective, the hydrological modelling community has yet to agree on benchmarks for the comparison of model performance (Seibert et al., 2018).

3.4 Results

3.4.1 Digital terrain analysis and model simulation

DECIPHeR was set up for GB, covering a total catchment area of 154 763 km² for 1366 gauges and 365 principal basins. Principal basin area ranged from 7.87 to 9935 km², with a median of 137 km². Using the HRU classifiers specified in Sect. 3.2, the number of HRUs contained within each principal basin ranged from 17 to 8978, with a median of 123 HRUs. HRU area ranged from 0.0025 to 14.33 km², with a median HRU area of 0.65 km².

In total, 13 660 000 55-year time-series flow simulations were produced. One simulation over the 55-year time period for the largest river basin (9935 km²), with 8978 HRUs, takes approximately 15 min to run on a standard CPU, outputting simulated discharge for all the 98 gauges that lie within the Thames at Kingston river basin. For the smallest river basin that has 17 HRUs and one river gauge, a single simulation over the 55-year time period on a standard CPU takes less than a second.

3.4.2 Overall model performance

Our first assessment of model performance is the overall model performance for the four performance metrics calculated from the 10 000 simulated daily flow time series produced for each gauge. Figure 6 shows the percentage of catchments that met the stricter and weaker performance thresholds defined in Table 3 from the entire ensemble of 10 000 model simulations and from the top 1 % behavioural ensemble generated from the combined ranking of the four metrics. Our results show that most catchments are able to meet both the performance thresholds. The vast majority of catchments (92 %) gain an NSE score greater than zero (i.e. better than mean climatology), and 80 % of the catchments gain an NSE score greater than 0.5. The model does well in reproducing low flow volumes and the slope of the flow duration curve, with most gauges (98 % and 96 % respectively) meeting the stricter performance threshold.

Figure 6 Percentage of catchments for each metric that meet the weaker and stricter performance thresholds for the entire ensemble of 10 000 model simulations and from the top 1 % behavioural ensemble of 100 model simulations generated from the combined ranking of the four metrics.

RRBIAS (bias in runoff ratio) evaluates the model's ability to reproduce water balance in the catchment; the current implementation of the model has to maintain mass balance, while many of the observed flow data for many of these catchments do not maintain mass balance due to inter-catchment groundwater flows, anthropogenic influences such as surface and ground water abstractions, or data errors (this is further discussed in Sect. 4.4.4). Consequently, RRBIAS is a more difficult metric for the model to capture, and this is reflected by the fact that 75 % of the catchments meet the weaker threshold and just over 62 % meet the stricter threshold.

These numbers decrease slightly for the behavioural ensemble as expected due to trade-offs between the four metrics, but the overall trends remain the same.

3.4.3 Spatial model performance

To analyse model performance spatially across GB, the four evaluation metrics for the best simulation (as defined by the combined rank across all four metrics) for each catchment are summarized in Fig. 7.

Figure 7 Model performance for the best simulation (as defined by the combined rank across all four metrics) for each evaluation metric (a) NSE (–), (b) bias in runoff ratio (%), (c) bias in low flow volume (%), and (d) bias in slope of the flow duration curve between the 30th and 70th percentile (%).

For NSE, model performance is variable across the country, but generally, better model performance is found in the wetter catchments in the north and west of GB, with poorer model performance in drier catchments in the south and east. Model performance is poor in groundwater-dominated areas, particularly in the underlying chalk regions in the south-east. This region has particularly low runoff coefficients (see Fig. 4d) and does not maintain mass balance with large water losses. Consequently, results for RRBIAS show that the model tends to overestimate flows in the south-east. While bias in the runoff ratio shows that the model is generally overestimating flows, biases in the low flow volume is a more mixed picture, with the model underestimating low flows in some locations, particularly in the Midlands and north-eastern Scotland. From Fig. 4d, these areas are characterized by particularly low flow duration curve slopes, suggesting strongly damped flow responses with high base flow. Flow in the Midlands region is heavily regulated by reservoirs which sustain low flows and could be a potential reason for overestimating low flows in this area. The bias in slope of the flow duration curve shows that DECIPHeR does well at reproducing the flow variability but tends to underestimate the slope in Scotland and northern Wales, suggesting that the hydrographs in these catchments are too smooth and not sufficiently flashy.

3.4.4 Relationship between model performance and catchment characteristics

To further analyse and understand the reasons for good or poor model performance, relationships between key catchment characteristics and model performance were further explored. Firstly, the catchments were grouped according to key catchment characteristics based on discharge, runoff coefficient and base flow index. The 5th, 50th and 95th percentiles of NSE and RRBIAS were calculated from the ensemble of runs for all catchments within each group to explore relationships between model performance and catchment characteristics (see Table 4). The relationship between runoff coefficient, wetness index and RRBIAS was also analysed to further explore the importance of water gains and/or losses on model performance.

Table 4 Summary statistics of DECIPHeR performance metrics for GB, with catchments grouped by runoff coefficient and base flow index. Percentiles are taken from the behavioural ensemble from all catchments within each group. The column N indicates the number of catchments in each group. Cells are formatted in typeface according to the thresholds outlined in Sect. 4.3.3: bold italic for the stricter threshold, bold for the weaker threshold and italic where it does not meet either of the thresholds.

There is a clear link between model performance and catchments with a low runoff coefficient. Table 4 highlights poor model performance in catchments where observed runoff coefficients are less than 0.2. In this group, the model always overpredicts (as shown by the RRBIAS result) and consequently leads to poor NSE scores. Figure 8 shows that for many catchments where the model overpredicts flows (and particularly for catchments with a runoff coefficient less than 0.2), observed potential evapotranspiration estimates are not high enough to account for water losses culminating in an overestimation of flows. This is unsurprising given that currently the model maintains water balance and cannot lose or gain water beyond the “natural” conceptualizations of precipitation, discharge and evaporation dynamics. Consequently, we are either missing a process (such as water loss due to inter-catchment groundwater flows or anthropogenic impacts) or the data are wrong.

Figure 8 Scatter plot of wetness index (mean annual precipitation divided by mean annual potential evapotranspiration), runoff coefficient (mean annual discharge divided by mean annual precipitation) and bias in runoff ratio for each GB catchment evaluated in this study. Any points above the horizontal dashed line are where runoff exceeds total rainfall inputs in a catchment, and any points below the curved line are where runoff deficits exceed total potential evapotranspiration in a catchment.

Poorer model performance is also found in groundwater-dominated catchments with a high base flow index (high BFI) catchments (Table 4); however, the results also show that we can also gain very good simulations in these types of catchments (5th percentile has a NSE score of 0.83). Hence, the challenge is to better understand water losses and/or gains in groundwater catchments to improve the representation of groundwater dynamics in the model.

3.4.5 Simulated flow time series

Finally, we examined the simulated flow time series for six example catchments with different characteristics. Figure 9 shows the observed discharge, observed precipitation and the 5th–95th percentile uncertainty bounds of the behavioural simulations for six catchments with different characteristics (see Table 5) for a representative 2-year period of the 55-year time series simulated. The 5th–95th percentile uncertainty bounds are generated from the likelihood-weighted distribution of the top 1 % of the model simulations using the GLUE framework (Beven, 2006).

Figure 9 Observed discharge and uncertainty bounds for the behavioural simulations (5th and 95th percentile of the likelihood-weighted simulated discharge) for six catchments with different characteristics (shown in Table 5). The plots show a 2-year period (2010–2012) from the 55-year time series simulated.

Table 5 Catchment characteristics and model performance for the six catchments shown in Fig. 9. Base flow index is a measure of the proportion of the river runoff that can be classified as base flow and is derived from Marsh and Hannaford (2008). Water balance is calculated as mean annual rainfall minus mean annual discharge and potential evapotranspiration (as actual evapotranspiration is not available). NSE and RRBias for the best ranked simulation according to the combined score described in Sect. 3.3.3 are shown for each catchment alongside the NSE and RRBias derived from the mean of the behavioural ensemble.

Our results show that the model can capture a range of different hydrological dynamics from wetter catchments in the north-west (Fig. 9a) to drier catchments in the south-east (Fig. 9b). While model performance for groundwater catchments can be very good (Fig. 9c and Table 5), it also shows that we need to incorporate additional model capability to simulate the dynamics of groundwater-dominated catchments. Where we have a very low (for Great Britain) runoff coefficient, this is assumed to involve water losses into a more regional groundwater storage not expressed at the outlet and not yet represented in this version of the model (Fig. 9d). While the catchments shown in Fig. 9a–d are relatively un-impacted by human influences, the catchments shown in Fig. 9e and f are heavily impacted by human influences and highlight the challenge of simulating flows nationally across catchments with diverse hydrological behaviour.

4 Outlook and ongoing developments

4.1 National-scale model evaluation

This is the first study to comprehensively benchmark hydrological model performance across GB. We calculated four evaluation metrics for 10 000 model simulations for 1366 GB gauges to provide an initial benchmark of model performance. DECIPHeR generally performs well for the flow time series evaluated in this study, with better results in the west and north in wet catchments as compared to drier catchments in the south and east. This is a common finding for hydrological models, with many studies finding poor model performance and greatest water balance errors in drier catchments (Gosling and Arnell, 2011; McMillan et al., 2016; Newman et al., 2015; Pechlivanidis and Arheimer, 2015). These results are also reflected in other GB model evaluation studies. For example, Coxon et al. (2014) applied FUSE to 24 GB catchments and found the best model performances in wet catchments compared to dry, chalk catchments, Rudd et al. (2017) evaluated G2G for low flows across 61 GB catchments and found positive bias in low flow volumes in small catchments in the south-east of England, and Crooks et al. (2010) evaluated the Probability Distributed Moisture model (PDM) across 120 GB catchments and found poorer model performance in groundwater-dominated, drier catchments.

Poor model performance in these catchments is partially due to some of the metrics chosen in this study; for example, percent bias is most sensitive to small absolute biases in the driest catchments when compared to other metrics such as absolute bias. However, positive bias in the runoff ratio could be caused by a number of factors such as underestimation of potential evapotranspiration (there are other UK gridded potential evapotranspiration products which estimate much higher potential evapotranspiration), inter-catchment groundwater flows and/or human influences such as water abstraction. Population density is much higher in the south and east compared to the north and west, so this regional disparity in model performance could also be explained by a greater rate of abstractions and managed watercourses which alter the flow time series. For example, 55 % of the effective rainfall in the Thames catchment is licensed for abstraction (Thames Water, 2017).

These results provide an initial test of DECIPHeR capabilities against a large sample of catchments, but this is only a first-order evaluation of model performance. A more rigorous evaluation would assess the model over different seasons (Freer et al., 2004), under changing climatic conditions (Fowler et al., 2016), for different hydrological extremes (Coron et al., 2012; Veldkamp et al., 2018; Zaherpour et al., 2018) and for multiple objectives simultaneously (Kollat et al., 2012) and incorporate input and flow data uncertainty (Coxon et al., 2014; Kavetski et al., 2006; McMillan et al., 2010; Westerberg et al., 2016).

4.2 Characterizing spatial heterogeneity and connectivity

The intended use of DECIPHeR is to determine how much spatial variability and complexity is required for a given set of modelling objectives. It can be run as a lumped model (one HRU), semi-distributed (multiple HRUs) or fully gridded (HRU for every single grid cell). In this paper DECIPHeR was applied across 1366 GB gauges, with catchment masks, 5 km input grids and three classes of accumulated area and slope as classifiers for the hydrological response units, resulting in a total of 133 286 HRUs. Future work needs to consider the appropriate spatial complexity and hydrologic connectivity needed to represent relevant processes (Andréassian et al., 2004; Blöschl and Sivapalan, 1995; Boyle et al., 2001; Chaney et al., 2016; Clark et al., 2015; Metcalfe et al., 2015; Wood et al., 1988). While this work highlights the clear potential of a computationally efficient large-scale modelling framework that can run large ensembles, a balance is required to ensure computational efficiency when running large ensembles that also maintains sufficient spatial complexity to represent different hydrological processes.

4.3 Hypothesis testing and model parameterization

To demonstrate the modelling framework, we implemented a single model structure, provided in the open-source model code, in all HRUs across GB and did not experiment with different model structures in different parts of the landscape. This provides a good benchmark of DECIPHeR's ability at the national scale across GB, but the results suggest that different model structures are needed to represent a greater heterogeneity of hydrological responses beyond the conceptual dynamics currently implemented in this simple model (as shown in Fig. 9). We can gain new process understanding of regional differences in catchment behaviour by testing different model representations (Atkinson et al., 2002; Bai et al., 2009; Perrin et al., 2001). Future work will concentrate on adding modules to DECIPHeR to enhance performance across national and continental scales, with a focus on improved representation of groundwater dynamic and human influences to address poor model performance in catchments with a low runoff coefficient. Furthermore, we have ensured that the code is open-source and well-documented so that the hydrological community can contribute new and/or different conceptualizations of the processes shown in this paper.

It is challenging to parameterize a hydrological model across large scales. Here we simply applied the same parameter set across each catchment. Using this basin-by-basin approach has the disadvantage of producing a “patchwork quilt” of parameter fields, with discontinuities in parameter values across catchment boundaries. This is only effective for gauged catchments (Archfield et al., 2015). Ongoing work aims to address these issues by implementing the multiscale parameter regionalization (MPR) technique for DECIPHeR across GB. This technique links model parameters to geophysical catchment attributes through transfer functions applied at the finest possible resolution (Samaniego et al., 2010). The coefficients of the transfer functions are then calibrated, and parameters are upscaled to produce spatially consistent fields of model parameters at any resolution across the entire model domain. The MPR technique has been applied elsewhere, proving that it can produce seamless parameter fields across large domains and produce scale-invariant parameters (Kumar et al., 2013; Mizukami et al., 2017; Samaniego et al., 2017), which is ideal for a flexible framework such as DECIPHeR.

5 Conclusions

DECIPHeR is a new flexible modelling framework which can be applied from small catchments to the continental scale for complex river basins resolving small-scale spatial heterogeneity and connectivity. The model is underpinned by a flexible, computationally efficient framework with a number of novel features.

1. Spatial variability and connectivity. It has the ability to modify spatial variability and connectivity in the model via the specification of hydrological response units with different topographic, landscape and input layers.

2. Model structures and parameterizations. It has the ability to experiment with different model structures and parameterizations in different parts of the landscape.

3. Computationally efficient. Grouping of hydrologically similar points in the landscape into hydrological response units enables faster run times.

4. Automated build. This allows easy application over large scales.

5. Open-source. The open-source model code is implemented in Fortran, with a user manual to help researchers and/or practitioners to use the model.

This paper describes the modelling framework and its key components and demonstrates the model's ability to be applied a large model domain. DECIPHeR is shown to be computationally efficient and perform well over large samples of gauges. This work highlights the potential for catchment- to continental-scale predictions by making use of available big datasets, advances in flexible modelling frameworks and computing power.

Appendix A: Analytical solution for kinematic subsurface flow

A1 Introduction

This appendix provides an analytical solution to the equations which were solved numerically by Beven and Freer (2001) in Dynamic TOPMODEL to route subsurface flows. Here, we use calculus to integrate the relevant equations through time, as opposed to the finite volume scheme described by Beven and Freer (2001), for better computational speed and increased numerical stability.

This development starts from the kinematic wave description of flow (a partial differential equation) and integrates that partial differential equation along the flow direction to obtain an ordinary differential equation in time. We then integrate that ordinary differential equation in time to get an analytical solution which gives the flow and storage at the end of each time step, in relation to the conditions at the start of the time step, and the inflow from both upslope and from drainage. Each hydrological response unit (HRU) in the model may be comprised of one or more sets of spatially contiguous cells. We use the term “spatial element” (SE) to refer to one of these contiguous sets of cells within an HRU. This is the same scale as referred to by Beven and Freer (2001) as a “group of elements”.

By first integrating the kinematic wave equation in space, we have effectively chosen to model the flow at this scale using a non-linear reservoir so that there is no wave travelling in space within a spatial element. A wave-like behaviour at larger scales is mimicked by having the groups of elements in a type of cascade (linked by the weighting matrix). This approach of integrating in space is the same as that selected by Beven and Freer (2001), using their finite volume approach.

A2 Flow in a spatial element

Assume the that spatial element (SE) has area A and that x is the distance measured along the flow direction of the SE. Define Q as the downslope flow rate (L³ T⁻¹) at some point x. Assume that the flow is kinematic, i.e. that Q depends only on S (L), the local storage deficit per unit area and the SE geometry. The drainage input from above is assumed to be r (L T⁻¹). Assume the width of the SE is w(x), at distance x (L). At any point x in the SE we can write a partial differential equation for Q:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial Sw}{\partial t}}={\displaystyle \frac{\partial Q}{\partial x}}-rw.\end{array}$$

This is a kinematic wave equation describing the subsurface flow at point x within an SE. Note that both S and r have been multiplied by w, the width of the SE, so that they can be compared with Q, which is the total flow through the SE, at distance x.

To simplify the problem, we will now average over the entire SE, along the flow direction, x, from the upslope end (x=0) to the downslope end (x=L). This will produce an equation describing how $\stackrel{\mathrm{‾}}{S}$, the SE average of S, changes with time:

$$\begin{array}{}\text{(A2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{L}}{\int }_{\mathrm{0}}^L{\displaystyle \frac{\partial Sw}{\partial t}}\mathrm{d}x={\displaystyle \frac{\mathrm{1}}{L}}{\int }_{\mathrm{0}}^L{\displaystyle \frac{\partial Q}{\partial x}}\mathrm{d}x-{\displaystyle \frac{\mathrm{1}}{L}}{\int }_{\mathrm{0}}^{L}rw\mathrm{d}x,\text{(A3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{L}}\underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{\partial Sw}{\partial t}}\mathrm{d}x={\displaystyle \frac{Q\left(L,t\right)-Q\left(\mathrm{0},t\right)}{L}}-r{\displaystyle \frac{\mathrm{1}}{L}}{\int }_{\mathrm{0}}^{L}w\mathrm{d}x.\end{array}$$

The variables Q(0, t) and Q(L, t) refer to flows at the upslope and downslope ends of the SE [L³ T⁻¹].

If we assume S and w are uncorrelated as x varies and let $W=\frac{\mathrm{1}}{L}{\int }_{\mathrm{0}}^{L}w\mathrm{d}x$, then

$$\begin{array}{}\text{(A4)}& {\displaystyle }W{\displaystyle \frac{\partial \frac{\mathrm{1}}{L}{\int }_{\mathrm{0}}^{L}S\mathrm{d}x}{\partial t}}={\displaystyle \frac{Q(L,t)-Q(\mathrm{0},t)}{L}}-rW.\end{array}$$

Dividing by W,

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}={\displaystyle \frac{Q(L,t)-Q(\mathrm{0},t)}{LW}}-r.\end{array}$$

Note that A=LW is the area of the SE, so we can now define $q=Q/A$ as flow per unit plan area (L T⁻¹), which is the same dimension as that used by Beven and Freer (2001):

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}=q(L,t)-\left(q\right(\mathrm{0},t)+r).\end{array}$$

In Eq. (6), q(0,t) and r are assumed to be known, and q(L,t), the outflow from the SE, is assumed to be a function of the mean deficit $\stackrel{\mathrm{‾}}{S}$. Thus the SE is being modelled as a non-linear reservoir, where $\stackrel{\mathrm{‾}}{S}$ is the state variable, the input is $q(\mathrm{0},t)+r$ and the outflow is $q(L,t)=f\left(\stackrel{\mathrm{‾}}{S}\right(t\left)\right)$. Note that the inflow is now assumed to be applied as a spatially uniform flux within the SE, rather than being applied at x=0. There is no representation of motion within the SE. Motion at larger scales is represented by the cascading of flow from one reservoir to another.

Note that in the following equations, Q is equivalent to Q_SAT (Eq. 4). Because no motion within the SE is represented, Q_IN and Q_UZ (Eq. 4) can be lumped together into a single term, here called r.

A2.1 Analytical solutions for an exponential conductivity profile

There are several parsimonious descriptions of the vertical profile of saturated hydraulic conductivity which are hydrologically plausible. Here we consider the standard exponential profile and a profile truncated at finite depth. In each case we find the analytical solution for both $\stackrel{\mathrm{‾}}{S}$ and q(L,t) as functions of time. Analytical solutions are also possible for the parabolic and linear profiles given in Ambroise et al. (1996).

Defining $u=\left(q\left(\mathrm{0},t\right)+r\right)$ and $q=q\left(L,t\right)$,

$$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle }q=q_{\mathrm{0}}\exp \left(-\stackrel{\mathrm{‾}}{S}/m\right),\text{(A8)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}=q-u.\end{array}$$

If we substitute Eq. (7) into Eq. (8), and integrate Eq. (8) from $\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)$ at t=0 up to $\stackrel{\mathrm{‾}}{S}\left(t\right)$, we obtain the intermediate result

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\exp \left(\stackrel{\mathrm{‾}}{S}/m\right)={\displaystyle \frac{q_{\mathrm{0}}}{u}}+\exp \left(\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)/m\right)\left(\mathrm{1}-{\displaystyle \frac{q_{\mathrm{0}}\exp \left(-\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)/m\right)}{u}}\right)\\ \text{(A9)}& {\displaystyle }& {\displaystyle }\exp \left(-{\displaystyle \frac{ut}{m}}\right).\end{array}$$

From this we can get expressions for both $\stackrel{\mathrm{‾}}{S}\left(t\right)$ and q(t):

$$\begin{array}{}\text{(A10)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{S}\left(t\right)=m\log \left[{\displaystyle \frac{q_{\mathrm{0}}}{u}}+\left({\displaystyle \frac{\mathrm{1}}{\exp \left(-\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)/m\right)}}-{\displaystyle \frac{q_{\mathrm{0}}}{u}}\right)\exp \left(-{\displaystyle \frac{ut}{m}}\right)\right],\text{(A11)}& {\displaystyle }& {\displaystyle }q\left(t\right)={\left[{\displaystyle \frac{\mathrm{1}}{u}}+\left({\displaystyle \frac{\mathrm{1}}{q\left(\mathrm{0}\right)}}-{\displaystyle \frac{\mathrm{1}}{u}}\right)\exp \left(-{\displaystyle \frac{ut}{m}}\right)\right]}^{-\mathrm{1}}.\end{array}$$

In the special case where u=0, we instead obtain

$$\begin{array}{}\text{(A12)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{S}\left(t\right)=m\log \left(\exp \left(\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)/m\right)+{\displaystyle \frac{q_{\mathrm{0}}t}{m}}\right),\text{(A13)}& {\displaystyle }& {\displaystyle }q\left(t\right)={\left[{\displaystyle \frac{\mathrm{1}}{q\left(\mathrm{0}\right)}}+{\displaystyle \frac{t}{m}}\right]}^{-\mathrm{1}}.\end{array}$$

Note that parameters m and q₀ in these equations are equivalent to SZM and Q₀ as defined in the paper.

A2.2 Exponential truncated smoothly at S_max (Beven and Freer, 2001; Eq. 9)

$$\begin{array}{}\text{(A14)}& {\displaystyle }& {\displaystyle }q=\left\{\begin{array}{ll}{q}_{\mathrm{0}}\cos \mathit{\beta }[\exp \left(-\cos \mathit{\beta }\stackrel{\mathrm{‾}}{S}/m\right)& \\ -\exp \left(-\cos \mathit{\beta }{S}_{max}/m\right)]& \stackrel{\mathrm{‾}}{S}\le S_{max}\\ \mathrm{0}& \stackrel{\mathrm{‾}}{S}>S_{max}\end{array}\right.,\text{(A15)}& {\displaystyle }& {\displaystyle }q=\left\{\begin{array}{ll}{q}_{\mathrm{1}}\exp \left(-\cos \mathit{\beta }\stackrel{\mathrm{‾}}{S}/m\right)-q_{\mathrm{2}}& \stackrel{\mathrm{‾}}{S}\le S_{max}\\ \mathrm{0}& \stackrel{\mathrm{‾}}{S}>S_{max}\end{array}\right.,\end{array}$$

where q₁=q₀cos β and $q_{\mathrm{2}}=q_{\mathrm{0}}\cos \mathit{\beta }\exp \left(-\cos \mathit{\beta }{S}_{max}/m\right)$.

Let's look first at the case where $\stackrel{\mathrm{‾}}{S}\le S_{max}$

$$\begin{array}{}\text{(A16)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}=q_{\mathrm{1}}\exp \left(-\stackrel{\mathrm{‾}}{S}/\left(m/\cos \mathit{\beta }\right)\right)-\left(q_{\mathrm{2}}+u\right).\end{array}$$

If we let m₂=m/cos β and $u_{\mathrm{2}}=q_{\mathrm{2}}+u$, then we can rewrite this as

$$\begin{array}{}\text{(A17)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}=q_{\mathrm{1}}\exp \left(-\stackrel{\mathrm{‾}}{S}/m_{\mathrm{2}}\right)-u_{\mathrm{2}}.\end{array}$$

This is now exactly the same form as the exponential profile above, so the solution is formally identical: we just used q₁ instead of q₀, m₂ instead of m and u₂ instead of u. The resulting equations are

$$\begin{array}{}\text{(A18)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{S}\left(\mathrm{t}\right)=m_{\mathrm{2}}log\left[{\displaystyle \frac{q_{\mathrm{1}}}{u_{\mathrm{2}}}}+\left({\displaystyle \frac{\mathrm{1}}{\exp \left(-\frac{\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)}{m_{\mathrm{2}}}\right)}}-{\displaystyle \frac{q_{\mathrm{1}}}{u_{\mathrm{2}}}}\right)\exp \left(-{\displaystyle \frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}}\right)\right],\text{(A19)}& {\displaystyle }& {\displaystyle }q\left(t\right)={\left[{\displaystyle \frac{\mathrm{1}}{u_{\mathrm{2}}}}+\left({\displaystyle \frac{\mathrm{1}}{q\left(\mathrm{0}\right)}}-{\displaystyle \frac{\mathrm{1}}{u_{\mathrm{2}}}}\right)\exp \left(-{\displaystyle \frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}}\right)\right]}^{-\mathrm{1}}.\end{array}$$

This solution collapses to the standard exponential result if cos β=1 and $S_{max}=\mathrm{\infty }$.

Note that provided $S_{max}<\mathrm{\infty }$, q₂>0, so u₂>0, and there is no need to consider the case of zero forcing.

The deficit cannot go beyond S_max as a result of outflow; however deficits larger than S_max can arise through evaporation (this prepares it for future developments; evaporation is not currently included in the conceptualization of the saturated zone). Here we consider the case where $\stackrel{\mathrm{‾}}{S}>S_{max}$, so q=0, but u>0, so the deficit is decreasing:

$$\begin{array}{}\text{(A20)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}}=-u.\end{array}$$

This can be integrated to give

$$\begin{array}{}\text{(A21)}& {\displaystyle }\stackrel{\mathrm{‾}}{S}\left(t\right)=\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)-ut.\end{array}$$

If $\stackrel{\mathrm{‾}}{S}\left(t\right)<S_{max}$, then we switch to the $\stackrel{\mathrm{‾}}{S}\le S_{max}$ solution partway through the computational interval. We use the equation for $\stackrel{\mathrm{‾}}{S}$ (t) when $\stackrel{\mathrm{‾}}{S}\le S_{max}$ because in that case the q(t) equation will lead to division by zero if it is started at q(0)=0.

A2.3 Extra note on computational issues

If u₂ is very small but not zero, numerical problems can arise in the calculation of $\stackrel{\mathrm{‾}}{S}\left(t\right)$ because of loss of significance when subtracting two numbers of very different magnitudes. This can lead to calculating the logarithm of zero during calculation of $\stackrel{\mathrm{‾}}{S}\left(t\right)$.

This can be avoided by making a Taylor series expansion of $\stackrel{\mathrm{‾}}{S}\left(t\right)$ for small non-zero values of u₂. We obtain

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{S}\left(t\right)\cong m_{\mathrm{2}}\log \\ \text{(A22)}& {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{q_{\mathrm{1}}}{u_{\mathrm{2}}}}+\left({\displaystyle \frac{\mathrm{1}}{\exp \left(-\frac{\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)}{m_{\mathrm{2}}}\right)}}-{\displaystyle \frac{q_{\mathrm{1}}}{u_{\mathrm{2}}}}\right)\left(\mathrm{1}-{\displaystyle \frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{!}}}{\left({\displaystyle \frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}}\right)}^{\mathrm{2}}\right)\right].\end{array}$$

If we expand and then neglect terms in $u_{\mathrm{2}}^{\mathrm{2}}$, we obtain

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{S}\left(t\right)\cong m_{\mathrm{2}}\log \\ \text{(A23)}& {\displaystyle }& {\displaystyle }\left[\left({\displaystyle \frac{\mathrm{1}}{\exp \left(-\frac{\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)}{m_{\mathrm{2}}}\right)}}+{\displaystyle \frac{q_{\mathrm{1}}t}{m_{\mathrm{2}}}}\right)-{\displaystyle \frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}}\left({\displaystyle \frac{\mathrm{1}}{\exp \left(-\frac{\stackrel{\mathrm{‾}}{S}\left(\mathrm{0}\right)}{m_{\mathrm{2}}}\right)}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{q_{\mathrm{1}}t}{m_{\mathrm{2}}}}\right)\right].\end{array}$$

We use this solution in cases where $\frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}\ll \mathrm{1}$, currently implemented as $\frac{u_{\mathrm{2}}t}{m_{\mathrm{2}}}<{\mathrm{10}}^{-\mathrm{10}}$.