In this paper, we introduce a new numerically robust distributed rainfall–runoff model for computationally efficient simulation at high spatio-temporal resolution: the distributed simple dynamical systems (dS2) model. The model is based on the simple dynamical systems approach as proposed by

Hydrological models are essential tools for applications ranging from sensitivity analysis to impact assessment and forecasting. Generally, the aim of rainfall–runoff models is to simulate streamflow given precipitation input. Depending on factors such as the research aim and the climatological/geological setting, different model structures or process representations might be preferred. This, in combination with the inherent complexity and heterogeneity of (sub)surface hydrological processes, has led to the development of numerous different hydrological models over the past decades, each with its own focus. Examples of such rainfall–runoff models and modelling tools include, amongst many others,
SWAT

Hydrological models are often classified from more conceptual to more process-based models. Conceptual models have fewer processes explicitly parametrized, and as a result their limited number of parameters makes them easier to calibrate. In conceptual models, catchments are often represented by a series of buckets or storages, which mimic processes with different response times. Their typical scale of application is that of small to medium (mesoscale) catchments, generally in a lumped fashion. Process-based models, on the other hand, contain many more parameters. These models are often applied in a distributed fashion, where many of the parameter values are based on maps of vegetation and soil properties. However, often conceptual parameters remain that require calibration or tweaking. Even though process-based models should give a better representation of the physical reality, many models can easily be beaten in performance by a simple neural network

Whereas conceptual models often perform satisfactorily at the daily resolution at the lumped basin scale, they lack the ability to explicitly simulate spatially distributed processes that might be needed to accurately predict streamflow response at larger scales. In a study over a large number of basins in France,

The need to improve spatially explicit information in hydrological models aligns with the increasing availability of high-resolution continental-scale forcing datasets. These include for instance merged radar data

Conceptual models have been applied in a (semi-)distributed manner to account for spatially distributed input data: a lumped model is applied at each individual grid cell, and water is most often transferred to the outlet using a post-processing function that accounts for catchment routing. This method of running the model subsequently for each individual grid cell is, however, not necessarily the most computationally efficient way to deal with spatially distributed data, but often the result of historical developments. Whereas increased computational power has driven the application of distributed models at increasingly fine (spatial) resolution

An example of a conceptual model that, in spite of its extreme simplicity, has shown good performance for discharge simulation at the scale of smaller catchments and at fine temporal (hourly) resolution is the simple dynamical systems (SDS) approach introduced by

Here, we present a flexible and computationally efficient distributed implementation and extension of the simple dynamical systems approach that can be used to investigate the spatially distributed hydrological response using data at high spatial and temporal resolution: the distributed simple dynamical systems (dS2) model. Our aim was not only to create a model that is capable of computationally efficient simulations of discharge in mesoscale basins at high spatio-temporal resolutions but also to develop a model code that is easy to use, read and modify. Therefore, the model is written in Python, although there are more computationally efficient languages available. We believe, however, that this is the optimum between model speed and code adjustability required in most hydrological model studies. In this model, the catchment is divided into smaller sections using a regular grid, and discharge is simulated according to the SDS approach for each pixel. This distributed implementation allows the concept to be applied to bigger catchments, and also allows for spatial heterogeneity, both in the forcing and in the parameters. Since the original concept consists essentially of only one differential equation, it can be applied in a computationally efficient fashion, vectorizing all cells in the catchment. Snow and routing modules are added to the model to allow for application in snow-dominated regions, and to transport the water from each pixel to the catchment outlet via the drainage network. This efficient distributed implementation lowers the computational burden for high-spatial- and/or high-temporal-resolution studies, and it opens doors for extensive uncertainty studies. We will first introduce the model concept and describe the technical application. After that, we discuss the parameter sensitivity, and finally we show an application of the model.

The simple dynamical systems approach proposed by

Originally, the SDS approach was introduced as a lumped approach to simulate discharge in small catchments of approximately 10 km

To respect the original scale of development and to capture spatial variability, we have developed a distributed implementation of the simple dynamical systems approach.
Our distributed implementation builds on the original concept as proposed by

Efficient distributed implementation of the simple dynamical systems approach, which is solved for each pixel. The left-hand side of panel

As previously mentioned, the sensitivity function is required to translate changes in storage to changes in discharge.

Parametrization of the discharge sensitivity, including the effects of the three parameters.

In our distributed implementation, we allow precipitation to fall as snow. Snow is treated as a separate storage (see Fig.

In Fig.

Simulation results with and without the snow conceptualization. Top panel shows the model input; bottom panel shows the model output. The four hydrographs represent runs with different snow parameter values (see legend, in mm h

A simple, efficient routing conceptualization was added to the distributed model to transport water from each grid cell to the outlet of the catchment. In most studies applying the simple dynamical systems approach, a constant delay factor was added to the generated runoff time series, in order to account for the delay caused by the river network. Since our distributed implementation has pixels at distinctly different locations, we need to account for attenuation caused by the river network.
The routing concept presented here is based on the already existing width function concept

The routing concept visualized:

The width function in Fig.

One of the main aims of this model was computational efficiency. However, we do not claim that we have built the most computationally efficient conceptual hydrological model, since we sacrificed some calculation time for the user-friendly and hugely popular Python programming language. By choosing this language, we encourage users of the model to adapt, improve and change all components of the model, in order to find answers to their research questions. The model is written in Python 3.6 and largely uses the NumPy library. NumPy uses C libraries to ensure fast computations over entire arrays, something the base functions of Python do not offer. Due to the simplicity of the simple dynamical systems approach, we need to solve only one equation (besides the snow conceptualization). NumPy allows functions to be vectorized: receiving and outputting an array of values, while applying the same function to each individual value. This is computationally more efficient than the step-by-step application of the same function to each element in the array.

To get an idea of the computational efficiency of the model, example runtimes are presented in Fig.

Runtimes for a simulation of 3 months with an hourly time step with varying number of pixels, run on a single core of a normal desktop (Intel Core i7-6700, 16GB RAM). The grey bars represent model runtimes where the numerical solver is called for each individual pixel and are keyed to the right-hand scale, making the difference in runtimes appear smaller than it really is. The dotted line represents the number of pixels when simulating Europe at 5

When applying the model to very large basins and/or running for very long periods, the random access memory (RAM) can become a limiting factor, as storing data in the RAM is by far most efficient. To prevent exceeding the RAM, the user can define the maximum amount of RAM the model is allowed to use, and dS2 will chunk the input and output data accordingly. For example, running the model for a basin with 1500 pixels for a period of 4 years on a hourly time step would require approximately 800 MB. Additionally, to reduce the time spent reading and writing the files on the disk, we rely on a special data format: NumPy memmap. This data file directly maps arrays to the hard drive without storing any metadata, to ensure fast reading and writing. The downside is that this format does not have any metadata, meaning the shape and data type need to be known prior to reading the file. In order to store the model output in a more common format, the model can transform the NumPy memmap data to the more commonly used NetCDF format, including metadata.

When solving a differential equation, one would ideally opt for an analytical solution. However, many equations (including Eq.

Since the simple dynamical systems approach allows for non-linear reservoirs, the simulated response can vary over several orders of magnitude within a single time step. This non-linearity can cause numerical errors, which can be prevented by reducing the time step.
Several options are available, where the Cash–Karp method is the typical textbook approach to explicitly solving differential equations

Three steps in the adaptive time stepping scheme. Note that in extreme cases both

The adaptive time stepping scheme operates in three steps. In the first step, the discharge for the next time step is calculated using the explicit Runge–Kutta 4 scheme. For both time steps, the discharge sensitivities are calculated based on the following equation:

During periods with low discharge and relatively high evaporation, the numerical solver can result in negative discharge amounts. To prevent this from happening, we define a threshold for discharge (which thus corresponds to a threshold storage level), below which no evaporation is allowed to occur. This mimics evaporation reduction during periods with low storage volumes, as discharge and storage are directly linked. For the pixels where the discharge is below this threshold, the evaporation rate is set to

The concept of this model is based on the water balance, yet it does not explicitly solve the water balance as most hydrological models do. The water balance is indirectly solved by calculating changes in storage. To check whether the model still respects the water balance, the closure of the water balance is calculated using the following equation:

Closure of the water balance per time step. Panel

To estimate the volume of water discharged during the time step, one needs a discharge value that is representative for the volume of water discharged during that time step. The easiest way to calculate this is to take the average of the discharges at time steps

We investigated the response of the model to the five main parameters (

Response surface plots for all main parameter combinations with the Kling–Gupta efficiency (KGE) as the performance metric for a synthetic experiment. The white dot in the middle of each graph represents the location where KGE is equal to 1. Each plot consists of 4900 model runs – where the model has 200 cells and was run for 1 year on an hourly time step – and took 1 h and 15 min to run on a normal desktop (Intel Core i7-6700, 16GB RAM).

Furthermore, we also analysed the parameter sensitivity according to

Sobol' parameter sensitivity for a synthetic experiment, shown for different performance statistics. The main effect presents the sensitivity induced by that parameter alone, and the total effect quantifies the sensitivity when parameter interactions are taken into account as well.
The bottom row

It is clearly visible that parameter sensitivity depends on the performance metric used. The NSE and KGE show roughly the same response, as is expected. The NSE calculated on the logarithmic discharge values shows a very different response, with

The most interesting patterns are visible in the bottom three subplots, where the three components of the KGE are presented. The correlation coefficient (

In this section, we apply the model to the Thur, a mesoscale basin (1700 km

For this application, we focus on three sub-basins in the Thur to validate the model at three different spatial extents (see Fig.

Application of dS2 to the Thur Basin in Switzerland. Panel

In Fig.

To elaborate this, we zoom in to a single event in Fig.

As stated in the Introduction, the aim of this paper is not to present the best or most comprehensive conceptual rainfall–runoff model, but to develop a model that can more easily be used to perform studies that require a large number of runs, such as uncertainty or sensitivity studies. As a result, this model concept has several limitations, either linked to the original simple dynamical systems approach or linked to the (distributed) implementation of this concept.

A limitation of the simple dynamical systems approach is that it assumes that the storage–discharge relation is the same on the rising limb and on the falling limb of the hydrograph. Studies have shown that hysteresis exists in multiple basins, especially those dominated by a variable contributing area

Surface runoff as a result of either infiltration excess or saturation excess is not explicitly accounted for. Surface runoff induced by oversaturated soils is implicitly accounted for in the discharge sensitivity function, which will reach high sensitivity values with high discharge (i.e. storage). During precipitation events with high intensity, surface runoff can also occur when the soil's infiltration capacity cannot accommodate the precipitation intensity. The model might therefore underestimate some discharge peaks resulting from infiltration excess overland flow. However, we expect that the influence of this process is minor in humid catchments. Not accounting for this process relates back to the purpose of this model, where we do not focus on correctly representing catchments' internal fluxes but instead focus on simulating the total discharge in a computationally efficient manner. It does, however, disqualify this concept for catchments with extreme precipitation events and the corresponding overland flow processes.

Furthermore, this model relies on the user to provide accurate actual evapotranspiration values as input for the model.
As can be seen in Eq. (

The vectorized implementation limits the ability to transfer water between neighbouring pixels. This subsurface lateral flow is mostly driven by gravity and is therefore expected to be most important in regions with large elevation differences. However, at the current recommended spatial–temporal scale of application (

As proposed by

Finally, the current technical implementation of the model implies two other limitations. First, the routing scheme only induces a time lag on the discharge peak, and it does not account for diffusion of the wave as it travels through the river network. For smaller basins, the effect of diffusion is only minor, but, when simulating larger basins such as the Rhine, diffusion of the discharge wave will play a more important role. We recommend using a more advanced routing scheme when applying the model to larger basins. The model can output a NetCDF file with runoff generated in each pixel per time step, which can easily be connected to other routing software such as SOBEK. Secondly, the current version of the model does not (yet) support multi-threading, meaning that the model currently only runs on a single thread. Spreading the computational load over multiple threads is expected to further reduce the time required to run the model, especially in large basins. To make the model more efficient in single-thread mode, one can split the entire basin into several smaller sections and simultaneously run the model for each section. After the simulation is complete, the user can combine the several resulting pixel-wise runoff maps and apply a routing scheme to the entire basin.

The distributed simple dynamical systems (dS2) model is a new rainfall–runoff model that aims to simulate discharge in a computationally efficient manner. This model is based on the simple dynamical systems approach proposed by

Synthetic examples demonstrate that, although dS2 solves the water balance implicitly, the model is able to accurately close the water balance. The response surface plots for all parameter combinations show that there are some correlations between the parameters, especially for the three discharge sensitivity parameters. However, we also showed that the parameters clearly influence different performance metrics, which provides the opportunity to reduce the calculation time of optimization algorithms. Finally, we have applied the model in the Thur Basin as a case study to validate the performance of the model. Our model is able to correctly simulate discharge, both at the local scale (e.g. the Rietholzbach catchment in Switzerland,

The dS2 model has several unique strengths compared to other rainfall–runoff models: (1) it is computationally efficient even at high temporal and spatial resolutions; (2) it only has five parameters to calibrate; (3) two parameters have a clearly defined influence on the discharge time series, making them easy to identify; and (4) the Python model code is open source and easily adjustable. The computational efficiency of this model creates the opportunity to answer different research questions. Since the model is able to simulate regions in relatively short amounts of time, performing sensitivity and uncertainty analyses with a large number of samples becomes feasible. Since the model is distributed, the sensitivity and uncertainty analyses can be performed both on a spatial and temporal basis. Furthermore, this model can be a valuable tool for educational purposes, to explain and directly show the effects of modifying parameters in distinct groups (e.g. hydrological response, snow, routing and evaporation). Overall, this makes dS2 a valuable addition to the already large pool of conceptual rainfall–runoff models, both for researchers and for practitioners interested in large sample studies.

The model code is open source and is archived at the 4TU repository:

The supplement related to this article is available online at:

JB, LAM and AJT conceived the idea behind dS2. JB developed the dS2 model, with valuable input from all co-authors. JB performed the synthetic experiments and case study; LAM and AJT helped with interpreting the results. JB prepared the manuscript, with significant contributions from all co-authors.

The authors declare that they have no conflict of interest.

We would like to thank MeteoSwiss and Massimiliano Zappa for providing the forcing data for the case study.

Funding for this research was provided by the Netherlands Directorate-General for Public Works and Water Management (Rijkswaterstaat).

This paper was edited by Richard Neale and reviewed by two anonymous referees.