2.1 Reference GR4 model

GR4 (Perrin et al., 2003) is a lumped bucket-type 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.

Table 1Details of the equations of the GR4 model and discrete and continuous formulations. The discrete formulations are the continuous equations integrated individually over the modelling time step using the operator splitting technique while continuous equations correspond to the terms of the water balance differential equation of each store. * The values of UH₂ are calculated using Eq. (17) in Perrin et al. (2003). Please note that the two discrete formulations use either the unit hydrograph equations or the Nash cascade formulation.

Download Print Version | Download XLSX

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₉ in Fig. 1a) is routed by a nonlinear routing store (R in Fig. 1a). The complementary component (10 % of P_r, Q₁ in Fig. 1a) directly reaches the outlet. The groundwater exchange term (F) is added or removed from the routing store and from the Q₁ 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₁, x₂, x₃ and x₄) 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.

Table 2Meaning of the free and fixed parameters (from Perrin et al., 2003, except for U_t and nres).

Download Print Version | Download XLSX

2.2 A state-space formulation for the GR4 model

To create this state-space 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 state-space 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₄ 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⁻¹) 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₄ is the time to peak of the hydrograph.

The Nash cascade parameters are calculated depending on x₄ 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ösi-Nagy (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₄. By fixing the parameters in this way, only the x₄ 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₄ 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).

Figure 2Impulse response with x₄=2 time steps for the unit hydrograph of GR4 (dotted line) and the Nash cascade with nres=11 stores and $k=\frac{\mathrm{11}-\mathrm{1}}{x_{\mathrm{4}}}$ (dashed line).

Download

Using this formula, the x₄ 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₄ parameter makes it possible to reduce the possibility of an infinite impulse response.

2.4 Continuous differential equations of the state-space 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 state-space representation of the Nash cascade is the same as the one proposed by Szöllösi-Nagy (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.

Table 3Temporal transformations of the GR4 parameters (Ficchí et al., 2016).

Download Print Version | Download XLSX

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⁻¹). We decided to keep the interception out of the state-space representation because it is not represented by a store in the reference GR4J and we wanted to avoid introducing an additional difference between the state-space 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₉ and Q_r are, respectively, the inflow and the outflow of the routing store and F is the inter-catchment 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 state-space 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 state-space 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 sub-step 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 sub-step rather than single-step 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 single-step integration to those obtained with the adaptive sub-step 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 single-step 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.

Figure 3Location of the 240 flow gauging stations used for the tests and their associated catchments. The River Azergues at Châtillon is used as an example for the results (Sect. 4.1).

Download

3.2 Catchment set and data

To compare the performance and behaviour of the reference and the discrete and continuous state-space 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 split-sample test (Klemeš, 1986). These three versions are the reference model, a discrete state-space model (with a Nash Cascade but solved using OS) and a continuous state-space model. Comparing the reference and discrete state-space models allows us to measure the impact of replacing the unit hydrograph with a Nash cascade. Comparing the discrete and continuous state-space 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 2-year warm-up 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 low-flow 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 one-hundredth 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).

Figure 4Performance comparisons obtained in validation among the reference (with unit hydrograph), the discrete state-space (with Nash cascade) and the continuous state-space daily GR4 on 240 catchments, focusing on high (a), intermediate (b) and low (c) flows after calibration with the $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ (i.e. focusing on intermediate flow). The large points represent the mean performance and the smaller ones represent the outliers. The 5, 25, 50, 75 and 95 percentiles are represented by the box plots.

Download

Figure 5Simulated hydrograph of the River Azergues in the first half of 2012 during the validation period. The reference GR4 model (output in blue), the GR4 discrete state-space solution (output in green) and the continuous state-space solution (output in red) were calibrated with $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ as the objective function.

Download

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 split-sample 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 state-space 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 square-rooted flows. The performances of the reference model and the continuous state-space 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 state-space 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.

Figure 6Scatter plots of the four free parameters of the different versions of the models obtained by calibration with $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ as an objective function on the basins of the data set. Parameter comparison between unit hydrograph and Nash cascade is in black and parameter comparison between discrete and continuous state-space parameters is in red. The values of x₁, x₂ and x₃ are similar for the models (the line represents the y=x line). The x₄ values are higher in the discrete state-space model than for the other model versions.

Download

Figure 7Daily inputs in the routing store of the River Azergues in the first half of 2012. The models are calibrated with the $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ as the objective function. The peaks are lower with the discrete state-space GR4 (green lines) and occur sooner with the continuous state-space GR4 (red lines).

Download

Figure 8Daily routing store filling of the River Azergues in the first half of 2012. The reference GR4 (blue line) and the continuous state-space representation (red line) are calibrated with the $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ as the objective function.

Download

The study of the hydrographs provides complementary information. The reference GR4 model and the continuous state-space solution are very similar while the discrete state-space 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₄ parameter, which are systematically higher for the discrete state-space 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₄ values are due to errors caused by unsuitable solving is confirmed by the fact that the x₄ parameter values are similar for the three models at an hourly time step (not shown here).

Last, to understand the internal impact of the state-space formulation on the model, we analysed state variables and internal fluxes. Two differences are induced by the model's state-space formulation. First, the discrete Nash cascade output peaks are lower than the peaks of the unit hydrograph (Fig. 7). The peaks of the continuous state-space 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 state-space solution because the inputs in the routing store are too different for the discrete state-space solution. The peak levels are higher in the continuous state-space 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 high-flow conditions, the quantity of exchanged water and the outflow of the routing store in the discrete model are lower than those of the continuous state-space 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 state-space version due to the adaptive sub-step method. This is important to consider for some applications.

This computational time rise is essentially due to the adaptive sub-step algorithm. For example, in the River Azergues at Châtillon catchment, the mean number of sub-steps is 22 and it can reach 100 during some days. However, in Sect. 3.1 we argue that the adaptive sub-step method seems necessary to avoid numerical errors.

To conclude with these results, we can argue that the modifications brought by the continuous state-space 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 state-space representation of GR4 needs to be solved with a robust numerical technique.

4.2 Consistency of the state-space representation through time steps

The analysis of temporal consistency provides the most valuable result produced by the continuous state-space 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 state-space representation. The parameters of the state-space 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 state-space model.

Figure 9Scatter plots representing the four parameters of the reference (daily and hourly) GR4 models obtained by calibration with $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ as the objective function. The solid line represents the y=x regression and the dashed lines the transformation relations of Table 3.

Download

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₃ parameter follow the relations proposed by Ficchí et al. (2016) (the dashed lines). The high values of x₁ are underestimated compared to the theoretical relation as are the low values of the x₂ parameter. There is also an issue with the unit hydrograph parameter (x₄ 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₁, x₂ and x₄ 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₁ and the low values of x₂ are due to temporal inconsistencies in the interception calculation. The case of the x₄ parameter is more problematic. The differences in the x₄ values probably stem from the discretization of the unit hydrograph at different time steps.

In the continuous state-space 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₁ and the values of x₂ slightly diverge from the x=y line.

Figure 10Scatter plots representing the four parameters of the continuous state-space (daily and hourly) GR4 models obtained by calibration with $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ as the objective function. The solid line represents the y=x line.

Download

Figure 11Performance comparisons obtained in validation among the reference (with unit hydrograph), the discrete state-space (with Nash cascade) and the continuous state-space hourly GR4 on 240 catchments, focusing on high (a), intermediate (b) and low (c) flows after calibration with the $C_{\mathrm{2}M}\left(\sqrt{Q}\right)$ (i.e. focusing on intermediate flow). The points represent the mean performance.

Download

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₄ values because the x₄ 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 sub-step integration approximates a continuous time input in the Nash cascade. The results obtained with the x₄ 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₄ time consistency can also be explained by the integration errors produced by the OS at a daily time step.

The outliers in x₃ 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 non-sensitivity of the x₃ 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 state-space representation shows better temporal stability in the x₄ parameter values with similar performance.