the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Continuous statespace representation of a buckettype rainfallrunoff model: a case study with the GR4 model using statespace GR4 (version 1.0)
Guillaume Thirel
Charles Perrin
In many conceptual rainfall–runoff models, the water balance differential equations are not explicitly formulated. These differential equations are solved sequentially by splitting the equations into terms that can be solved analytically with a technique called “operator splitting”. As a result, only the solutions of the split equations are used to present the different models. This article provides a methodology to make the governing water balance equations of a buckettype rainfall–runoff model explicit and to solve them continuously. This is done by setting up a comprehensive statespace representation of the model. By representing it in this way, the operator splitting, which makes the structural analysis of the model more complex, could be removed. In this statespace representation, the lag functions (unit hydrographs), which are frequent in rainfall–runoff models and make the resolution of the representation difficult, are first replaced by a socalled “Nash cascade” and then solved with a robust numerical integration technique. To illustrate this methodology, the GR4J model is taken as an example. The substitution of the unit hydrographs with a Nash cascade, even if it modifies the model behaviour when solved using operator splitting, does not modify it when the statespace representation is solved using an implicit integration technique. Indeed, the flow time series simulated by the new representation of the model are very similar to those simulated by the classic model. The use of a robust numerical technique that approximates a continuoustime model also improves the lag parameter consistency across time steps and provides a more timeconsistent model with timeindependent parameters.
 Article
(834 KB)  Fulltext XML
 BibTeX
 EndNote
1.1 On the need for an adequate mathematical and computational hydrological model
Hydrological modelling is a widely used tool to manage rivers at the catchment scale. It is used to predict floods and droughts as well as groundwater recharge and water quality. In a review on the different existing hydrological models, Gupta et al. (2012) determined that all the existing models follow three modelling steps:

establish a conceptual representation of reality,

represent this conceptualization in a mathematical model,

set up a computational model to be used on a computer.
In terms of conceptual representation, many models exist and conceptualize the hydrological processes in the catchment differently, resulting in models with various levels of complexity. In this study, we will focus on the buckettype models, which are among the simplest. These models, such as Variable Infiltration Capacity (VIC) (Wood et al., 1992), Hydrologiska Bryåns Vattenbalansavdelning (HBV) (Bergström and Forsman, 1973) and Sacramento (Burnash, 1995), describe various conceptualizations of the hydrological processes at the catchment scale. Their parsimony (they usually need few input data and use few parameters) make them very useful for research as well as in operational applications thanks to their robustness and good performance (Michel et al., 2006).
In the context of this study, buckettype models are advantageous because, even if the concepts are often well documented, this is not the case of the mathematical and the computational models. In the models' documentations, the water balance equations that would govern the models are rarely explicitly formulated (Clark and Kavetski, 2010). The authors of the models often specify the discrete time equations, i.e. the result of the analytical or numerical temporal integration of the governing water balance equations. The problem is that the temporal resolution of the differential governing equations is part of the computational model. As a consequence, when the discrete time equations are the only ones available, the real mathematical model does not appear clearly. In addition, the descriptions of the numerical method used to solve the water balance equations and to obtain these discrete equations are rarely detailed.
However, several studies in the last decade (see for example Clark and Kavetski, 2010; Kavetski and Clark, 2010; Schoups et al., 2010) point out that the numerical solutions implemented to solve the differential equations that govern the models are sometimes poorly adapted. Clark and Kavetski (2010) showed that the use of the explicit Euler scheme (which is frequent for this type of model) can introduce significant errors in the simulated variables compared to more stable numerical schemes. Moreover, other studies prove that poorly adapted numerical treatment causes discontinuities and local optima in the parameter hyperspace (Kavetski et al., 2003; Kavetski and Kuczera, 2007; Schoups et al., 2010). This results in problems efficiently calibrating the models and in uncertainty on parameter values.
Another numerical approximation is commonly applied for buckettype models: the operator splitting (OS) technique (Kavetski et al., 2003). The aim is to split a differential equation into more simple equations that can be solved analytically in order to reduce inaccuracies in the numerical treatment. In the case of hydrological modelling, operator splitting results from the sequential calculation of processes such as runoff, evaporation and percolation (Schoups et al., 2010). Kavetski et al. (2003), Clark and Kavetski (2010) and Schoups et al. (2010) identified several widely used models in which the differential equations are solved using this type of treatment, e.g. VIC (Wood et al., 1992), Sacramento (Burnash, 1995) and GR4J (Perrin et al., 2003). However, even if OS may reduce numerical errors, Fenicia et al. (2011) cite several limitations to its use in hydrology. Indeed, it is physically unsatisfying to separate the different processes in time because, in reality, they are concomitant. In addition, it creates numerical splitting errors that are difficult to identify.
According to different studies, an inadequate numerical treatment like OS can lead to inconsistencies in flux simulations (see for example the study conducted by Michel et al., 2003, on an exponential store). It may also create inconsistencies in the model state variables (Clark and Kavetski, 2010; Kavetski and Clark, 2010). This results in the model inaccurately simulating flows.
For these reasons, it is important to use a robust numerical treatment to better estimate the other uncertainties (for example, parameter uncertainty).
1.2 Scope of this study
The first step to improve the numerical treatment of rainfall–runoff models is to properly separate the mathematical model from the computational model (Kavetski and Clark, 2010; Gupta et al., 2012). This article proposes a method to do this by setting up a continuous statespace representation of a rainfall–runoff model. A statespace representation is a matricial function of a system that depends on input, output and state variables. At all times, the system is described by the values of its state variables (referred to as “states” in this article). In the case of rainfall–runoff models, inputs can be potential evapotranspiration and precipitation and output can be the flow at the outlet of the catchment. The soil water content or the amount of water in the hydrographic network are physical examples of possible state variables. The level of the buckettype model stores is a conceptual example of possible state variables. This statespace representation will give the governing equations to be solved over time. This resolution will be proceeded by using an OS technique to be used as a comparison point and by using a more robust numerical technique, i.e. implicit Euler with an adaptive substep number. The model solved by implicit Euler will be called continuous statespace because it approximates a continuous model. By opposition, the OS statespace representation will be named as discrete.
In addition to a clearer mathematical model, we hope that the statespace representation will gain stability due to the direct implementation of the time step in the numerical resolution. We thus hope to obtain more stable parameter values across time steps (Young and Garnier, 2006).
To illustrate the methodology proposed, the widely used GR4J model (Perrin et al., 2003) will be taken as an example. Indeed, this model is currently implemented using the OS technique. A statespace representation will be set up, following the GR4J's conceptualization of the hydrological processes as much as possible. Its behaviour, both with a discrete and a continuous solving, will be compared to the current formulation of the GR4J model on a wide range of French catchments with different time steps (day and hour), in terms of performance and parameters.
Hereafter, the notation GR4 will be used to refer to the structure of the GR4J model (J stands for journalier, i.e. daily; Perrin et al., 2003), which is transformed and used at different time steps. This is a lumped buckettype model described in its current form (Sect. 2.1) and in its statespace form (Sect. 2.2). A discussion on the Nash cascade introduced in the GR4 statespace form is given in Sect. 2.3. The continuous differential equations of the statespace form are described in Sect. 2.4. The adaptations needed to change the model time step will be described in Sect. 2.5.
2.1 Reference GR4 model
GR4 (Perrin et al., 2003) is a lumped buckettype daily rainfall–runoff model with four free parameters. It is widely used for various hydrological applications in France (Grouillet et al., 2016; van Esse et al., 2013) and in other countries (Dakhlaoui et al., 2017; Seiller et al., 2017). It has shown good performances on a wide range of catchments (Coron et al., 2012). The equations of the reference GR4J model (Perrin et al., 2003) are the result of the integration of the water balance equations at a discrete time step (here the daily or hourly time step).
The version of GR4 used here is slightly different from the one presented by Perrin et al. (2003) because the two unit hydrographs were replaced by a single one placed before the flow separation (Fig. 1a, Mathevet, 2005). This simplification of the model does not substantially change the resulting simulated flows.
The equations of the model are given by Perrin et al. (2003) and listed in Table 1. GR4 represents the rainfall–runoff relationship at the catchment scale using an interception function, two stores, a unit hydrograph and an exchange function (see Fig. 1a). The model structure can be split into water balance and routing operators.
The water balance operators evaluate effective rainfall (i.e. the part of rainfall that will reach the catchment outlet) by estimating several quantities: actual evaporation, storage within the catchment and groundwater exchange. It involves an interception function and a production (soil moisture accounting) store (S in Fig. 1a). The interception corresponds to a neutralization of rainfall by potential evapotranspiration. The remaining rainfall (P_{n}), if any, is split into a part going into the production store (P_{s} in Fig. 1a) and a complementary part (P_{n}−P_{s} in Fig. 1a) that is directed to the routing component of the model. The quantity of rainfall that feeds the production store depends on the level of water in the store at the beginning of the time step. In case there is remaining energy for evapotranspiration after interception (E_{n} in Fig. 1a), some water is evaporated from the production store at an actual rate depending on the level of the production store (E_{s} in Fig. 1a). The higher the level is at the beginning of the time step, the closer E_{s} is to E_{n}. Thus, the production store represents the evolution of the catchment moisture content at each time step. The last water balance operator is a groundwater exchange term (F in Fig. 1a, positive or negative), which acts on the routing part of the model.
The routing function of the model is fed with the rainfall that does not feed the production store (P_{s}−P_{n}) plus a percolation term (Perc in Fig. 1a) from the production store, which generally represents a small amount of water. The total amount (P_{r} in Fig. 1a) is lagged by a symmetric unit hydrograph and then split into two flow components. The main component (90 % of P_{r}, Q_{9} in Fig. 1a) is routed by a nonlinear routing store (R in Fig. 1a). The complementary component (10 % of P_{r}, Q_{1} in Fig. 1a) directly reaches the outlet. The groundwater exchange term (F) is added or removed from the routing store and from the Q_{1} component.
The simulated flow at the catchment outlet (Q in Fig. 1a) is the sum of the outputs of the two flow components (Q_{r} and Q_{d} in Fig. 1a).
Four free parameters (called x_{1}, x_{2}, x_{3} and x_{4}) are used to adapt the model to the variety of catchments. Their meanings are given in Table 2.
As mentioned in the introduction, the governing water balance equations of the model are solved using OS. By considering that inputs to the store are added at the beginning of the time step as Dirac functions (Michel, 1991), it becomes possible to find analytical expressions of the model processes when equations are integrated over the time step. Consequently, the model processes are treated sequentially.
2.2 A statespace formulation for the GR4 model
To create this statespace representation, it is important to identify the different model state variables. In the GR4 model, two obvious states are the levels of the production and routing stores. The main challenge to describe the statespace formulation is to deal with the unit hydrograph. The discrete form used in GR4 corresponds to a convolution product in the state space as implemented in SUPERFLEX (Kavetski and Fenicia, 2011). This convolution product makes the mathematical resolution of the model that is necessary for the continuous version that will be introduced in Sect. 2.4 more complex. Here we chose to replace this unit hydrograph with a series of linear stores in order to simplify this resolution. The use of stores is also convenient because it creates a model that is only composed of stores.
Different combinations of linear stores were tested and the choice was made to replace the unit hydrograph with a Nash cascade (Nash, 1957). It is implemented at the same location in the model structure as the unit hydrograph (Fig. 1b). The Nash cascade is a chain of linear stores that empty into each other. It has two parameters to govern the shape of the outflow response, namely the number of stores and the outflow coefficient, which is identical for all stores. In our case, we decided to fix the number of stores and to only consider the outflow coefficient as a free parameter. This choice will be discussed in the following section (Sect. 2.3). With this type of model, the outflow of the last store has a similar shape to a unit hydrograph.
2.3 Parameterization of the Nash cascade
As introduced in the previous section, the Nash cascade has two parameters, namely the number of stores and the outflow coefficient. The number of stores can only take integer values, which is an issue for automatic calibration because it introduces threshold effects. As a consequence, the number of stores was not optimized automatically and the outflow coefficient is the preferential parameter to calibrate.
To obtain a response that is equivalent to the GR4 unit hydrograph response, we attempted to determine whether a relationship exists between the Nash cascade parameters and the GR4 x_{4} parameter. To manage this, the determination of the Nash cascade parameter is based on the comparison of the impulse response of the Nash cascade and the response of the unit hydrograph.
The impulse response of the Nash cascade is (Nash, 1957)
where h_{Nash}(t) is the impulse response of the Nash cascade at time t, nres is the number of stores, k is the outflow coefficient (t^{−1}) and Γ(nres) corresponds to the gamma function of nres.
The impulse response of the GR4 symmetrical unit hydrograph is (Perrin et al., 2003)
where h_{UH}(t) is the impulse response of the unit hydrograph at time t and x_{4} is the time to peak of the hydrograph.
The Nash cascade parameters are calculated depending on x_{4} in such a way that the time to peak and the peak flow would be the same for the two impulse responses. According to SzöllösiNagy (1982), the time to peak of the Nash cascade is equal to
and the peak flow is equal to
Using Eq. (2), the time to peak of the GR4 unit hydrograph is equal to
and the peak flow to
So, from these values the following system can be deduced:
which can be transformed into
A nres=11 is the best integer approximation to solve the second equation of Eq. (8). The outflow coefficient is deduced from this number of stores and from x_{4}. By fixing the parameters in this way, only the x_{4} parameter has to be calibrated. This method allows a direct comparison between the parameters of the Nash cascade and the parameter of the unit hydrograph. For a given x_{4} parameter, the unit hydrograph and the Nash cascade impulse responses have the same time to peak and the same peak flow (see the dotted and the dashed curve in Fig. 2).
Using this formula, the x_{4} parameters of the two models are equivalent and it can be argued that their meaning is nearly identical.
Fixing the number of stores in the Nash cascade also provides another advantage. Indeed, one of the potential issues that arises when replacing the unit hydrograph with a Nash cascade was the equifinality with the routing store. Given that the recession curve of the cascade is theoretically infinite, it could have the same function as the routing store. Calculating the parameters of the cascade regarding the x_{4} parameter makes it possible to reduce the possibility of an infinite impulse response.
2.4 Continuous differential equations of the statespace model
Once the model is only represented by stores, a differential equation can be written for each store (details are provided in Table 1). For the production and routing stores, the equations were built by adding all the processes that affect the stores. For example, the differential equation for the production store is the sum of the differential equations of evaporation, rainfall and the percolation (respectively, E_{s}, P_{s} and Perc in Fig. 1). This means that all the processes that are a function of this state are treated simultaneously, unlike the initial model version in which the processes are treated sequentially. The statespace representation of the Nash cascade is the same as the one proposed by SzöllösiNagy (1982).
The resulting model is composed of the differential equations governing the states' evolution (here represented as a vector in the Eq. (9), taking into account nres stores in the Nash cascade):
The notation $\dot{S}$ stands for $\frac{\mathrm{d}S}{\mathrm{d}t}$, the derivative of S against time t and the different elements of this equation are specified in Table 1.
The output equation to calculate the instantaneous output flow (q(t) in Eq. 10) completes the model:
The different elements in Eqs. (9) and (10) are shown in Table 1.
The input, state variable and output values are as follows.

Inputs. E_{n} and P_{n} are the potential evapotranspiration (after the interception) and the precipitation amounts after the interception phase (mm t^{−1}). We decided to keep the interception out of the statespace representation because it is not represented by a store in the reference GR4J and we wanted to avoid introducing an additional difference between the statespace and the reference models.

Output. Q is the output flow; it corresponds to the integration of q(t) (Eq. 10) over the time step.

State variables. S, R and S_{h,k} are respectively the levels of the production store, the routing store and the Nash cascade store number k (with $k\in \mathit{\{}\mathrm{1},\mathrm{\cdots},\text{nres}\mathit{\}}$) in millimetres.

Fluxes. P_{s} and E_{s} are, respectively, the rainfall added to the production store and the evapotranspiration extracted from the production store. Perc is the outflow from the production store. P_{r} is the amount of water that reaches the model routing operators. Q_{Sh,k} is the outflow of the Nash cascade store number k (with $k\in \mathit{\{}\mathrm{1},\mathrm{\cdots},\text{nres}\mathrm{1}\mathit{\}}$). Q_{uh} is the outflow of the Nash cascade store number nres (this notation is used to be consistent with the discrete model). Q_{9} and Q_{r} are, respectively, the inflow and the outflow of the routing store and F is the intercatchment groundwater exchange. Q_{d} is the outflow of the complementary flow component.
The parameter meanings are explained in Table 2. The model is constructed to ensure that the parameters (${x}_{\mathrm{1}},\mathrm{\cdots},{x}_{\mathrm{4}}$ in the equations) have the same meaning in the continuous model and in the discrete GR4. The statespace formulation was sought to be as close as possible to the original model's formulation, to keep the same general model structure. We expect similar results to be obtained by the different tested model versions.
2.5 Hourly model
The GR4 model was first designed for daily time step modelling and it was adapted for the hourly time step (GR4H, Mathevet, 2005; Ficchí et al., 2016). The structure and the equations are similar in GR4H (hourly) and in GR4J (daily). The hourly versions of the GR4 models used here are the same as the ones showed in Fig. 1.
The adaptation to the time step is handled by a change in the parameter values, which depend on time. Ficchí et al. (2016) gave the theoretical relationships to transform the GR4 free parameter values as a function of the time step length (Table 3). The fixed percolation coefficient (ν in Table 1) is also time dependent.
The continuous statespace GR4 model used for the hourly time step is exactly the same as the one used at the daily time step, with no change in the percolation coefficient. The time step change is not managed by a change in parameter values but by the numerical integration. For the daily time step, the model is integrated on Δt=1 day while; for the hourly time step, it is integrated on Δt=1 h.
3.1 Numerical integration of the model
The integration of Eq. (9) (necessary to adapt the model to discrete input data) cannot be made analytically. It is therefore necessary to implement a numerical method to solve this integration.
Following the recommendation in Clark and Kavetski (2010), an implicit Euler algorithm is used to perform this numerical integration. Our choice was to set up an adaptive substep algorithm (Press et al., 1992) to avoid the majority of numerical errors. The implicit equation is solved using a secant method when necessary.
The choice of using an adaptive substep rather than singlestep implicit method (as recommended by Clark and Kavetski, 2010) is a result of several tests that are not shown here. We compared the modelling results with a singlestep integration to those obtained with the adaptive substep algorithms and found some differences in resulting flows (in particular for high flows). The differences found this way were not negligible. In this case, we can say that the stability of the implicit singlestep integration is not sufficient to sufficiently reduce the integration errors.
For both hourly and daily time steps, the inputs are considered as constant during the time step. Even if this assumption is a simplification of the truth, we chose to keep it constant to simplify the calculation and not to introduce treatment differences between hourly and daily time step models.
3.2 Catchment set and data
To compare the performance and behaviour of the reference and the discrete and continuous statespace GR4 model versions, a large data set of 240 catchments across France was set up (Fig. 3). Testing the models on many catchments will help obtain general conclusions (Andréassian et al., 2006; Gupta et al., 2012).
The data set was built by Ficchí et al. (2016) to test GR4 at different time steps. In this article, we only used daily and hourly data. The climate data of the SAFRAN daily reanalysis (Quintana Seguì et al., 2008; Vidal et al., 2010) are used as input data (precipitation and temperature). Precipitation and temperature are spatially aggregated on each catchment since the GR4 models are lumped. The hourly precipitation data were obtained by disaggregating the daily SAFRAN precipitation using the subdaily distribution of rain gauge measurements. Potential evapotranspiration at the daily time step was calculated from the SAFRAN temperature using the Oudin formula (Oudin et al., 2005) and hourly spread with a Gaussian distribution. Full details on this data set are available in Ficchí et al. (2016).
Hourly observed flows are available at each catchment outlet and come from the Banque HYDRO (http://www.hydro.eaufrance.fr/, last access: 18 April 2018, French Ministry of the Environment). For daily modelling, hourly measurements are aggregated at the daily time step. Their availability covers the 2003–2013 period.
The catchments were selected to have less than 10 % precipitation falling as snow, to avoid requiring a snow model.
3.3 Testing methodology
Three versions of the model were assessed on the 240 catchments following a splitsample test (Klemeš, 1986). These three versions are the reference model, a discrete statespace model (with a Nash Cascade but solved using OS) and a continuous statespace model. Comparing the reference and discrete statespace models allows us to measure the impact of replacing the unit hydrograph with a Nash cascade. Comparing the discrete and continuous statespace models allows us to measure the impact of a nearly continuous numerical integration. For every catchment, the observed flow data period was divided into a calibration period (the first half) and a validation period (the second half). A 2year warmup period was used for each catchment, before both the calibration and validation periods. The calibration was made automatically with an algorithm used in Coron et al. (2017) and based on the work of Michel (1991).
The objective function used for calibration is the Kling–Gupta efficiency (KGE^{′}; Kling et al., 2012). This objective function is often used in hydrology and assesses different components of the error made by the model (mean bias, variance bias, correlation). In addition, to target different flow levels, mathematical transformations are applied (Pushpalatha et al., 2012). The logarithm is applied to analyse the errors in lowflow conditions (KGE^{′}(log (Q))); no transformation is applied to preferentially analyse the error on high flows (KGE^{′}(Q)) and the root square of the flow is used as a compromise representing the error on intermediate flows (${\text{KGE}}^{\prime}\left(\sqrt{Q}\right)$). In the case of logarithm transformation, following the recommendations made by Pushpalatha et al. (2012), a small quantity which corresponds to onehundredth of the catchment mean flow is added to avoid troubles with null flows. These three transformations represent three distinct objective functions. The models were calibrated separately and successively on the three objective functions. To avoid strongly negative values of the KGE^{′} criterion, we used the C_{2M} formulation, which restricts the variation range into $[\mathrm{1};\mathrm{1}]$ (see Mathevet et al., 2006).
The results of the calibrations were also analysed in terms of performance in validation on the three evaluation criteria (i.e. C_{2M}(Q), C_{2M}(log (Q)) and ${C}_{\mathrm{2}M}\left(\sqrt{Q}\right)$). Given the large number of catchments, it is possible to draw a conclusion on the global difference in performance among the three studied model versions. This avoids a discrepancy due to specific catchment conditions. In addition to the performance analysis, the simulated hydrographs were visually analysed to detect discrepancies in the flow simulation. An analysis of the time series of internal fluxes and state variables also provided further insights to interpret the difference among the model versions. Last, the differences in parameter values among the models was analysed. It is important to verify that the parameter values are similar and do not take outlier values that would compensate for model inconsistencies.
A second test was carried out in order to analyse the time step dependency of the models. The splitsample test was performed at the hourly time step and the parameter values were compared to those obtained at the daily time step. With the reference model, the calibrated parameter values were compared to those theoretically obtained using the equations in Table 3. With the continuous statespace model, we verified the stability of the parameters. This stability is very important for designing a model that is not dependent on its time step.
4.1 Comparison of tested models at the daily time step
Figure 4 shows that performances are globally similar among the different versions of the model with a calibration using the C_{2M} on squarerooted flows. The performances of the reference model and the continuous statespace solution are also similar after calibration with the two other transformations of the flow in the objective function (not shown). In the case of the discrete statespace solution, the model does not seem to be able to reproduce high flows well but performs better on low flows than the two other models when the objective function used is the C_{2M} with logarithmic transformation.
The study of the hydrographs provides complementary information. The reference GR4 model and the continuous statespace solution are very similar while the discrete statespace solution simulates lower peak flows (see example hydrograph in Fig. 5). This behaviour can be explained because solving the 11 linear stores introduces errors that propagate and amplify across the Nash cascade.
To extend the analysis on the similarity of the models, we compared the parameter values obtained by calibration. As shown in Fig. 6, the parameters have the same range of values. We can still note differences in the values of the x_{4} parameter, which are systematically higher for the discrete statespace model. These differences in the values are probably due to the differences in response shape between the Nash cascade and the unit hydrograph (see Sect. 2.3) and to the errors produced by OS solving of the Nash cascade. The assumption that the differences in x_{4} values are due to errors caused by unsuitable solving is confirmed by the fact that the x_{4} parameter values are similar for the three models at an hourly time step (not shown here).
Last, to understand the internal impact of the statespace formulation on the model, we analysed state variables and internal fluxes. Two differences are induced by the model's statespace formulation. First, the discrete Nash cascade output peaks are lower than the peaks of the unit hydrograph (Fig. 7). The peaks of the continuous statespace representation are more similar with the reference but the peaks occur sooner. The second difference between the models concerns the levels of the routing store (Fig. 8). Here we only compared the reference GR4 to the continuous statespace solution because the inputs in the routing store are too different for the discrete statespace solution. The peak levels are higher in the continuous statespace representation, even sometimes higher than the maximum capacity of the routing store. The reason for this is that we shifted from the discrete model in which the processes are treated sequentially to a continuous model in which all the processes are solved simultaneously. In the discrete model, the exchanges are first calculated based on the routing level at the beginning of the time step, then the output of the unit hydrograph is added and last the outflow of the routing store is calculated. Due to this sequential treatment, in highflow conditions, the quantity of exchanged water and the outflow of the routing store in the discrete model are lower than those of the continuous statespace representation. Given that most of the time the exchange parameter is negative, the lower outflow of the routing store is compensated for by less water loss with the groundwater exchange in the complementary flow branch. This can explain why the simulated flows are similar despite these internal differences.
Moreover, by analysing the differences between the two models, it is also important to take into account the computational time. Indeed, running the original model version is on average 3 times faster than the continuous statespace version due to the adaptive substep method. This is important to consider for some applications.
This computational time rise is essentially due to the adaptive substep algorithm. For example, in the River Azergues at Châtillon catchment, the mean number of substeps is 22 and it can reach 100 during some days. However, in Sect. 3.1 we argue that the adaptive substep method seems necessary to avoid numerical errors.
To conclude with these results, we can argue that the modifications brought by the continuous statespace representation, although they modify the model's internal fluxes, do not degrade the model's performance, but only slightly modify the model's internal fluxes. It is important to underline that the OS solving of a Nash cascade creates more errors than a discrete unit hydrograph. To be equivalent to the reference model, the statespace representation of GR4 needs to be solved with a robust numerical technique.
4.2 Consistency of the statespace representation through time steps
The analysis of temporal consistency provides the most valuable result produced by the continuous statespace representation. The work of Ficchí et al. (2016) resulted in a GR4 model that is nearly consistent across time steps. However, to adapt the model, they chose to include the time step variations in a theoretical transformation between the free parameter values and the percolation fixed coefficient (Table 3) at different time steps. In this section, we only compare the reference GR4 with the continuous solution of the statespace representation. The parameters of the statespace representation discrete solution show the same behaviour as the reference GR4 ones so it was chosen not to show them. This proves that all the improvements shown in this section are only due to the continuous resolution of the statespace model.
In Fig. 9, the free parameter values obtained by calibration at the hourly time step are compared to those obtained at the daily time step using the reference GR4 version. The dashed lines represent the regression obtained by the theoretical relations reported in Table 3. One can note that the calibrated parameters (the dots in Fig. 9) are quite different between the two time steps but it is important to note that the values of the x_{3} parameter follow the relations proposed by Ficchí et al. (2016) (the dashed lines). The high values of x_{1} are underestimated compared to the theoretical relation as are the low values of the x_{2} parameter. There is also an issue with the unit hydrograph parameter (x_{4} in Fig. 9) for which calibrated hourly parameter values are systematically lower than the values it would have by following the transformation. Kavetski et al. (2011) and Littlewood and Croke (2008) encountered the same issue with the lag parameter of their models.
The values of x_{1}, x_{2} and x_{4} are inconsistent compared to the values expected using the theoretical transformations. Regarding the work of Ficchí (2017), we can argue that the changes in the high values of x_{1} and the low values of x_{2} are due to temporal inconsistencies in the interception calculation. The case of the x_{4} parameter is more problematic. The differences in the x_{4} values probably stem from the discretization of the unit hydrograph at different time steps.
In the continuous statespace model, the time step is taken into account in the temporal numerical integration of the model. For this reason, in theory there is no need to adapt the values of the parameters. This is confirmed in Fig. 10, where the values of calibrated parameters remain approximately constant despite the time step change. Only the high values of x_{1} and the values of x_{2} slightly diverge from the x=y line.
This result is useful in building a model that can adapt its time step resolution depending on the given conditions. The results are particularly interesting for the case of x_{4} values because the x_{4} values are constant between the two time steps, resolving the issue encountered by Littlewood and Croke (2008), Kavetski et al. (2011) and Ficchí et al. (2016) with lag parameters. As explained in the work of Littlewood and Croke (2013), this improvement can be explained by the fact that the adaptive substep integration approximates a continuous time input in the Nash cascade. The results obtained with the x_{4} parameter here tend to confirm this earlier work on a wide range of catchments. However, in addition to the input errors, the lack of x_{4} time consistency can also be explained by the integration errors produced by the OS at a daily time step.
The outliers in x_{3} values that occur in Fig. 10 are also present in Fig. 9. No explanations relating to physical characteristics of these catchments or simulation performance were found. We assume that these outlier values are due to the nonsensitivity of the x_{3} parameter for these catchments.
Finally, to verify stability, we also need to compare the performance of the models at the hourly time step. Figure 11 shows that, as at the daily time step, the performance is similar for the different versions.
Thus, the continuous statespace representation shows better temporal stability in the x_{4} parameter values with similar performance.
The objective of this study was to present a version of a buckettype rainfall–runoff model with a robust numerical resolution of the governing water balance equations by setting up a continuous statespace representation. The methodology is based on (i) identifying the state variables, (ii) writing their differential equations, (iii) replacing certain components of the model with more easily described components in terms of differential equations (namely replacing the unit hydrograph with a Nash cascade here), and (iv) solving these equations with a robust numerical integration technique. Finally, all the fluxes that form the water balance equation governing a state are solved simultaneously while they are solved sequentially in OS models. As stated by Fenicia et al. (2011), this is more physically satisfying.
This work was presented using the example of the GR4 model. The new version was created to be as close as possible to the initial model but a single modification was implemented: a Nash cascade substitutes the model's unit hydrograph.
When analysing the results and the output flows, it was shown that the new formulation, when solved with a robust numerical technique, has a limited impact on performance. However, the analysis of the parameter values and of the internal fluxes of the model shows that some discrepancies occur when running the model. The peak flow of the Nash cascade occurs sooner than the peak flow of the unit hydrograph. The amount of water in the routing store and exchanged by the groundwater exchange function is also higher for the statespace representation, particularly during highflow periods.
Nonetheless, the continuous statespace representation simulates flows that are very similar to the flows simulated by the original GR4 version and performs equally well. It also seems to provide greater stability in the parameter values, particularly regarding different modelling time steps. Moreover, the use of the Nash cascade rather than the unit hydrograph improves (when solved with implicit Euler) the lag parameter value stability with time steps. This improved stability can make it easier to calibrate the model with a given data set and to apply it at a finer time step for which no discharge data are available. It can also allow the use of a model that runs at a finer time step in highflow periods and a larger time step in lowflow periods.
Furthermore, the comparison between the discrete and continuous statespace model shows that the benefits provided by the continuous statespace representation are a result of the use of a robust numerical integration technique. Indeed, solving the statespace representation using OS introduces errors that impact the simulated flow values and do not result in parameter stability. Thus, the real benefit of the use of the Nash cascade is to simplify the numerical solving application.
The performance obtained with the continuous statespace model is not better than that of the original model. In addition, because the number of substeps sometimes needs to be high, the computational time is longer with the continuous statespace representation of the model. Consequently, the use of this representation would be helpful for particular applications such as timevariable modelling. It might also be useful for certain data assimilation techniques (typically variational methods) because all the components are represented as states and the governing equations are clearly defined.
In addition, it could also be advantageous to find a way to adapt the number of stores of the Nash cascade to the catchment studied.
Although it is necessary to adapt the Nash cascade to different unit hydrograph shapes, this article suggests a sufficiently general methodology to erase OS in hydrological buckettype modelling and can be transposed to other models.
The Fortran code used in this article can be freely downloaded from GitHub at https://github.com/HYDROgroupIrsteaAntony/GR4Statespaceversion1.0 (last access: 18 April 2018). The statespace model can be tested on an example catchment data set with already calibrated model parameters. The full reference for this code can be found in the references (Santos, 2017); it is referenced with the following DOI: https://doi.org/10.5281/zenodo.1118183.
This work is part of LS' PhD work; he made the technical development and the analysis and wrote the paper. GT and CP are the PhD supervisors; they supervised this work and the paper writing.
The authors declare that they have no conflicts of interest.
The first author's PhD grant was provided by Irstea. We thank Météo France for providing the SAFRAN climatic data used in this work. We also would like to thank Martyn Clark for his advice in settingup the differential equations, Nicolas Le Moine for sharing his ideas to replace the unit hydrographs and on the numerical integration, Fabrizio Fenicia for his advice on numerical integration and PaulHenry Cournède for his analysis of the mathematical adequacy of the model. Finally, we give special thanks to Andrea Ficchí for his work on the database and for the discussions on the temporal stability of the GR4 model.
We thank the topical editor, Jeffrey Neal, for his monitoring of the review
process and his relevant reviewers choice. We also acknowledge the two
reviewers, Barry Croke, and the anonymous reviewer for their very interesting
and complementary remarks.
Edited by: Jeffrey Neal
Reviewed by: Barry Croke and one anonymous referee
Andréassian, V., Hall, A., Chahinian, N., and Schaake, J.: Large sample basin experiments for hydrological model parametrization, chap. Introduction and Synthesis: Why should hydrologists work on a large number of basin data sets?, IAHSAISH P., 307, 1–5, 2006. a
Bergström, S. and Forsman, A.: Development of a conceptual and deterministic rainfallrunoff model, Nord. Hydrol., 4, 147–170, 1973. a
Burnash, R. J. C.: The NWS river forecast system – catchment modeling, chap. 10, in: Computer Model of Watershed Hydrology, Water Resources Publications, Highlands Ranch, Colorado, USA, 311–366, 1995. a, b
Clark, M. P. and Kavetski, D.: Ancient numerical daemons of conceptual hydrological modeling: 1. Fidelity and efficiency of time stepping schemes, Water Resour. Res., 46, W10510, https://doi.org/10.1029/2009wr008894, 2010. a, b, c, d, e, f, g
Coron, L., Andréassian, V., Perrin, C., Lerat, J., Vaze, J., Bourqui, M., and Hendrickx, F.: Crash testing hydrological models in contrasted climate conditions: An experiment on 216 Australian catchments, Water Resour. Res., 48, W05552, https://doi.org/10.1029/2011WR011721, 2012. a
Coron, L., Thirel, G., Delaigue, O., Perrin, C., and Andréassian, V.: The suite of lumped GR hydrological models in an R package, Environ. Modell. Softw., 94, 166–177, https://doi.org/10.1016/j.envsoft.2017.05.002, 2017. a
Dakhlaoui, H., Ruelland, D., Tramblay, Y., and Bargaoui, Z.: Evaluating the robustness of conceptual rainfallrunoff models under climate variability in northern Tunisia, J. Hydrol., 550, 201–217, https://doi.org/10.1016/j.jhydrol.2017.04.032, 2017. a
Fenicia, F., Kavetski, D., and Savenije, H. H. G.: Elements of a flexible approach for conceptual hydrological modeling: 1. Motivation and theoretical development, Water Resour. Res., 47, W11510, https://doi.org/10.1029/2010wr010174, 2011. a, b
Ficchí, A.: An adaptive hydrological model for multiple timesteps: Diagnostics and improvements based on fluxes consistency, PhD thesis, Université Pierre et Marie Curie, Université Pierre et Marie Curie, Paris, France, 2017. a
Ficchí, A., Perrin, C., and Andréassian, V.: Impact of temporal resolution of inputs on hydrological model performance: An analysis based on 2400 flood events, J. Hydrol., 538, 454–470, https://doi.org/10.1016/j.jhydrol.2016.04.016, 2016. a, b, c, d, e, f, g, h
Grouillet, B., Ruelland, D., Vaittinada Ayar, P., and Vrac, M.: Sensitivity analysis of runoff modeling to statistical downscaling models in the western Mediterranean, Hydrol. Earth Syst. Sci., 20, 1031–1047, https://doi.org/10.5194/hess2010312016, 2016. a
Gupta, H. V., Clark, M. P., Vrugt, J. A., Abramowitz, G., and Ye, M.: Towards a comprehensive assessment of model structural adequacy, Water Resour. Res., 48, W08301, https://doi.org/10.1029/2011wr011044, 2012. a, b, c
Kavetski, D. and Clark, M. P.: Numerical troubles in conceptual hydrology: Approximations, absurdities and impact on hypothesis testing, Hydrol. Process., 25, 661–670, https://doi.org/10.1002/hyp.7899, 2010. a, b, c
Kavetski, D. and Fenicia, F.: Elements of a flexible approach for conceptual hydrological modeling: 2. Application and experimental insights, Water Resour. Res., 47, W11511, https://doi.org/10.1029/2011wr010748, 2011. a
Kavetski, D. and Kuczera, G.: Model smoothing strategies to remove microscale discontinuities and spurious secondary optima in objective functions in hydrological calibration, Water Resour. Res., 43, W03411, https://doi.org/10.1029/2006wr005195, 2007. a
Kavetski, D., Kuczera, G., and Franks, S. W.: Semidistributed hydrological modeling: A “saturation path” perspective on TOPMODEL and VIC, Water Resour. Res., 39, 1246, https://doi.org/10.1029/2003wr002122, 2003. a, b, c
Kavetski, D., Fenicia, F., and Clark, M. P.: Impact of temporal data resolution on parameter inference and model identification in conceptual hydrological modeling: Insights from an experimental catchment, Water Resour. Res., 47, W05501, https://doi.org/10.1029/2010wr009525, 2011. a, b
Klemeš, V.: Operational testing of hydrological simulation models, Hydrolog. Sci. J., 31, 13–24, https://doi.org/10.1080/02626668609491024, 1986. a
Kling, H., Fuchs, M., and Paulin, M.: Runoff conditions in the upper Danube basin under ensemble of climate change scenarios, J. Hydrol., 424–425, 264–277, https://doi.org/10.1016/j.jhydrol.2012.01.011, 2012. a
Littlewood, I. G. and Croke, B. F. W.: Data timestep dependency of conceptual rainfall streamflow model parameters: an empirical study with implications for regionalisation, Hydrolog. Sci. J., 53, 685–695, https://doi.org/10.1623/hysj.53.4.685, 2008. a, b
Littlewood, I. G. and Croke, B. F. W.: Effects of data timestep on the accuracy of calibrated rainfall–streamflow model parameters: practical aspects of uncertainty reduction, Hydrol. Res., 44, 430–440, https://doi.org/10.2166/nh.2012.099, 2013. a
Mathevet, T.: Quels modèles pluiedébit globaux au pas de temps horaire? Développements empiriques et comparaison de modèles sur un large échantillon de bassins versants, PhD thesis, Ecole Nationale du Génie Rural, des Eaux et des Forêts, Paris, France, 2005 (in French). a, b
Mathevet, T., Michel, C., Andréassian, V., and Perrin, C.: A bounded version of the NashSutcliffe criterion for better model assessment on large sets of basins, IAHSAISH P., 307, 211–219, 2006. a
Michel, C.: Hydrologie appliquée aux petits bassins versants ruraux, Tech. rep., Cemagref, Antony, 320 pp., 1991 (in French). a, b
Michel, C., Perrin, C., and Andreassian, V.: The exponential store: a correct formulation for rainfall runoff modelling, Hydrolog. Sci. J., 48, 109–124, https://doi.org/10.1623/hysj.48.1.109.43484, 2003. a
Michel, C., Perrin, C., Andréassian, V., Oudin, L., and Mathevet, T.: Has basinscale modelling advanced beyond empiricism?, IAHSAISH P., 307, 108–116, 2006. a
Nash, J. E.: The form of the instantaneous unit hydrograph, Int. Assoc. Sci. Hydrol. Publ., 45, 114–121, 1957. a, b
Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V., Anctil, F., and Loumagne, C.: Which potential evapotranspiration input for a lumped rainfall–runoff model?, J. Hydrol., 303, 290–306, https://doi.org/10.1016/j.jhydrol.2004.08.026, 2005. a
Perrin, C., Michel, C., and Andréassian, V.: Improvement of a parsimonious model for streamflow simulation, J. Hydrol., 279, 275–289, https://doi.org/10.1016/s00221694(03)002257, 2003. a, b, c, d, e, f, g, h, i, j, k
Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical recipes in C, 2nd Edn., Press Syndicate of the University of Cambridge, 1992. a
Pushpalatha, R., Perrin, C., Moine, N. L., and Andréassian, V.: A review of efficiency criteria suitable for evaluating lowflow simulations, J. Hydrol., 420–421, 171–182, https://doi.org/10.1016/j.jhydrol.2011.11.055, 2012. a, b
Quintana Seguì, P., Le Moigne, P., Durand, Y., Martin, E., Habets, F., Baillon, M., Canellas, C., Franchisteguy, L., and Morel, S.: Analysis of NearSurface Atmospheric Variables: Validation of the SAFRAN Analysis over France, J. Appl. Meteorol. Clim., 47, 92–107, https://doi.org/10.1175/2007JAMC1636.1, 2008. a
Santos, L.: HYDROgroupIrsteaAntony/GR4Statespaceversion1.0: First release of GR4Statespaceversion1.0, Zenodo, available at: https://doi.org/10.5281/zenodo.1118183 (last access: 18 April 2018), 2017. a
Schoups, G., Vrugt, J. A., Fenicia, F., and van de Giesen, N. C.: Corruption of accuracy and efficiency of Markov chain Monte Carlo simulation by inaccurate numerical implementation of conceptual hydrologic models, Water Resour. Res., 46, W10530, https://doi.org/10.1029/2009wr008648, 2010. a, b, c, d
Seiller, G., Roy, R., and Anctil, F.: Influence of three common calibration metrics on the diagnosis of climate change impacts on water resources, J. Hydrol., 547, 280–295, https://doi.org/10.1016/j.jhydrol.2017.02.004, 2017. a
SzöllösiNagy, A.: The discretization of the continuous linear cascade by means of state space analysis, J. Hydrol., 58, 223–236, https://doi.org/10.1016/00221694(82)900361, 1982. a, b
van Esse, W. R., Perrin, C., Booij, M. J., Augustijn, D. C. M., Fenicia, F., Kavetski, D., and Lobligeois, F.: The influence of conceptual model structure on model performance: a comparative study for 237 French catchments, Hydrol. Earth Syst. Sci., 17, 4227–4239, https://doi.org/10.5194/hess1742272013, 2013. a
Vidal, J.P., Martin, E., Franchisteguy, L., Baillon, M., and Soubeyroux, J.M.: A 50year and highresolution atmospheric reanalysis over and France with the Safran system, Int. J. Climatol., 30, 1627–1644, https://doi.org/10.1002/joc.2003, 2010. a
Wood, E. F., Lettenmaier, D. P., and Zartarian, V. G.: A landsurface hydrology parameterization with subgrid variability for general circulation models, J. Geophys. Res., 97, 2717–2728, https://doi.org/10.1029/91JD01786, 1992. a, b
Young, P. and Garnier, H.: Identification and estimation of continuoustime, databased mechanistic (DBM) models for environmental systems, Environ. Modell. Softw., 21, 1055–1072, https://doi.org/10.1016/j.envsoft.2005.05.007, 2006. a