# Conceptual rainfall-runoff model with a two-parameter, infinite characteristic time transfer function

Vincent Guinot, Marie Savéan, Hervé Jourde, Luc Neppel

► To cite this version:

Vincent Guinot, Marie Savéan, Hervé Jourde, Luc Neppel. Conceptual rainfall-runoff model with a two-parameter, infinite characteristic time transfer function. Hydrological Processes, 2015, 29 (22), pp.4756-4778. 10.1002/hyp.10523. hal-01143608

HAL Id: hal-01143608

https://hal.science/hal-01143608

Submitted on 18 Apr 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Conceptual rainfall-runoff model with a two-parameter,
infinite characteristic time transfer function

Vincent Guinot* (1, 2), Marie Savéan (1), Hervé Jourde (1), Luc Neppel (1)

*: Corresponding author. Tel +33 (0)4 67 14 90 56, Fax +33 (0)4 67 14 47 74, e-mail: vincent.guinot@um2.fr

1. Université Montpellier 2, UMR HSM, Place Eugène Bataillon, 34095 Montpellier Cedex 5,
France

2. Inria, Team LEMON, Bât 5, CC05 017, 860 rue Saint-Priest, 34095 Montpellier Cedex 5,
France

Abstract

A two-parameter transfer function with an infinite characteristic time is proposed for con-
ceptual rainfall-runoff models. The large time behaviour of the unit response is an inverse
power function of time. The infinite characteristic time allows long term memory effects to be
accounted for. Such effects are observed in mountainous and karst catchments. The governing
equation of the model is a fractional differential equation in the limit of long times. Although
linear, the proposed transfer function yields discharge signals that can usually be obtained
only using non-linear models. The model is applied successfully to two catchments, the Dud
Koshi mountainous catchment in the Himalayas and the Durzon karst catchment in France. It
compares favourably to the linear, non-linear single reservoir models and to the GR4J model.
With a single reservoir and a single transfer function, the model is capable of reproducing
hysteretic behaviours identified as typical of long term memory effects. Computational effi-
ciency is enhanced by approximating the infinite characteristic time transfer function with a
sum of simpler, exponential transfer functions. This amounts to partitioning the reservoir into
several linear subreservoirs, the output discharges of which are easy to compute. An efficient
partitioning strategy is presented to facilitate the practical implementation of the model.

# 1 Introduction

The hydrological response of a number of natural systems, such as karst and mountainous catch-
ments, is well-known to involve multiple time scales. Such catchments typically respond to the
precipitation signal in the form of very fast and sharp discharge peaks, followed with long and
slowly decreasing base discharge signals. Besides, the higher discharge peaks are not necessarily
triggered by the higher precipitation intensities (see e.g. Latron et al., 2008), which is a clear
indication of long-term memory effects. As far as karst catchments are concerned, such effects
are attributed to the dual role (storage and propagation) of the epikarst and conduits (Jukic and
Denic-Jukic, 2009). Catchments with steep slopes may exhibit a bimodal runoff response to rainfall
events, a feature that is attributed to parallel responses of surface and subsurface compartments
(Graeff et al., 2009). The consequence is a variable response time to the recharge signal (Delbart
et al., 2014). Accounting for multiple transfer time scales simultaneously thus appears as a highly
desirable feature for rainfall-runoff models (Terzić et al., 2012).

Three classical modelling approaches to reconstruct such behaviours are (i) increasing the num-
ber of reservoirs, (ii) introducing a degree of non-linearity in the response of the reservoir(s) in the
form of threshold/power discharge laws, hysteretic behaviours, etc., (iii) increasing the complexity
of the transfer functions by introducing additional variables or time-dependent transfer functions,
in the form of unit hydrographs (Chow et al., 1988; Dooge, 1959; Jakeman et al., 1990; Nash, 1957;
Sherman, 1932) or model state-dependent functions. Most models available from the literature
combine several of these approaches (Birkel et al., 2010; Campbell and Sullivan, 2002; Chen and
Goldscheider, 2014; Edijatno et al., 1999; Fleury et al., 2007; Graeff et al., 2012; Hartmann et al,
2012, 2013a-b; Jukic and Denic-Jukic, 2009; Padilla and Pulido-Bosch, 2008; Perrin et al., 2003;
Tritz et al., 2011; Yue and Hashino, 2000). The price to pay is a substantial increase in the number
of model parameters, which is well-known to improve model calibration but does not necessarily help in increasing the predictive power of the model (Perrin et al., 2001).

All these approaches share the common property of yielding finite characteristic times, thus suggesting a limited potential for taking long-term memory effects into account.

This paper presents an alternative approach, whereby the unit response of the reservoir model is an inverse power function of time. Experimental evidence for such behaviours is suggested by Majone et al. (2004) for karst catchments. A salient feature of such a response is that it allows for infinite characteristic times. While inverse power responses may be obtained using non-linear reservoir equations, the proposed model is linear. Its flexibility makes it applicable to a wide range of situations.

The idea of modelling transfer processes with infinite characteristic times or lengths is not new. Infinite characteristic times are usually a signature for trapping phenomena in the transport/signal transfert process (Metzler and Klafter, 2000). They are obtained under the assumption of heavy-tailed Probability Density Functions (PDFs) for signal and/or matter transfer times and distances. An important class of heavy-tailed PDF is the family of PDFs that become equivalent to inverse power laws for large times. Random variables obeying such PDFs verify the generalized central limit theorem. Analyzing such PDFs within the Continuous Time Random Walk (CTRW) formalism (Klafter et al., 1987) leads to fractional differential governing equations. Transport processes obeying such governing equations are termed “anomalous” as opposed to “normal”, finite characteristic time-based processes. Anomalous transport processes obeying such laws have been identified in turbulent flows (Richardson, 1926), the migration of electric charges in amorphous solids (Scher and Montroll, 1975; Pfister and Sher, 1977), biological transfer processes (Thurner et al., 2003), epidemics dynamics (Da Silva et al., 2013), also see the review by Metzler and Klafter (2000). As far as environmental processes are concerned, anomalous transport models have been used to model the migration of contaminants in porous media (see e.g. Levy and Berkowitz, 2003; Gao et al., 2009) and complex geological settings (Meerschraert et al., 1999; Sun et al., 2014). To our best knowledge, the idea of incorporating heavy-tailed, inverse power law probability density functions to model the propagation of a hydrological signal in a lumped conceptual model has not been proposed before. The purpose of this paper is to examine how beneficial such an approach can be to catchment modelling.

In Section 2, the principle of a unit response that behaves as an inverse power of time for asymptotically long times is introduced. The governing equation of the reservoir model is shown to obey a fractional differential equation for asymptotically long time scales. Section 3 shows how such a transfer function can be easily implemented in the structure of a classical reservoir model. Section 4 illustrates the application of the proposed model to two different catchments: the Dud Koshi catchment (Himalayas mountain chain) and the Durzon karst catchment. In both applications, the model outperforms classical linear and non-linear approaches based on a finite transfer time. Section 5 is devoted to conclusions.

# 2 Single reservoir model

## 2.1 Structure and general governing equations

Consider the single linear reservoir model with the following structure (Figure 1):

- a single input in the form of a prescribed inflow rate $r(t)$ as a function of time, independent of the height (or specific volume) $h$ stored in the model. In the field of surface catchment modelling, the inflow is the net rainfall. In the field of subsurface hydrological modelling, the term recharge is more often used,

- a single output $q(t)$, that is not necessarily a function of $h$, but depends on the transfer time of the inflow signal through the reservoir. A Dirac-shaped inflow signal entering the reservoir leaves the reservoir after a time (called waiting time or transfer time hereafter) $t$. The Probability Density Function (PDF) of the transfer time, or unit response of the reservoir, is denoted by $w(t)$ hereafter, with the constraint

$$ \int_{0}^{+\infty} w(t) dt = 1 \qquad (1) $$

and the characteristic time is defined as

$$ T_c = \frac{\int_{0}^{+\infty} tw(t) dt}{\int_{0}^{+\infty} w(t) dt} = \int_{0}^{+\infty} tw(t) dt \qquad (2) $$
The outflowing discharge is related to the inflow by the convolution (hence the term “convolution kernel” for $w(t)$)

$$q(t) = r * w(t) = \int_{0}^{t} r(t) w(t - \tau) d\tau \quad (3)$$

and the specific volume $h$ obeys the balance equation

$$\frac{dh}{dt} = r - q \quad (4)$$

Applying the Laplace transform to the above two equations yields

$$\hat{q} = \hat{r}\hat{w} \quad (5a)$$

$$\hat{h} = \frac{1 - \hat{w}}{s} \hat{r} \quad (5b)$$

where the $\hat{\bullet}$ operator denotes the Laplace transform and $s$ is the Laplace variable. The consequences of a finite and infinite characteristic time are examined in the next sections.

## 2.2 Finite characteristic time transfer function

Characteristic time scales for the unit response $w(t)$ are typical of Markovian processes, whereby the behaviour of the system at a given time is independent of its behaviour at previous times. Markovian processes yield exponential-like behaviours for $w(t)$ (Klafter et al., 1987). The simplest possible example of a Markovian system with a linear response is the linear reservoir model (Maillet, 1906). The unit response $w(t)$ and characteristic time $T_c$ of the linear reservoir model are

$$w(t) = k \exp(-kt) \quad (6a)$$

$$T_c = \frac{1}{k} \quad (6b)$$

where $k$ is the discharge constant of the reservoir. Placing linear reservoirs in series or in parallel may modify the unit response of the system, but in any case the characteristic time is finite. Applying the Laplace transform to (6a) gives

$$\hat{w} = \frac{k}{k + s} \quad (7)$$

and (5b) becomes

$$(k + s) \hat{h} = \hat{r} \quad (8)$$

Reverting from the Laplace variable to the classical time variable yields the well-known linear reservoir equation:

$$\frac{dh}{dt} = r - kh \quad (9)$$

## 2.3 Infinite characteristic time transfer function

Assume now that the unit response $w$ obeys the following asymptotic power law:

$$w(t) \sim \frac{A}{t^{1+\alpha}}, \quad 0 < \alpha < 1 \quad (10)$$

where the exact value of the constant $A$ does not need to be known at this stage. The constraint (1) imposes that $A$ depends not only on $\alpha$ but also on the behaviour of $w$ at small times. The coefficient $\alpha$ cannot be negative, otherwise the integral in (1) does not converge. It cannot be larger than unity because $\alpha \ge 1$ yields a finite characteristic time scale. The assumption $0 < \alpha < 1$ is thus retained in what follows. Asymptotic behaviours in the form (10) are customary in the field of fractional dynamics and anomalous transport (Klafter et al., 1987; Metzler and klafter, 2000) and have been observed in a variety of natural processes. They usually express the influence of trapping phenomena inducing asymptotically long transfer times. They are characteristic of non-Markovian processes, with long-term memory effects that makes the behaviour of a system dependent of all its past states. The asymptotic behaviour (10) is supported by time series analyses on karst catchments (Majone et al., 2004).

The following formula is proposed for $w$ in the present application:
$$w(t) = \frac{\alpha \tau_0^\alpha}{(\tau_0 + t)^{1+\alpha}}, 0 < \alpha < 1 \qquad (11)$$

Although many other functions would be acceptable, equation (11) is interesting because it fulfils the two necessary conditions (1, 10), while remaining fairly simple. It involves only two parameters, the time scale $\tau_0$ and the exponent $\alpha$. Figure 2 shows the typical behaviour of the convolution kernel (11) compared to that of the linear reservoir (6a). The long term memory effect of the kernel (11) is clearly visible, with $w(t)$ decreasing much slower than the exponential response (6a).

**Asymptotic governing equation for large times.** The governing equations are written for a time scale such that the asymptotic behaviour (10) is valid. In this case, one has (see Appendix A for a proof)

$$\hat{w} \sim_{s \to 0} 1 - \beta s^{\alpha}, \quad \beta = \frac{A}{\alpha} \Gamma(1-\alpha) \qquad (12)$$

where the Gamma function is the extension of the factorial to real numbers. In the limit of small $s$, Equation (5b) becomes

$$\hat{h} = \beta s^{\alpha-1} \hat{r} \qquad (13)$$

Applying the inverse Laplace transform to the above equation gives the following fractional differential equation for large times

$$\frac{d^{1-\alpha}h}{dt^{1-\alpha}} = \beta r, \quad t \to +\infty \qquad (14)$$

where $\frac{d^{1-\alpha}}{dt^{1-\alpha}}$ is the non-integer fractional derivative (to be distinguished from the fractal derivative) defined in Appendix B, with lower integration bound $a = 0$. An alternative formulation is

$$\frac{dh}{dt} = r - q, \quad q = r - \beta \frac{d^{\alpha}r}{dt^{\alpha}}, \quad t \to +\infty \qquad (15)$$

It should be kept in mind that the formulations (14, 15) are only asymptotic versions of the balance equations for asymptotically large times. These equations are not to be used as governing equations for small time scales, where the approximation (10) does not hold.

# 3 Practical implementation

## 3.1 Implementation using local operators

### 3.1.1 General

The purpose is to find a computationally efficient way to implement the governing equation (4). Substituting equation (11) into (3) yields

$$q(t) = \alpha \tau_0^\alpha \int_0^t \frac{r(t-\tau)}{(\tau_0 + \tau)^{1+\alpha}} d\tau \qquad (16)$$

The computation of the integral (16) is CPU-intensive because the number of required arithmetic operations is proportional to the simulated time. The consequence is that simulating $N$ time steps requires a number of operations proportional to $N^2$. To give but an example, multiplying the simulated period by 10 yields a CPU time multiplied by 100. Using such a model over long time series would quickly become impracticable. This issue stems directly from the non-local character of the convolution (16), a property that shared with the fractional derivatives appearing the governing equation (14) in the limit of infinite times (see also Appendix B for the non-local character of the fractional derivative).

This computational burden can be eliminated by replacing the convolution (16) with a set of local operations, that is, operations involving only the current state of the system and not all its past states. Bearing this in mind, an alternative approach is proposed. It is based on the following remark: the convolution (16) does not need to be computed exactly for all times. In practice, the simulation time will never exceed an upper limit $T$. Therefore, the convolution (16) needs to be computed accurately only up to $t = T$. For this reason, it is proposed that the convolution kernel (11) be approximated with a set of functions that characterize a local behaviour. The simplest possible function is that of the linear reservoir (6a). The following approximation is proposed:
$$w(t) \approx \sum_{i=1}^{R} W_i(t), \quad W_i(t) = \frac{\theta_i}{k_i} \exp(-k_i t) \qquad (17)$$

where the $\theta_i$ are constant coefficients, the determination of which is explained hereafter. The outflow discharge is obtained as

$$q(t) \approx \sum_{i=1}^{R} \theta_i q_i(t), \quad q_i(t) = \int_{0}^{t} r(t-\tau) \frac{1}{k_i} \exp(-k_i \tau) d\tau \qquad (18)$$

This amounts to partitioning the reservoir into R linear sub-reservoirs placed in parallel (Figure 3). Note that this idea is not new. A similar approach is used in (Jakeman et al., 1990) for the representation of the unit hydrograph. The outflowing discharge q in (16) is the sum of the R discharges $q_i$ ($i = 1, \dots, R$) flowing out of the linear sub-reservoirs weighted by $\theta_i$. Each of these linear sub-reservoirs has a finite characteristic time $T_{c,i} = \frac{1}{k_i}$ (see equation (6b)). The behaviour of this system is local in that each discharge $q_i$ is a function of the current state of the sub-reservoir i alone and does not involve its past states. The governing equations for the model then become

$$h(t) = \sum_{i=1}^{R} \theta_i h_i(t), \quad q(t) = \sum_{i=1}^{R} \theta_i q_i(t) \qquad (19a)$$

$$q_i = k_i h_i \qquad (19b)$$

$$\frac{dh_i}{dt} = r(t) - k_i h_i \qquad (19c)$$

$$\sum_{i=1}^{R} \theta_i = 1 \qquad (19d)$$

where $k_i$ and $h_i$ are respectively the discharge constant and the specific volume for the ith sub-reservoir, while h is the specific volume for the entire catchment model.

Mass conservation stemming from equations (19c-19d), the only requirement on the model (19c) is that the discharge constants $k_i$ allow for a sufficiently accurate approximation of the convolution kernel $w(t)$ between $t = 0$ and $t = T$. How the partition should be carried out is explained in the next subsection.

### 3.1.2 Partitioning strategy

The accuracy of the approximation (17) is strongly conditioned by (i) the number R of linear sub-reservoirs, (ii) the choice of the discharge coefficients $k_i$, $i = 1, \dots, R$. The partitioning method is based on the following considerations.

1. Assuming that the discharge constants $k_i$, $i = 1, \dots, R$ are predefined, the only unknowns in the approximation (17) are the fractions $\theta_i$, $i = 1, \dots, R$. The $\theta_i$ can be determined uniquely by enforcing the approximation (17) for a set of R predefined times $t_j$, $j = 1, \dots, R$:

$$\sum_{i=1}^{R} \frac{\theta_i}{k_i} \exp(-k_i t_j) = \frac{\alpha \tau_0^\alpha}{(\tau_0 + t_j)^{1+\alpha}} \quad \forall t_j, j = 1, \dots, R \qquad (20)$$

and solving the system (20) for the $\theta_i$, $i = 1, \dots, R$.

2. A desirable behaviour is that the approximation (17) yield positive discharges $q_i(t)$ for all times. For this, the coefficients $\theta_i$ must all be positive, thus yielding a positive, monotonically decreasing approximation (17).

3. Condition 2 can be easily satisfied provided that the coefficients $k_i$ are sufficiently different from each other, so that the system (20) is correctly conditioned. An easy way of enforcing this is that the predefined times $t_j$ stemming from Consideration 1 are sufficiently different from each other and that the coefficients $k_i$ are related to these predefined times. The following option has proved successful for a wide range of $\alpha$:

$$w(t_j) = \left(\frac{R-j+1}{R}\right)^a w(0) \implies t_j = \left[\left(\frac{R}{R-j+1}\right)^{\frac{a}{1+\alpha}} - 1\right] \tau_0 \qquad (21a)$$

$$k_i = \frac{b}{\tau_0 + t_i} \qquad (21b)$$
Numerical experiments indicate that $R \in [10, 20]$, $a \in [5, 10]$ and $b \in [2, 3]$ yield accurate results. Three approximation examples of the law (11) with (17) are shown and commented in Appendix C.

## 3.2 Practical aspects

### 3.2.1 Need for an upper volume threshold

The partitioning method (17, 21a, 21b) implies that some of the $k_j$ will be extremely small. For instance, for $(R, a, b) = (15, 7.5, 2.5)$, one will have $t_{15} = 7.6 \times 10^5 \tau_0$ and $k_{15} = \frac{3.3 \times 10^{-6}}{\tau_0}$. Considering that $\tau_0$ is often of the order of a few days, this makes $t_{15}$ range from $10^3$ to $10^4$ years. This is far beyond the length of any realistic reservoir simulation, thus providing a very convincing justification for the approximation (17).

This, however, raises an issue. Considering an average inflow $\bar{r}$, the average height stored in the sub-reservoir $i$ is

$$ \bar{h}_i = \frac{\bar{r}}{k_i} \qquad (22) $$

For the above set of parameters, $\bar{r} = 1\,\text{m/yr}$ yields $\bar{h}_{15} = 830\,\text{m}$ for $\tau_0 = 1\,\text{day}$ and $\bar{h}_{15} = 2500\,\text{m}$ for $\tau_0 = 3\,\text{days}$. Clearly, such storage heights are unrealistic. It is thus necessary to bound the permissible specific volumes $h_i$ of the sub-reservoirs with an upper threshold $h_{\text{max}}$ that prevent an unrealistic storage volume in the slower linear sub-reservoirs. If the specific volume $h_i$ is larger than $h_{\text{max}}$, the inflow $r$ to the sub-reservoir $i$ is set to zero and the inflow is routed directly to the outlet, see Section 4.1 for the practical implementation of the inflow model.

This mechanism may bear several physical interpretations depending on the type of catchment under consideration. The instantaneous transfer of excess water to the outlet of the catchment may be seen as a model for Dunne-type runoff, a runoff generated after long and moderate rainfall events that lead to the progressive saturation of soils and the further impossibility for water to infiltrate. In the field of karst catchment modelling, it may be seen as the piston flow resulting from almost instantaneous pressure transfer when the conduit network comes to saturation. This mechanism is implemented in models such as the GR4J model (Perrin et al., 2003) and the Vensim model (Fleury et al., 2007) for karst catchment modelling.

### 3.2.2 Incorporating non-linearity: a straightforward approach

In classical reservoir models, two widespread approaches to non-linearity are threshold functions and power discharge laws. They may be viewed as modifications of the linear reservoir law. For instance, a power discharge law can be rewritten as

$$ q(h) = k \left( \frac{h}{h_{\text{max}}} \right)^{B-1} h \qquad (23) $$

where $B$ is a positive power (usually larger than unity) and $h_{\text{max}}$ is used for scaling purposes. This non-linearity may be viewed as the result of a volume-dependent characteristic time

$$ T_C(h) = \frac{1}{k} \left( \frac{h_{\text{max}}}{h} \right)^{B-1} \qquad (24) $$

$T_c$ being a decreasing function of the volume stored in the reservoir. This is in agreement with the assumption that the density of connected flow paths increases with the degree of saturation of the reservoir. This formula is easily generalized to the proposed transfer functions (19b) by setting

$$ q_i(h) = k_i \left( \frac{h}{h_{\text{max}}} \right)^{B-1} h_i \qquad (25) $$

For $B = 1$, the classical linear behaviour is recovered, while non-linearity is achieved for $B \neq 1$. Using the average specific volume $h$ as defined in (19a) applies the same time scaling to all the discharge constants of all the sub-reservoirs.

Note that the non-linear formula (23) yields output signals $q(t)$ in the form of inverse power functions of time, as does equation (11). However, there are two fundamental differences between the laws (11) and (23). Firstly, the non-linear nature of (23) makes the interpretation of the behaviour of the system in terms of unit response impossible. Secondly, (24) always yields finite characteristic times (with the only exception of an empty reservoir: $B > 1$, $h = 0$), while (11) always yields infinite characteristic times, even in the linear case.
### 3.2.3 Model initialization

A transfer function with an infinite characteristic time also raises the issue of model initialization. Formal and empirical sensitivity analyses of the response of reservoir models (Mazzilli et al., 2012) indicate that the initialization bias may exert a significant influence on the calibration and validation of conceptual models. The decay time of the initialization bias is proportional to the characteristic time $T_c$ of the transfer function. The warm-up period of the model (that is, the time needed to achieve a set of internal variables that are independent from the initial condition) is typically a few times $T_c$. In the case of an infinite characteristic time, it is theoretically impossible to eliminate completely the initialization bias because the length of the warm-up period is infinite.

A simple way to avoid using a long warm-up period consists in assuming that the long term effects of the inflow are strongly smoothed out over time and that the inflow may be replaced with its time average. A realistic approximation of the initial specific volume in the sub-reservoir $i$ is therefore given by the steady state solution

$$h_i(t=0) = \min \left( h_{\max}, \frac{\bar{r}}{k_i} \right) \qquad (26)$$

## 4 Application examples

### 4.1 Inflow model

In the application examples the model is also compared to the linear and non-linear reservoir with governing equation (23). The structure of the model is shown in Figure 4. Although Figure 4 is drawn for a subreservoir, the structure and functioning is exactly the same for the single linear and non-linear reservoir models.

The inflow model used in the application examples is very simple, especially in the way evaporation is considered. It is acknowledged that much more sophisticated models have been proposed in the literature (Fazal et al., 2005; Hartmann et al., 2013a; Jukic and Denic-Jukic, 2009). However, comparative studies have shown that very simple evapotranspiration models may give as satisfactory results as complex ones in the field of rainfall-runoff modelling (Oudin et al., 2005). It is stressed that the purpose of the present study is not to provide the best possible hydrological model but to illustrate the possibilities offered by the proposed convolution kernel (11).

The inflow model functions as follows.

- Each subreservoir receives the same net inflow rate $r(t)$. The net inflow $r(t)$ is computed as the difference between the precipitation rate $P$ multiplied by an infiltration coefficient $C_{inf}$ and the potential evapotranspiration rate PET.

- In the case of the Dud Koshi catchment presented in Section 4.2, snow and ice melt fluxes $M_s$ and $M_i$ are involved. They are added to the inflow without multiplication by the infiltration coefficient.

- The rate PET is subtracted from the reservoir only if the average depth $h = \sum_{i=1}^{R} h_i$ is above a minimum threshold $h_{min}$. For $h < h_{min}$, the water is assumed to be stored too deep in the system (lower part of the unsaturated zone and/or in the saturated zone) for plant uptake and/or evaporation.

- If the depth $h_i$ exceeds $h_{max}$ during the computational time step $\Delta t$, the subreservoir is not allowed to take up any more water. The net inflow is routed directly to the outlet of the subreservoir in the form of a spilling discharge $s_i$. As explained in Section 3.2, this direct routing may be seen as a model for Dunne-type runoff or pressure transfer-induced piston flow in karst aquifers.

$$r_i = \begin{cases} C_{inf}P + M_i + M_s & \text{if } h_i < h_{min} \\ C_{inf}P + M_i + M_s - \text{PET} & \text{if } h_{min} \le h_i < h_{max} + (C_{inf}P + M_i + M_s - \text{PET}) \Delta t \\ 0 & \text{if } h_i > h_{max} + (C_{inf}P + M_i + M_s - \text{PET}) \Delta t \end{cases} \qquad (27a)$$

$$s_i = \begin{cases} 0 & \text{if } h_i < h_{\max} + (C_{\inf}P + M_i + M_s - \text{PET}) \Delta t \\ C_{\inf}P + M_i + M_s - \text{PET} & \text{if } h_i = h_{\max} + (C_{\inf}P + M_i + M_s - \text{PET}) \Delta t \end{cases} \qquad (27b)$$
## 4.2 The Dud Koshi Himalayan catchment

### 4.2.1 Catchment description

The Dud Koshi catchment (3720 km²) is located in the central part of Himalayas in Nepal (Figure 5). It flows into the Koshi River, a tributary of the Ganges. The elevations rank from 700 m at Rabuwabazar in the south to the Mount Everest in the north. The relief is characterized by steep slopes and high orographic gradient. Three processes explain the Dud Koshi discharge: precipitation direct runoff, snow melt and ice melt, with glacier covering about 14% of the catchment area. The first one prevails during the monsoon season from June to September, the monsoon precipitations represent about 80% of the annual precipitation. The snow melt and the glacier melt lead to the inter season discharges from March to May and October to November. 20% of the annual precipitations come from the westerlies fluxes during the winter season. Atmospheric forcings are characterized by a strong spatial variability. The boundary of the snow cover area varies from 2000 m in winter to 4000 m during the monsoon. The maximum annual precipitation rate is observed at the elevation z = 2000 m with a maximum of about 2500 mm/year.

### 4.2.2 Data

A detailed description of the data can be found in (Savean et al., submitted), only an outline is given here. Forcings and discharges data at Rabuwabazar from 2001 to 2006 are used in this study, at a daily time step. The discharge time series are provided by the Department of Hydrology and Meteorology (DHM) of Nepal. The daily precipitation time series are computed from rain gauges measurement provided by DHM and EV-K2 CNR network, using a cokriging interpolation method. They are corrected by a piecewise constant multiplication factor that is adjusted for each hydrological year so as to preserve the annual hydrological balance. The underlying assumption of this correction is that the annual storage variations over the catchment are negligible compared to the integral of the fluxes.

The daily spatial air temperature, the net radiation, the relative humidity and the wind velocity time series are computed from NCEP/NCAR reanalysis data using a disaggregation procedure based on elevation. The Penman-Monteith method is used to estimate the potential evapotranspiration. The snow melt and ice melt fluxes have been estimated with the classical empirical degree-day model (Hock, 2003; Pokhrel et al., 2014).

### 4.2.3 Simulation results

The proposed model is compared to the linear and non-linear models, as well as to the GR4J model (Perrin, 2002; Perrin et al., 2003), a parsimonious rainfall-runoff model that functions at a daily time step. The GR4J model is considered standard in the field of catchment hydrology. It has two reservoirs and four parameters. The detailed description of the model structure and functioning, given in Perrin (2002) and Perrin et al. (2003), will not be recalled here for the sake of conciseness. It is important to stress, however, that the model has two reservoirs with different sizes $X_1$ and $X_3$, each of which obeys a non-linear discharge law in the form (23), and two unit hydrograph-based routing functions, with finite supports $X_4$ and $2X_4$. The parameter $X_2$ allows for mass balance corrections via exchanges through the catchment boundaries. It may thus be considered that the model incorporates four different time scales, two of which are fixed via $X_1$ and two of which are variable according to the considerations in Section 3.2.2. The GR4J model thus appears as a flexible tool to simulate the behaviour of catchments with complex structures.

The same calibration/validation methodology is used for the four models. Years 1 and 2 are used for model warm-up. Years 3 and 4 are used for calibration. Year 5 is used for validation. Calibration is carried out by a systematic exploration of the model parameter space and successive refinements around the global optimum. The objective function is the Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970). The authors are aware of the limitations associated with the NSE indicator (see e.g. Criss and Winston, 2008) and that many other objective functions could be used (Criss and Winston, 2008; Hogue et al., 2000; Krause et al., 2005; Legates and McCabe, 1999; Schaeffli and Gupta, 2007) for proposing themselves various generalized forms of objective functions (Guinot et al., 2011). The purpose here being to illustrate the flexibility of the model, the NSE is used only as a model benchmarking indicator. For each model, two different calibration options are used: calibration with an emphasis on peak flows, with the NSE computed using the discharges, and calibration with an emphasis on low flows, with the NSE computed using the transformed discharges $q^{1/4}$. There again, many other options could be used for low-flow model
fitting, but the NSE of the transformed flow variable has been used by a majority of authors (see the review by Pushpalatha et al. (2012)).

Table 1 shows the optimal calibration parameter sets obtained for the four models. The simulated discharges are shown in Figures 6 to 9. The origin of time $t = 0$ is taken at the beginning of the calibration period.

Figure 6 shows the simulation results for the linear model (6a). Calibrating for low flows gives more weight to the dry periods, with a large characteristic time ($T_C$ = 25 days) that allows the slowly rising base discharge to be better reconstructed in the beginning of the wet seasons (Figure 6a). This, however, is achieved at the expense of the reconstruction of the peak flows, that are considerably smoothed out by the model (Figure 6a). Calibrating for the peak flows yields a smaller value for the characteristic time ($T_C$ = 1/k = 4 days). The discharge peaks are much better reproduced (Figure 6b) than in Figure 6a, but reducing the characteristic time implies a fast decreasing discharge and an almost complete dryout of the model during the low flow periods. This is clearly the indication that a model with a single characteristic time cannot account satisfactorily for the behaviour of the catchment.

Figure 7 shows the simulation results for the non-linear model (23). Making the characteristic time volume-dependent allow the simulation results to be improved dramatically, both for low flow (Figure 7a) and peak flow calibration (Figure 7b). In particular, the low flow calibration yields a much better representation of the receding part of the hydrographs during the dry seasons.

Figure 8 shows the simulation results for the GR4J model. The model also succeeds in reproducing the succession of flow peaks and base flows. A better fit is obtained when the models calibrated against low flows. Compared to the peak flow calibration, the low flow calibration expectedly draws the simulated hydrograph closer to the observed one during the dry season. However, the model remains satisfactorily capable of simulating the fast peaks observed during the wet season. Strikingly enough, the optimal parameter set (see Table 1) is obtained for $X_1$ = 0 mm, while the optimal values for $X_3$ = 475 mm lies far outside the [200 mm – 300 mm] 80% confidence interval reported in Perrin et al. (2003). In other words, optimal model behaviour is obtained for an instantaneous percolation process. This may be considered as the sign that the vadoze zone plays a negligible role in catchment behaviour. The peak flow-calibrated $X_2$ = -10.4mm also lies beyond the range of usual values. Only this unusually high value allows the model to restore a correct water balance.

Figure 9 shows the simulation results obtained with the proposed model (11). The better NSE is obtained for the low flow calibration, a good compromise being achieved in terms of hydrograph reconstruction between the wet and dry periods. The peak flow calibration consistently underestimates the base discharge during the dry seasons. The NSE is systematically better than that of the non-linear reservoir models for both the calibration and validation period, regardless of the calibration option used. Its performance is very similar to that of the GR4J model, both in the calibration and validation phases.

Figure 10 shows scatter plots of the simulated specific discharges versus specific volumes in the three models. The linear and non-linear models expectedly exhibit one-to-one ($h, q$) relationships (Figures 10a, b). The GR4J model exhibits a bit more dispersion, especially when calibrated against peak flows (red dots in Figure 10c). Quite strikingly, the respective positions of the ($h, q$) sets are opposite to those of the single reservoir models. The low flow-calibrated GR4J model achieves smaller discharges for a given volume that the peak flow-calibrated model. This is related to the values for $X_2$, that is larger in the low flow-calibrated model than in the peak flow-calibrated one. The proposed model (Figure 10d) spans a much wider range of internal states than the other three. The resulting hysteretic behaviour directly stems from the long term memory effects embedded in equation (16).

## 4.3 The Durzon karst catchment

### 4.3.1 Catchment description

The Durzon karst catchment is located in the Southern Massif Central (France), as shown in Figure 11. The average ground elevation is 750 m ASL. The unsaturated and saturated zone are located in a 400 m thick aquifer made of Jurassic limestones and dolomites. The bottom, Southern and North-Eastern boundaries of the catchment are made of a 200 m thick marl layer (Bruxelles, 2001) that acts as a zero-flux boundary. The other parts of the catchment boundary are inferred from the topography. The average discharge at the main catchment outlet, the Durzon spring (533 m ASL) over the 2002-2007 period is $1.4 \text{ m}^3\text{s}^{-1}$, with peak discharge of about $18 \text{ m}^3\text{s}^{-1}$ (Bruxelles, 2001). The recharge area has been estimated between $100 \text{ km}^2$ and $120 \text{ km}^2$ (Bruxelles, 2001; Jacob et al.,
2008; Ricard and Bakalowicz, 1996). In the present study, a value of 116.8 km² is used as in (Tritz et al., 2011).

### 4.3.2 Data

The data used for the daily time step simulations is described in detail in (Tritz et al., 2011). Only its main features are recalled hereafter.

The precipitation history at “La Blaquèrerie” is not used for the present study because it is too short, the station being set up only recently. The precipitation record is taken from a meteorological station located at “Le Caylar” (SE corner of the map in Figure 11). Although outside the catchment, the daily time step precipitation records at Le Caylar were found to exhibit a satisfactory correlation with those measured at “La Blaquererie” and were used successfully in a previous study (Tritz et al., 2011).

The Potential Evapotranspiration is interpolated using a one year period harmonic function from monthly average values computed using Thornthwaite's model. The sinusoidal interpolation was shown to predict more realistic evapotranspiration values than the monthly averages (Tritz et al., 2011).

The daily outlet discharge is known indirectly from stage measurements at the Durzon spring. The stage-discharge relationship was shown in (Tritz et al., 2011) to exhibit an uncertainty of at least ±3% in the predicted discharge.

### 4.3.3 Simulation results

As in the previous application, Years 1-2 in the time series are used for model warm-up, Years 3-4 serve calibration purposes and Year 5 is used for model validation. Table 2 gives an overview of the calibration and validation results for the four models. Figures 12 to 15 show the simulated time series.

It should be stressed that calibrating and validating a model against the Durzon data set is an extremely challenging task. The calibration period ($0 \le t \le 730$) includes 4 moderately increasing discharge peaks triggered by moderate rainfall events, followed by a steadily decreasing discharge over a long period (almost 18 months), despite rather regularly spaced rainfall inputs. In contrast, the validation period ($t > 730$) includes an unusually high flow peak triggered by a very mild rainfall episode. This alone suggests a strong influence of past events and internal state in the response of the catchment, with storage presumably playing an essential role.

Figure 12 shows the results obtained with the linear reservoir model. Calibrating for low flows (Figure 12a) and peak flows (Figure 12b) yields very similar responses. Both calibration options yield dramatically underestimated peak discharges. This is not surprising in that the calibration period includes very few discharge peaks. The 18 months long dry spell contributes to bias the calibration process in favour of low flows and slow discharge release.

Figure 13 shows the simulation results for the non-linear model. The low flow calibration gives the linear model as the optimum ($B = 1$), see Figure 13a. The peak flow calibration yields a slightly improved representation of the receding part of the flow peaks compared to the linear model, but during the calibration period only.

Figure 14 shows the discharges simulated by the GR4J model. As for the Dud Koshi catchment, the model shows a good ability to reproduce the various time scales involved in catchment behaviour. It brings a tremendous improvement over the linear and non-linear models. The combination of the fast routing, unit hydrograph-based transfer functions (with $X_1=0.5$ day) and the non-linear behaviour of the lower reservoir is obviously responsible for this. Indeed, the non-linear transfer function alone does not allow the behaviour of the catchment to be reproduced (compare Figures 13 and 14). This may be seen as a confirmation that multiple time scales operating simultaneously lead to a more flexible model than a single time scale, even variable. It is visible in Table 2 that the optimal values of the parameters $X_1$ to $X_3$ are very different from the classically admitted orders of magnitude. The 80% confidence intervals are reported to be $X_1 \in [100\text{mm}, 1200\text{mm}]$, $X_2 \in [-5\text{mm}, +3\text{mm}]$ and $X_3 \in [20\text{mm}, 300\text{mm}]$ (Perrin et al., 2003).

Figure 15 shows the discharges simulated using the proposed long term memory model. The model gives very similar results when calibrated against low flows and peak flows. The consistency between the two parameter sets in Table 2 is remarkable. Among the four models, the proposed long term memory model is the one with the better NSE values during the calibration period. It does not perform as well as the GR4J model in the validation period. The main reason for this is that the GR4J model is more successful in reproducing the peak flows around $t = 770$ d and
$t \approx 1000$ d. In contrast, the proposed model is better in reproducing the base flow during the dry periods.

Figure 16 shows the scatter plots for the simulated $(h, q)$ variables in the four models. As in the Dud Koshi application, the first two models exhibit one-to-one $(h, q)$ relationships (Figure 14a, b). In contrast, both the GR4J model and the proposed model exhibit hysteretic behaviours (Figure 14c, d). The ability of the GR4J model to span a wider range of internal states than the linear and non-linear models is obvious. This also the case with the proposed model, that embeds many different time scales operating simultaneously. The satisfactory performance of the proposed model can again be explained by its ability to span a wider range of internal state configurations than the other two single reservoir models.

# 5 Discussion

The application examples in Section 4 tend to indicate that the proposed model can be helpful in the field of catchment modelling. However, one may be entitled to questioning the physical relevance of a hydrological model with an infinite characteristic time. Indeed, the concept of an infinite memory leads to wonder (i) whether an infinite characteristic time will not lead to a strong initialization bias, making the model worthless in practice and (ii) how realistic an infinite characteristic time transfer function is when all hydrosystems have been existing and will exist only for a finite time. Moreover, one may ask (iii) whether using the asymptotic version of the governing equations (14, 15) may be appropriate, as is classically done in fractional dynamics. These questions are addressed in separate subsections hereafter.

## 5.1 Sensitivity to initial conditions

For a linear model, the sensitivity of the outflowing discharge to the initial condition is by definition the PDF $w(t)$. For the linear model, the sensitivity decreases exponentially with time, following equation (6a). For the proposed model, the sensitivity to the initial condition obeys equation (11), with the asymptotic behaviour (10). For the non-linear model (23), the sensitivity of the output discharge to the initial conditions is examined as follows. The governing equation for this model is:

$$ \frac{dh}{dt} = -kh^B, \quad h(t=0) = h_0, \quad B > 1 \qquad (28) $$

with analytical solution

$$ h(t) = \frac{1}{[h_0^{1-B} + (B-1)kt]^{\frac{1}{B-1}}} \qquad (29a) $$

$$ q(t) = \frac{k}{[h_0^{1-B} + (B-1)kt]^{\frac{B}{B-1}}} \qquad (29b) $$

The sensitivity of $q$ to the initial condition $h_0$ is given by

$$ \frac{\partial q}{\partial h_0} = \frac{kB}{[h_0^{1-B} + (B-1)kt]^{\frac{2B-1}{B-1}}} h_0^{-B} \underset{t \to \infty}{\sim} \frac{kB}{[(B-1)kt]^{\frac{2B-1}{B-1}}} h_0^B \qquad (30) $$

Note that, this model being non-linear, the sensitivity of $q$ to the initial condition is dependent on the history of the inflows $r(\tau)$, $0 \le \tau \le t$. Comparing the sensitivity of a non-linear model to that of linear ones is thus a biased exercise.

For all models, the decay of the sensitivity of $q$ to $h_0$ with time is assessed using the dimensionless sensitivity factor:

$$ A_s(t) = \frac{\frac{\partial q}{\partial h_0}(t)}{\frac{\partial q}{\partial h_0}(0)} \qquad (31) $$

Table 3 gives the formulae for $A_s$ for the linear, non-linear and proposed models. The numerical values obtained at the end of the model warm-up period are also provided for the Dud Koshi and Durzon catchments, for both the low flow and peak flow calibration parameter sets. As shown by the table, the bias introduced by the initial condition for the proposed model decreases to less than 0.25% of its initial value after 2 years of warm-up.
## 5.2 Physical relevance of an infinite memory transfer function

Applying models with infinite memory to hydrosystems that have existed and are bound to exist for a finite time may also appear questionable, if not nonsense. Several answers may be given to such an objection.

Firstly, widely accepted non-linear models yield responses with infinite characteristic times under certain conditions, without triggering any special reluctance within the hydrologic community. Although the concept of a characteristic time for a non-linear process is questionable, it is stressed that the analytical solution (29b) for the non-linear model (28) also has an infinite characteristic time for $B > 1$.

Secondly, bearing in mind that the model is to be operated over a finite time interval, the unit response of the model does not strictly need to have an infinite characteristic time. The important point is that it behave over the simulation time interval as if its characteristic time were infinite. This remark is illustrated by the following PDF

$$W(t) = \begin{cases} \alpha\tau_0^\alpha (\tau_0 + t)^{-1-\alpha} & \text{if } t \le T \\ \alpha\tau_0^\alpha (\tau_0 + T)^{-1-\alpha} & \text{if } T < t \le T + \frac{\tau_0 + T}{\alpha}, \quad 0 < \alpha < 1 \\ 0 & \text{if } t > T + \frac{\tau_0 + T}{\alpha} \end{cases} \quad (32)$$

Over the time interval $[0, T]$, the PDF $W(t)$ is exactly equal to the transfer function $w(t)$ proposed in equation (11): It takes the value $w(T)$ over $[T, T + \frac{\tau_0 + T}{\alpha}]$ and is zero beyond this time so as to achieve a unit integral. A model using $W$ instead of $w$ will produce exactly the same results over the time frame $[0, T]$, although $W$ and $w$ have respectively finite and infinite characteristic times. In other words, only using a model over an infinite simulation period would allow the infinite memory concept to be validated or invalidated. The finite character of any simulation period makes this issue irrelevant.

Thirdly, the practical implementation of the proposed model is based on a weighted sum of exponentials, see Appendix C. It follows that the approximate PDF (17) has a finite characteristic time $T_c = \sum_i \theta_i/k_i$.

## 5.3 Relevance of a fractional differential governing equation

The third issue tackled in this discussion is whether using the asymptotic version (14) may turn out more profitable and convenient than computing the convolution (16), even approximately. Two arguments militate in favour of a negative answer.

Firstly, the fractional differential equation (14) and its equivalent (15) are asymptotic versions of the exact governing system (4, 16). They are valid only for very large times and should not be used at time scales for which the approximation (10) does not hold.

Secondly, the unit response of a system governed by equation (15) is physically unacceptable. Consider a unit, instantaneous recharge in the form of a Dirac pulse, $r(t) = \delta(t)$. The resulting unit response $q(t)$ is obtained from equation (15):

$$\begin{aligned} q(t) &= r - \beta \frac{d^\alpha r}{dt^\alpha} \\ &= \delta(t) - \beta \frac{d^\alpha \delta}{dt^\alpha} \\ &= \delta(t) - \frac{\beta}{\Gamma(-\alpha)} \int_0^t \frac{\delta(\tau)}{(t-\tau)^{1+\alpha}} d\tau \\ &= \delta(t) - \frac{\beta}{\Gamma(-\alpha)} \frac{1}{t^{1+\alpha}} \end{aligned} \quad (33)$$

where the extension (41) of Cauchy's formula to non-integers has been used. The unit response is made of two components: an instantaneous transfer of the input Dirac to the outlet of the catchment, followed by a negative, inverse power law function of time. In other words, after the input as been transferred instantaneously to the outlet, water is being taken back into the model through the outlet at a rate that is initially infinite. Besides, the integral of $q(t)$ does not converge at $0$ for $0 < \alpha < 1$, which means that mass conservation cannot be guaranteed.

# 6 Conclusions

In this paper, a transfer function with an infinite characteristic time has been proposed for conceptual rainfall-runoff modelling. With only two parameters, the function (11) proves very flexible. The asymptotic behaviour of the unit response function of the model is an inverse power of time. Such a response can usually be obtained only by making the discharge law of the model non-linear.
The present model, however, is linear. A model with a single reservoir is able to reproduce correctly the measured discharge signals of two very different catchments: the Dud Koshi catchment in Nepal and the Durzon karst catchment in France. In the long time limit, the model obeys a fractional differential equation.

The proposed model has five parameters: the infiltration coefficient $C$, the maximum reservoir size $h_{max}$, the evapotranspiration depth threshold $h_{min}$ and the two parameters $(\alpha, \tau_0)$ in the transfer function (11). Leaving the model structure unchanged, replacing the function (11) with more classical linear (6a) or non-linear (23) transfer functions is seen to yield significantly poorer model performance. This indicates that the proposed transfer function plays a key role in model performance and that the good quality of the results cannot be attributed to the inflow model alone. The proposed model exhibits similar performance to that of the four-parameter GR4J model, a classical lumped conceptual model with two reservoirs and two unit hydrograph-based routing functions. It is stressed that the GR4J model and the proposed one have very different structures. Moreover, the production functions of the two models are not comparable because evapotranspiration, that represents a substantial fraction of the mass balance, is not modelled in the same way.

The approximation of the convolution (16) with a set of local-in-time, first-order differential equations (19c) proves decisive in terms of computational efficiency. Had this approximation not been implemented, the CPU costs would have been proportional to the square of the simulation time, which would have made the model impractical from an operational point of view.

The long term memory induced by the convolution kernel (11) is best illustrated by Figures 10 and 16. The linear and non-linear reservoir models expectedly show a one-to-one correspondence between the specific volume and the outflowing discharge, while the proposed model exhibits a hysteretic behaviour. In other words, the output discharge for the proposed model does not obey an equation of state, while classical models do. This is a clear indication that the model incorporates long term memory effects. The shape of the hysteresis paths in the $(h,q)$ space are radically different in the Dud Koshi and the Durzon models, a further illustration of the flexibility of the proposed transfer function with only two parameters. This is to be compared with the four-parameter hysteretic transfer function proposed in (Tritz et al., 2011), where the shape of the hysteretic cycle is specified explicitly and does not incorporate as much flexibility as the present transfer function.

The proposed transfer function appears very promising in terms of response versatility. Many reasons can be given for this.

- The memory effects induced by the infinite characteristic time imply that the model response does not obey an equation of state (it is not only a function of the internal state variables), thus allowing complex processes such as flowpath connecticivity/recharge area variations to be better represented.

- The practical implementation in the form of parallel sub-reservoirs clearly illustrates the multiple time scale basis of this model, a feature that has been identified as highly desirable in both karst and mountainous catchment modelling (Graeff et al., 2009; Terzić et al., 2012).

- The subreservoir upper volume threshold should not be considered only as a numerical artefact and may bear a physical meaning. This threshold, combined with the small discharge constant of the slower sub-reservoirs, may also be useful in modelling the occasional activation of high-altitude springs in mountainous/karst catchments (Petrella et al., 2009).

Research perspectives include

- studying the effect of replacing hysteretic discharge laws such as that proposed in (Tritz et al., 2011) and unit hydrograph models such as that of GR4J with the proposed model, with a subsequent reduction in the number of parameters,

- exploring the possibility to propose other formulae than (11) with the same asymptotic long time behaviour (10), but with a different short-term behaviour,

- studying theoretically and experimentally the effect of placing the proposed transfer function in cascade with finite or other infinite characteristic time transfer functions,

- investigating whether incorporating a non-linear behaviour as suggested in Section 3.2.2 can lead to a further improvement in model performance,
- benchmarking systematically the proposed model against a variety of conceptual hydrological models over a variety of catchments in order to better assess the usefulness and limitations of the proposed approach.

## Acknowledgements

The financial support of the Université Montpellier 2 and the French Ministry of Research to the second author's work is gratefully acknowledged. The authors wish to express their gratitude to the Reviewer for many constructive comments, for suggesting the comparison with GR4J and triggering the discussion in Section 5.

## Appendix A. Proof for the fractional differential equations

Proving equations (14, 15) stems from the equivalence (12) under the assumption

$$ w(t) \underset{t \to +\infty}{\sim} \frac{A}{t^{1+\alpha}}, \quad 0 < \alpha < 1 \qquad (34) $$

Deriving equation (12) from the above equivalence is not straightforward because the Laplace transform of $t^{-1-\alpha}$ does not exist for $0 < \alpha < 1$. This is because the integral $\int_a^b t^{-1-\alpha}dt$ diverges as $a$ tends to zero. An approximate expression is thus sought in the form of a power series expansion around $s = 0$. The constraint (1) yields the following condition

$$ \hat{w}(s=0) = \int_{0}^{+\infty} w(t) e^{-0t} dt = \int_{0}^{+\infty} w(t) dt = 1 \qquad (35) $$

Moreover, the derivative with respect to the Laplace variable obeys the following rule:

$$ \frac{d\hat{w}}{ds} = -\hat{tw} \qquad (36) $$

The Laplace transform of the function $tw$ exists because $\alpha < 1$, thus yielding a definite integral as the lower bound tends to zero:

$$ \hat{tw} \underset{s \to 0}{\sim} A\Gamma(1-\alpha) s^{\alpha-1} \qquad (37) $$

therefore

$$ \frac{d\hat{w}}{ds} \underset{s \to 0}{\sim} -A\Gamma(1-\alpha) s^{\alpha-1} \qquad (38) $$

Integrating the above equation with respect to time with the initial condition (35) leads to equation (12). The alternative formulation (15) is obtained by substituting the asymptotic equivalence (12) into equation (5a) and returning from the Laplace to the time variable.

## Appendix B. Fractional derivatives

There are many different ways of defining fractional derivatives (Oldham and Spanier, 1974; Podlubny, 1999). Only a summary is presented hereafter.

A first approach consists in extending Cauchy's formula for the nth-order integral of a function $f$ cancelling at $a$

$$ a I_f^n(t) = \frac{1}{(n-1)!} \int_a^t (t-\tau)^{n-1} f(\tau) d\tau \qquad (39) $$

to non-integer integration orders as

$$ a I_f^\alpha(t) = \frac{1}{\Gamma(\alpha)} \int_a^t (t-\tau)^{\alpha-1} f(\tau) d\tau \qquad (40) $$

and considering the $\alpha$th-order derivative as the $(-\alpha)$th-order integral:

$$ a D_f^\alpha(t) = a I_f^{-\alpha}(t) = \frac{1}{\Gamma(-\alpha)} \int_a^t \frac{f(\tau)}{(t-\tau)^{1+\alpha}} d\tau \qquad (41) $$
The Riemann-Liouville approach (Liouville, 1832; Riemann, 1953) considers the αth-order derivative as the nth-order derivative of the (n - α)th-order integral:

$$a D_f^\alpha (t) = \frac{d^n}{dt^n} a I_f^{n-\alpha} (t) = \frac{1}{\Gamma(n-\alpha)} \frac{d^n}{dt^n} \int_a^t (t-\tau)^{n-\alpha} f(\tau) d\tau \quad (42)$$

while the Caputo derivative (Caputo, 1967, 1969) considers the αth-order derivative as the (n - α)th-order integral of the nth-order derivative:

$$a D_f^\alpha (t) = \frac{1}{\Gamma(n-\alpha)} \int_a^t \frac{1}{(t-\tau)^{\alpha+1-n}} \frac{d^n f}{d\tau^n} d\tau \quad (43)$$

(the lower integration bound a is set to zero in this paper). The particular case n = 1 yields

$$a D_f^\alpha (t) = \frac{1}{\Gamma(1-\alpha)} \int_a^t \frac{1}{(t-\tau)^\alpha} \frac{df}{d\tau} d\tau \quad (44)$$

with the advantage that this integral converges for α < 1 provided that f is differentiable over the integration interval.

The Grünwald-Letnikov definition is (Grünwald, 1867)

$$a D_f^\alpha (t) = \frac{d^\alpha f}{d(t-a)^\alpha} = \lim_{N \to 0} \left( \frac{N}{t-a} \right)^\alpha \frac{1}{\Gamma(-\alpha)} \sum_{i=0}^{N-1} \frac{\Gamma(i-\alpha)}{\Gamma(i+1)} f \left( t - \frac{i}{n}(t-a) \right) \quad (45)$$

A last possibility is to define the derivative from the Laplace transform:

$$D_f^\alpha = s^\alpha \hat{f} \quad (46)$$

where ⌜ and ⌜ denote respectively the Laplace and inverse Laplace transforms of a function.

Fractional derivatives differ from integer-order derivatives by a number of important features:

- the fractional derivative is a non-local operator. The derivative $a D_f^\alpha (t)$ depends on the behaviour of f over the entire interval $[a, t]$, in contrast with the classical derivative, that involves only the behaviour of f at t,

- different definitions may yield different results. As an example, the Grunwald-Letnikov fractional derivative of a constant function yields a power function (Oldham and Spanier, 1974; Podlubny, 1999), while the Caputo fractional derivative of a constant function is zero. Besides, the value of the derivative depends on the lower integration bound a,

- in contrast with integer-order derivatives, fractional-order derivatives commute under very specific conditions. The requirements for commutation are not the same for the Riemann-Liouville and the Caputo definition (see e.g. Podlubny (1999), Section 2.4). Besides, differentiating successively a function to orders α and γ does not yield the (α + γ)th-order derivative.

## Appendix C. Approximation of the convolution kernel (11) with exponentials

The approximation of the proposed convolution kernel (11) with a set of decreasing exponential functions (17) is illustrated with three different values of the exponent α shown in Table 4. The Table also gives the values of a, b and R used in equations (21a, 21b). The results are shown in Figure 17. The top graphs illustrate the behaviour of the interpolated functions for early times $0 \le t \le 10\tau_0$, while the bottom graphs illustrate the long-term behaviour over 6 orders of magnitude, for $\tau_0 \le t \le 10^6\tau_0$. The dots represent the exact kernel (11), the solid lines represent the approximation (17). As shown in the Figure, the approximation is extremely satisfactory for all three values of α up to $t = 10^5\tau_0$. Beyond this range, oscillations appear in the interpolated profiles (see bottom graphs with logarithmic coordinates). The amplitude of these oscillations, however, is extremely small (less than $10^{-7}w(0)$). Besides, for the two catchments presented in Section 4, the typical value found for $\tau_0$ after calibration is larger than 4 days, which gives $10^5\tau_0 > 10^3$ years. This is far longer than any realistic rainfall-runoff simulation. It can be concluded that the approximation (17) of the kernel (11) is extremely satisfactory for rainfall-runoff modelling purposes. The number of functions needed (hence the number of sub-reservoirs in the model), R = 15, makes the CPU effort moderate.
References

Birkel C, Tezlaff D, Dunn SM, Soulsby C. 2010. Towards a simple dynamic process conceptualization in rainfall-runoff models using multi-criteria calibration and tracers in temperate, upland catchments. Hydrological Processes 24: 260-275.

Bruxelles L. 2001. Modélisation hydrologique distribuée sur le Gardon d'Anduze ; étude comparative de différents modèles pluie-débit, extrapolation de la normale à l'extrême et tests d'hypothèses sur les processus hydrologiques. Evolution morphogénétique, conséquences géologiques et implications pour l'aménagement. PhD thesis, University Aix-Marseille I, 266 pp. (in French).

Campbell CW, Sullivan SM. 2002. Simulating time-varying cave flow and water levels using the Storm Water Management Model. Engineering Geology 65, 133-139. Caputo, M., 1967. Linear model of dissipation whose Q is almost frequency-independent - II. Geophys. J.R. Astr. Soc. 13: 529-539.

Caputo M. 1969. Elasticita e dissipazione. Zanichelli, Bologna.

Chen Z, Goldscheider N. 2014. Modeling spatially and temporally varied hydraulic behavior of a folded karst system with dominant conduit drainage at catchment scale, Hochifengottesacker, Alps. Journal of Hydrology 514: 41-52.

Chow VT, Maidment DR, Mays LW. 1988. Unit Hydrograph. Applied Hydrology, Mc Graw-Hill, NY.

Criss RE, Winston WE. 2008. Do Nash values have value? Discussion and alternate proposals. Hydrological Processes 22: 2723-2725.

Delbart C, Valdes D, Barbecot F, Tognelli A, Richon P, Couchoux L. 2014. Temporal variability of karst aquifer response time established by the sliding-windows cross-correlation method. Journal of Hydrology 511: 580-588.

Da Silva MB, Macedo-Filho A, Albuquerque EL, Serva M, Lyra ML, Fulco UL. Critical properties of a superdiffusive epidemic process. Physical Review E 87: 062108.

Dooge JCI. 1959. A general theory of the unit hydrograph. Journal of Geophysical Research 64: 241-256.

Edijatno NON, Yang X, Makhlouf Z, Michel C. 1999. GR3J: a daily watershed model with three free parameters. Hydrological Sciences Journal 44: 263-277.

Fazal MA, Imaizumi M, Ishida S, Kawachi T, Tsuchihara T. 2005. Estimating groundwater recharge using the SMAR conceptual model calibrated by genetic algorithm. Journal of Hydrology 303: 56-78.

Fleury P, Plagnes V, Bakalowicz M. 2007. Modelling of the functioning of karst aquifers with a reservoir model: Application to Fontaine de Vaucluse (South of France). Journal of Hydrology 345: 38-49.

Gao G, Zhan H, Feng S, Huang G, Mao X. 2009. Comparison of alternative models for simulating anomalous solute transport in a large heterogeneous soil column. Journal of Hydrology 377: 391-404.

Graeff T, Zehe E, Blume T, Francke T, Schröder B. 2012. Predicting event response in a nested catchment with generalized linear models and a distributed watershed model. Hydrological Processes 26: 3749-3769.

Graeff T, Zehe E, Reusser D, Lück E, Schröder B, Wenk G, John H, Bronstert A. 2009. Process identification through rejection of model structures in a mid-mountainous rural catchment: observations of rainfall-runoff response, geophysical conditions and model inter-comparison. Hydrological processes 23: 702-718.

Grünwald AK. 1867. Uber 'begrenzte' derivationene und deren anwendung. Z. Angew. Math. Phys. 12: 441.
Guinot V, Cappelaere B, Delenne C, Ruelland D. 2011. Towards improved criteria for hydrological model calibration: theoretical analysis of distance- and weak form-based functions. Journal of Hydrology 401: 1-13.

Hartmann A, Lange J, Vivo Aguado A, Mizyed N, Smiatek G, Kunstmann H. 2012. A multi-model approach for improved simulations of future water availability at a large Eastern Mediterranean karst spring. Journal of Hydrology, 468-469: 130-138.

Hartmann A, Barbera JA, Lange J, Andreo B, Weiler M. 2013a. Progress in the hydrologic simulation of time variant recharge areas of karst systems - Exemplified at a karst spring in Southern Spain. Advances in Water Resources 54: 149-160.

Hartmann A, Weiler M, Wagener T, Lange J, Kralik M, Humer F, Mizyed N, Rimmer A, Barbera JA, Andreo B, Butscher C, Huggenberger P. 2013b. Process-based karst modelling to relate hydrodynamic and hydrochemical characteristics to sytem properties. Hydrol. Earth System Sciences 17: 3305-3321.

Hock R. 2003. Temperature index melt modelling in mountain areas, Journal of Hydrology 282: 104-115.

Hogue TS, Sorooshian S, Gupta H, Holz A, Braatz D. 2000. A multistep automatic calibration scheme for river forecasting models. Journal of Hydrometeorology 1: 524-542.

Jacob T, Bayer R, Chery J, Jourde H, Le Moigne N, Boy JP, Hinderer J, Luck B, Brunet P. 2008. Absolute gravity monitoring of water storage variation in a karst aquifer on the Larzac plateau (Southern France). Journal of Hydrology 359: 105-117.

Jakeman AJ, Littlewood IG, Whitehead PG. 1990. Computation of the instantaneous unit hydrograph and identifiable component flows with application to small upland catchments. Journal of Hydology 117: 275-300.

Jukic D, Denic-Jukic V. 2009. Groundwater balance estimation in karst by using a conceptual rainfall-runoff model. Journal of Hydrology 373: 302-315.

Klafter J, Blumen A, Shlesinger MF. 1987. Stochastic pathway to anomalous diffusion. Physical Review A 35: 3081-3085.

Krause P, Boyle DP, Bäse F. 2005. Comparison of different efficiency criteria for hydrological model assessment. Advances in Geosciences 5: 89-97.

Latron J, Soler M, Llorens P, Gallar F. 2008. Spatial and temporal variability of the hydrological response in a small Mediterranean research catchment (Vallcebre, Eastern Pyrenees). Hydrological Processes 22: 775-787.

Legates DR, McCabe GJ. 1999. Evaluating the use of “goodness of fit” measures in hydrologic and hydroclimatic model validation. Water Resources Research 35: 233-241.

Levy M, Berkowitz B. 2003. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. Journal of Contaminant Hydrology 64: 203-226.

Liouville J. 1832. Mémoire : sur le calcul des différentielles à indices quelconques. Journal de l'Ecole Polytechnique 13: 71-162.

Maillet E. 1906. La vidange des systèmes de réservoirs. Annales Ponts et Chaussées, Mémoires et Documents, 21.

Majone B, Bellin A, Borsato A. 2004. Runoff generation in karst catchments: multifractal analysis. Journal of Hydrology 294: 176-195.

Mazzilli N, Guinot V, Jourde H. 2012. Sensitivity analysis of conceptual model calibration to initialization bias. Application to karst spring discharge models. Advances in Water Resources: 42, 1-16.

Meerschraert MM, Benson DA, Baümer B. 1999. Multidimensional advection and fractional dispersion. Physical Review E, 59: 5026-5028.

Metzler R, Klafter J. 2000. The random walk's guide to anomalous diffusion: a fractional dynamics approach. Physics Reports 339: 1-77.
Nash JE. 1957. The form of the instantaneous unit hydrograph. IAHS Publications 45: 112-121.

Nash JE, Sutcliffe JV. 1970. River flow forecasting through conceptual models. 1. A discussion of principles. Journal of Hydrology 10: 282-290.

Oldham KB, Spanier J. 1974. The fractional calculus. Academic Press, reprint Dover Publ. Inc. 2002.

Oudin L, Hervieu F, Michel C, Perrin C, Andréassian V, Anctil F, Loumagne C. 2005. Which potential evapotranspiration input for a lumped rainfall-runoff model? Part 2—Towards a simple and efficient potential evapotranspiration model for rainfall-runoff modelling. Journal of Hydrology 303: 290-306.

Padilla A, Pulido-Bosch A. 2008. Simple procedure to simulate karstic aquifers. Hydrological Processes 22: 1876-1884.

Perrin, C. 2002. Vers une amélioration d'un modèle global pluie-débit au travers d'une approche comparative. La Houille Blanche 6/7: 84-91.

Perrin C, Michel C, Andréassian V. 2001. Does a large number of parameters enhance model performance? Comparative assessment of common catchment model structures on 429 catchments. Journal of Hydrology 242: 275-301.

Perrin C, Michel C, Andréassian V. 2003. Improvement of a parsimonious model for streamflow simulation. Journal of Hydrology 279: 275-289.

Petrella A, Capuano P, Carcione M, Celico F. 2009. A high-altitude temporary spring in a compartmentalized carbonate aquifer: the role of low-permeability faults and karst conduits. Hydrological Processes 23: 3354-3364.

Pfister G, Sher H. 1977. Time-dependent electrical transport in amorphous solids: AS2Se3. Physical Review B 15: 2062-2083.

Podlubny I. 1999. Fractional differential equations. Academic Press.

Pokhrel BK, Chevallier P, Andréassian V, Tahir AA, Arnaud Y, Neppel L, Bajracharya OR, Budhathoki KP. 2014. Comparison of two snowmelt modelling approaches in the Dud Koshi basin (Eastern Himalayas, Nepal). Hydrological Sciences Journal 59: 1507-1518.

Pushpalatha R, Perrin C, Le Moine N, Andréassian V. 2012. A review of efficiency criteria suitable for evaluating low-flow simulations. Journal of Hydrology, 420-421: 171-182.

Ricard J, Bakalowicz M. 1996. Connaissance, aménagement et protection des ressources en eau du Larzac septentrional, Aveyron, France. BRGM Report No. R38953 (in French).

Richardson LF. 1926. Atmospheric diffusion shown on a distance neighbour graph. Proceedings of the Royal Society, London, 709-737.

Riemann B. 1953. Versuch Einer Allgemeinen Auffasung Der Integration Und Differentiation, the Collection Works of Bernhard Riemann, Dover, New York, NY, USA.

Savéan M, Delclaux F, Chevallier P, Gonga-Saholiariliva N, Sharma RR, Neppel L, Arnaud Y, submitted. Modeling the Dudh Koshi River basin (Nepal) and discussion on precipitation inputs. Journal of Hydrology.

Shchaefli B, Gupta HV. 2007. Do Nash values have value? Hydrological Processes 21: 1075-1080.

Sher H, Montroll HW. 1975. Anomalous transit-time dispersion in amorphous solids. Physical Review B, 12: 2455-2477.

Sherman LK, 1932. Stream flow from rainfall by the unit graph method. Engineering News-Record 108: 501-505.

Sun H, Zhang Y, Chen W, Reeves DM. 2014. Use of a variable-index fractional-derivative model to capture transient dispersion in heterogeneous media. Journal of Contaminant Hydrology 157: 47-58.
Terzić J, Stroj A, Frangen T. 2012. Hydrogeological investigation of karst system properties by common use of diverse methods: a case study of Lička Jesenica springs in Dinaric karst of Croatia. Hydrological Processes 26: 3302-3311.

Thurner S, Wick N, Hanel R, Sedivy R, Huber L. 2003. Anomalous diffusion on dynamical networks: a model for interacting epithelial cell migration. Physica A 320: 475-484.

Tritz S, Guinot V, Jourde H. 2011. Modelling the behaviour of a karst system catchment using non-linear hysteretic conceptual model. Journal of Hydrology, 397: 250-262.

Yue S, Hashino M. 2000. Unit hydrographs to model quick and slow runoff components of streamflow. Journal of Hydrology 227: 195-206.
Tables

<table><caption>Table 1: Dud Koshi catchment. Simulation results for the single reservoir model. N.A.: Not Applicable.</caption><thead><tr><th rowspan="2">Model<br>Calibration<br>criterion</th><th colspan="2">Linear</th><th colspan="2">Non-linear</th><th colspan="2">GR4J</th><th colspan="2">Proposed model</th></tr><tr><th>Low flows</th><th>Peak<br>flows</th><th>Low flows</th><th>Peak<br>flows</th><th>Low flows</th><th>Peak<br>flows</th><th>Low flows</th><th>Peak<br>flows</th></tr></thead><tbody><tr><td>Figure</td><td>6a, 9a</td><td>6b, 9a</td><td>7a, 9b</td><td>7b, 9b</td><td>8a, 9b</td><td>8b, 9c</td><td>9a, 10b</td><td>9b, 10c</td></tr><tr><td>B (-) eq. (23)</td><td>1</td><td>1</td><td>5.0</td><td>6.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td></tr><tr><td>C<sub>inf</sub> (-) eqs. (27a,<br>27b)</td><td>0.6</td><td>0.7</td><td>0.75</td><td>0.72</td><td>N.A.</td><td>N.A.</td><td>0.65</td><td>0.75</td></tr><tr><td>h<sub>max</sub> (m) eqs. (26,<br>24, 27a, 27b)</td><td>10.</td><td>0.12</td><td>1.0</td><td>0.50</td><td>N.A.</td><td>N.A.</td><td>0.65</td><td>0.30</td></tr><tr><td>h<sub>min</sub> (m) eqs. (27a,<br>27b)</td><td>1.0</td><td>1.9 ×<br>10<sup>-2</sup></td><td>0.50</td><td>0.30</td><td>N.A.</td><td>N.A.</td><td>0.39</td><td>9 ×<br>10<sup>-2</sup></td></tr><tr><td>1/k (d) eq. (6a)</td><td>25.0</td><td>4.0</td><td>6.0</td><td>3.5</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td></tr><tr><td>α (-) eq. (11)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>0.5</td><td>0.8</td></tr><tr><td>τ<sub>0</sub> (d) eq. (11)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>9.0</td><td>4.0</td></tr><tr><td>GR4J: X<sub>1</sub> (mm)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>0.00</td><td>0.00</td><td>N.A.</td><td>N.A.</td></tr><tr><td>GR4J: X<sub>2</sub> (mm)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>-2.91</td><td>-10.4</td><td>N.A.</td><td>N.A.</td></tr><tr><td>GR4J: X<sub>3</sub> (mm)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>1270</td><td>475</td><td>N.A.</td><td>N.A.</td></tr><tr><td>GR4J: X<sub>4</sub> (d)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>0.94</td><td>1.05</td><td>N.A.</td><td>N.A.</td></tr><tr><td>NSE calibration</td><td>0.81</td><td>0.83</td><td>0.91</td><td>0.84</td><td>0.88</td><td>0.87</td><td>0.92</td><td>0.86</td></tr><tr><td>NSE validation</td><td>0.58</td><td>0.76</td><td>0.83</td><td>0.78</td><td>0.81</td><td>0.86</td><td>0.86</td><td>0.81</td></tr></tbody></table>

<table><caption>Table 2: Durzon catchment. Simulation results. N.A.: Not Applicable.</caption><thead><tr><th rowspan="2">Model<br>Calibration<br>criterion</th><th colspan="2">Linear</th><th colspan="2">Non-linear</th><th colspan="2">GR4J</th><th colspan="2">Proposed model</th></tr><tr><th>Low flows</th><th>Peak<br>flows</th><th>Low flows</th><th>Peak<br>flows</th><th>Low flows</th><th>Peak<br>flows</th><th>Low flows</th><th>Peak<br>flows</th></tr></thead><tbody><tr><td>Figure</td><td>11a, 14a</td><td>11b,<br>14a</td><td>12a, 14b</td><td>12b,<br>14b</td><td></td><td></td><td>13a, 14c</td><td>13b,<br>14c</td></tr><tr><td>B (-) eq. (23)</td><td>1</td><td>1</td><td>1.0</td><td>5.0</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td></tr><tr><td>C<sub>inf</sub> (-) eqs. (27a,<br>27b)</td><td>0.3</td><td>0.8</td><td>0.3</td><td>0.8</td><td>N.A.</td><td>N.A.</td><td>0.65</td><td>0.67</td></tr><tr><td>h<sub>max</sub> (m) eqs. (26,<br>24, 27a, 27b)</td><td>70</td><td>1.0</td><td>50</td><td>2.0</td><td>N.A.</td><td>N.A.</td><td>15.</td><td>15.0</td></tr><tr><td>h<sub>min</sub> (m) eqs. (27a,<br>27b)</td><td>35</td><td>0.1</td><td>3.0</td><td>1.0</td><td>N.A.</td><td>N.A.</td><td>1.58</td><td>1.2</td></tr><tr><td>1/k (d) eq. (6a)</td><td>100</td><td>200</td><td>125</td><td>110</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td></tr><tr><td>α (-) eq. (11)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>0.45</td><td>0.55</td></tr><tr><td>τ<sub>0</sub> (d) eq. (11)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>9.5</td><td>14.5</td></tr><tr><td>GR4J: X<sub>1</sub> (mm)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>0.00</td><td>0.00</td><td>N.A.</td><td>N.A.</td></tr><tr><td>GR4J: X<sub>2</sub> (mm)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>-165</td><td>-124</td><td>N.A.</td><td>N.A.</td></tr><tr><td>GR4J: X<sub>3</sub> (mm)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>4500</td><td>3296</td><td>N.A.</td><td>N.A.</td></tr><tr><td>GR4J: X<sub>4</sub> (d)</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>N.A.</td><td>2.13</td><td>2.10</td><td>N.A.</td><td>N.A.</td></tr><tr><td>NSE calibration</td><td>0.76</td><td>0.73</td><td>0.79</td><td>0.78</td><td>0.87</td><td>0.82</td><td>0.92</td><td>0.91</td></tr><tr><td>NSE validation</td><td>0.43</td><td>0.21</td><td>0.42</td><td>0.29</td><td>0.68</td><td>0.70</td><td>0.54</td><td>0.65</td></tr></tbody></table>
<table><caption>Table 3: Dimensionless sensitivity factor at the end of the warm-up period for the linear, non-linear and proposed infinite characteristic time models.</caption><thead><tr><th>Model</th><th>Formula</th><th>Dud<br>Koshi low<br>flows</th><th>Dud<br>Koshi<br>Peak<br>flows</th><th>Durzon<br>low flows</th><th>Durzon<br>peak<br>flows</th></tr></thead><tbody><tr><td>Linear</td><td>A<sub>s</sub> = exp(-kt)</td><td>2 &times; 10<sup>-13</sup></td><td>5 &times; 10<sup>-80</sup></td><td>6.7 &times; 10<sup>-4</sup></td><td>2.6 &times; 10<sup>-2</sup></td></tr><tr><td>Non-linear</td><td>A<sub>s</sub> = (<sup>h<sub>0</sub><sup>-1-B</sup></sup>&frasl;<sub>h<sub>0</sub><sup>-B</sup>+(B-1)kt</sub>)<sup>2B-1</sup>&frasl;<sub>B-1</sub></td><td>4.3 &times; 10<sup>-4</sup></td><td>4 &times; 10<sup>-4</sup></td><td>6.7 &times; 10<sup>-4</sup></td><td>0.11</td></tr><tr><td>Proposed<br>model</td><td>A<sub>s</sub> = (<sup>&tau;<sub>0</sub></sup>&frasl;<sub>&tau;<sub>0</sub>+t</sub>)<sup>1+&alpha;</sup></td><td>1.3 &times; 10<sup>-3</sup></td><td>8.4 &times; 10<sup>-5</sup></td><td>1.8 &times; 10<sup>-3</sup></td><td>2.2 &times; 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>Table 4: Approximation of the convolution kernel with a set of exponential functions. Sample parameter sets.</caption><thead><tr><th>Set name</th><th>&alpha;</th><th>a</th><th>b</th><th>R</th></tr></thead><tbody><tr><td>Set 1</td><td>0.1</td><td>10</td><td>2.5</td><td>15</td></tr><tr><td>Set 2</td><td>0.5</td><td>10</td><td>2.5</td><td>15</td></tr><tr><td>Set 3</td><td>0.9</td><td>10</td><td>2.5</td><td>15</td></tr></tbody></table>
Figure captions

– Figure 1. Single reservoir, single output model.

- Figure 2. Proposed convolution kernel (11). Comparison with the exponential unit response (6a) used in the linear reservoir equation. The grey-shaded area indicates the range of possible values for α. For the sake of comparison between the various curves, the axes are dimensionless. t* is the dimensionless time, t* = kt for the linear model (6a), t* = t/τ₀ for the proposed convolution (11). w* = w(t)/w(0) is the convolution kernel normalized by its value at t = 0.

- Figure 3. Partitioning an infinite characteristic time reservoir into a set of linear sub-reservoirs with finite characteristic times.

- Figure 4. Model structure. Definition sketch for a subreservoir.

– Figure 5. Dud Koshi catchment. Location map.

– Figure 6. Dud Koshi catchment. Simulation results for the linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.

- Figure 7. Dud Koshi catchment. Simulation results for the non-linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.

- Figure 8. Dud Koshi catchment. Simulation results for the proposed model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.

- Figure 9. Dud Koshi catchment. (*h*, *q*) scatter plots for the simulation period. (a) linear model, (b) non-linear model, (c) proposed model.

– Figure 10. Durzon karst catchment. Hydrogeological map. From Tritz et al. (2011).

– Figure 11. Durzon catchment. Simulation results for the linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.

- Figure 12. Durzon catchment. Simulation results for the non-linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.

- Figure 13. Durzon catchment. Simulation results for the proposed model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.

- Figure 14. Durzon catchment. (*h*, *q*) scatter plots for the simulation period. (a) linear model, (b) non-linear model, (c) proposed model.

- Figure 15. Approximation of the convolution kernel with a set of exponential functions for the parameter sets in Table 3. The dimensionless plot coordinates are defined as t* = t/τ₀ and w* = w(t)/w(0).
Figures

<figure><img src="image_2.png" /><figcaption>Figure 1: Single reservoir, single output model.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Proposed convolution kernel (11). Comparison with the exponential unit response (6a) used in the linear reservoir equation. The grey-shaded area indicates the range of possible values for $\alpha$. For the sake of comparison between the various curves, the axes are dimensionless. $t^*$ is the dimensionless time, $t^* = kt$ for the linear model (6a), $t^* = t/\tau_0$ for the proposed convolution (11). $w^* = w(t)/w(0)$ is the convolution kernel normalized by its value at $t=0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Partitioning an infinite characteristic time reservoir into a set of linear sub-reservoirs with finite characteristic times.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Model structure. Definition sketch for a subreservoir.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Dud Koshi catchment. Location map.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Dud Koshi catchment. Simulation results for the linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Dud Koshi catchment. Simulation results for the non-linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Dud Koshi catchment. Simulation results for the GR4J model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Dud Koshi catchment. Simulation results for the proposed model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Dud Koshi catchment. ($h, q$) scatter plots for the simulation period. (a) linear model, (b) non-linear model, (c) GR4J model, (d) proposed model. $R$ and $S$ are respectively the specific volumes stored in the upper and lower compartments of the GR4J model.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Durzon karst catchment. Hydrogeological map. From Tritz et al. (2011).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Durzon catchment. Simulation results for the linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Durzon catchment. simulation results for the non-linear reservoir model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Durzon catchment. simulation results for the GR4J model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Durzon catchment. Simulation results for the proposed model. (a) Calibration against low flows, (b) calibration against peaks, (c) forcings.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: Durzon catchment. (*h*, *q*) scatter plots for the simulation period. (a) linear model, (b) non-linear model, (c) proposed model.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17: Approximation of the convolution kernel with a set of exponential functions for the parameter sets in Table 3. The dimensionless plot coordinates are defined as $t^* = t/\tau_0$ and $w^* = w(t)/w(0)$.</figcaption></figure>