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

A two-parameter transfer function with an in(cid:28)nite characteristic time is proposed for conceptual rainfall-runo(cid:27) models. The large time behaviour of the unit response is an inverse power function of time. The in(cid:28)nite characteristic time allows long term memory e(cid:27)ects to be accounted for. Such e(cid:27)ects are observed in mountainous and karst catchments. The governing equation of the model is a fractional di(cid:27)erential 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 identi(cid:28)ed as typical of long term memory e(cid:27)ects. Computational e(cid:30)-ciency is enhanced by approximating the in(cid:28)nite 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 e(cid:30)cient partitioning strategy is presented to facilitate the practical implementation of the model.


Introduction
The hydrological response of a number of natural systems, such as karst and mountainous catchments, 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 eects.As far as karst catchments are concerned, such eects are attributed to the dual role (storage and propagation) of the epikarst and conduits (Juki c and Deni c-Juki c , 2009).Catchments with steep slopes may exhibit a bimodal runo response to rainfall events, a feature that is attributed to parallel responses of surface and subsurface compartments (Grae 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-runo models (Terzić et al., 2012).
Three classical modelling approaches to reconstruct such behaviours are (i) increasing the number 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 complexiity 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;Grae et al., 2012;Hartmann et al, 2012Hartmann et al, , 2013a-b-b; Juki c and Deni c-Juki c, 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 nite characteristic times, thus suggesting a limited potential for taking long-term memory eects 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 innite characteristic times.While inverse power responses may be obtained using non-linear reservoir equations, the proposed model is linear.Its exibility makes it applicable to a wide range of situations.
The idea of modelling transfer processes with innite characteristic times or lengths is not new.
Innite 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 heavytailed 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 dierential governing equations.Transport processes obeying such governing equations are termed anomalous as opposed to normal, nite characteristic time-based processes.Anomalous transport processes obeying such laws have been identied in turbulent ows (Richardson, 1926), the migration of electric charges in amorphous solids (Scher and Montroll, 1975;Pster 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 benecial 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 dierential 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 dierent 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 nite transfer time.Section 5 is devoted to conclusions.
2 Single reservoir model

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 inow rate r (t) as a function of time, independent of the height (or specic volume) h stored in the model.In the eld of surface catchment modelling, the inow is the net rainfall.In the eld 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 inow signal through the reservoir.A Dirac-shaped inow 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 ˆ+∞ 0 w (t) dt = 1 (1) and the characteristic time is dened as The outowing discharge is related to the inow by the convolution (hence the term convolution kernel for w (t)) and the specic volume h obeys the balance equation Applying the Laplace transform to the above two equations yields q = r w where the • operator denotes the Laplace transform and s is the Laplace variable.The consequences of a nite and innite characteristic time are examined in the next sections.

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 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 nite.
Applying the Laplace transform to (6a) gives and (5b) becomes Reverting from the Laplace variable to the classical time variable yields the well-known linear reservoir equation:

Innite characteristic time transfer function
Assume now that the unit response w obeys the following asymptotic power law: 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 α but also on the behaviour of w at small times.The coecient α cannot be negative, otherwise the integral in (1) does not converge.It cannot be larger than unity because α ≥ 1 yields a nite characteristic time scale.The assumption 0 < α < 1 is thus retained in what follows.Asymptotic behaviours in the form (10) are customary in the eld 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 inuence of trapping phenomena inducing asymptotically long transfer times.They are characteristic of non-Markovian processes, with long-term memory eects 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: Although many other functions would be acceptable, equation ( 11) is interesting because it fulls the two necessary conditions (1, 10), while remaining fairly simple.It involves only two parameters, the time scale τ 0 and the exponent α. Figure 2 shows the typical behaviour of the convolution kernel (11) compared to that of the linear reservoir (6a).The long term memory eect 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 where the Gamma function is the extension of the factorial to real numbers.In the limit of small s, Equation (5b) becomes h = βs α−1 r (13) Applying the inverse Laplace transform to the above equation gives the following fractional dierential equation for large times where d 1−α dt 1−α is the non-integer fractional derivative (to be distinguished from the fractal derivative) dened in Appendix B, with lower integration bound a = 0.An alternative formulation is 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

General
The purpose is to nd a computationally ecient way to implement the governing equation (4).
Substituting equation ( 11) into (3) yields 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 innite 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: where the θ i are constant coecients, the determination of which is explained hereafter.The outow discharge is obtained as This amounts to partitioning the reservoir into R linear sub-reservoirs placed in parallel (Fig- ure 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 outowing discharge q in (16) is the sum of the R discharges q i (i = 1, . . ., R) owing out of the linear sub-reservoirs weighted by θ i .Each of these linear sub-reservoirs has a nite characteristic time T c,i = 1 ki (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 where k i and h i are respectively the discharge constant and the specic volume for the ith subreservoir, while h is the specic 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 suciently 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.

Partitioning strategy
The accuracy of the approximation ( 17) is strongly conditioned by (i) the number R of linear subreservoirs, (ii) the choice of the discharge coecients k i , i = 1, . . ., R. The partitioning method is based on the following considerations.
1. Assuming that the discharge constants k i , i = 1, . . ., R are predened, the only unknowns in the approximation (17) are the fractions θ i , i = 1, . . ., R. The θ i can be determined uniquely by enforcing the approximation (17) for a set of R predened times t j , j = 1, . . ., R: and solving the system (20) for the θ i , i = 1, . . ., R.
2. A desirable behaviour is that the approximation (17) yield positive discharges q i (t) for all times.For this, the coecients θ i must all be positive, thus yielding a positive, monotonically decreasing approximation (17).
3. Condition 2 can be easily satised provided that the coecients k i are suciently dierent from each other, so that the system (20) is correctly conditioned.An easy way of enforcing this is that the predened times t j stemming from Consideration 1 are suciently dierent from each other and that the coecients k i are related to these predened times.The following option has proved successful for a wide range of α: Three approximation examples of the law ( 11) with ( 17) are shown and commented in Appendix C.
Considering that τ 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 justication for the approximation (17).This, however, raises an issue.Considering an average inow r, the average height stored in the sub-reservoir i is For the above set of parameters, r = 1m/yr yields h 15 = 830 m for τ 0 = 1 day and h 15 = 2500 m for τ 0 = 3 days.Clearly, such storage heights are unrealistic.It is thus necessary to bound the permissible specic volumes h i of the sub-reservoirs with an upper threshold h max that prevent an unrealistic storage volume in the slower linear sub-reservoirs.If the specic volume h i is larger than h max , the inow r to the sub-reservoir i is set to zero and the inow is routed directly to the outlet, see Section 4.1 for the practical implementation of the inow 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 runo, a runo generated after long and moderate rainfall events that lead to the progressive saturation of soils and the further impossibility for water to inltrate.In the eld of karst catchment modelling, it may be seen as the piston ow 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.

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 modications of the linear reservoir law.For instance, a power discharge law can be rewritten as where B is a positive power (usually larger than unity) and h max is used for scaling purposes.This non-linearity may be viewed as the result of a volume-dependent characteristic time 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 ow paths increases with the degree of saturation of the reservoir.This formula is easily generalized to the proposed transfer functions (19b) by setting For B = 1, the classical linear behaviour is recovered, while non-linearity is achieved for B = 1.Using the average specic volume h as dened 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 dierences 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 nite characteristic times (with the only exception of an empty reservoir: B > 1, h = 0), while (11) always yields innite characteristic times, even in the linear case.

Model initialization
A transfer function with an innite 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 signicant inuence 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 innite characteristic time, it is theoretically impossible to eliminate completely the initialization bias because the length of the warm-up period is innite.
A simple way to avoid using a long warm-up period consists in assuming that the long term eects of the inow are strongly smoothed out over time and that the inow may be replaced with its time average.A realistic approximation of the initial specic volume in the sub-reservoir i is therefore given by the steady state solution However, comparative studies have shown that very simple evapotranspiration models may give as satisfactory results as complex ones in the eld of rainfall-runo 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 oered by the proposed convolution kernel (11).
The inow model functions as follows.
Each subreservoir receives the same net inow rate r (t).The net inow r (t) is computed as the dierence between the precipitation rate P multiplied by an inltration coecient 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 uxes M s and M i are involved.They are added to the inow without multiplication by the inltration coecient.
The rate PET is subtracted from the reservoir only if the average depth h = R i=1 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 ∆t, the subreservoir is not allowed to take up any more water.The net inow 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 runo or pressure transfer-induced piston ow in karst aquifers.

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 uxes.
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 uxes have been estimated with the classical empirical degree-day model (Hock, 2003;Pokhrel et al., 2014).

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-runo model that functions at a daily time step.The GR4J model is considered standard in the eld 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 dierent 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 nite 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 dierent time scales, two of which are xed 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 exible 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 renements around the global optimum.The objective function is the Nash-Sutclie Eciency (NSE) (Nash and Sutclie, 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;Schaei and Gupta, 2007) for proposing themselves various generalized forms of objective functions (Guinot et al., 2011).The purpose here being to illustrate the exibility of the model, the NSE is used only as a model benchmarking indicator.For each model, two dierent calibration options are used: calibration with an emphasis on peak ows, with the NSE computed using the discharges, and calibration with an emphasis on low ows, with the NSE computed using the transformed discharges q 1/4 .There again, many other options could be used for low-ow model tting, but the NSE of the transformed ow variable has bee 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 ows 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 ows, that are considerably smoothed out by the model (Figure 6a).Calibrating for the peak ows 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 ow 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 ow (Figure 7a) and peak ow calibration (Figure 7b).In particular, the low ow 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 ow peaks and base ows.A better t is obtained when the model s calibrated against low ows.Compared to the peak ow calibration, the low ow 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% condence 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 ow-calibrated X 2 = −10.4mmalso 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 ow calibration, a good compromise being achieved in terms of hydrograph reconstruction between the wet and dry periods.The peak ow 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 specic discharges versus specic 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 ows (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 ow-calibrated GR4J model achieves smaller discharges for a given volume that the peak ow-calibrated model.This is related to the values for X 2 , that is larger in the low ow-calibrated model than in the peak ow-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 eects embedded in equation ( 16).(Bruxelles, 2001) that acts as a zero-ux 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 m 3 s −1 , with peak discharge of about 18 m 3 s −1 (Bruxelles, 2001).The recharge area has been estimated between 100 km 2 and 120 km 2 (Bruxelles, 2001;Jacob et al., 2008;Ricard and Bakalowicz, 1996).In the present study, a value of 116.8 km 2 is used as in (Tritz et al., 2011).

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.

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 ≤ t ≤ 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 ow peak triggered by a very mild rainfall episode.This alone suggests a strong inuence 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 ows (Figure 12a) and peak ows (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 ows and slow discharge release.
Figure 13 shows the simulation results for the non-linear model.The low ow calibration gives the linear model as the optimum (B = 1), see Figure 13a.The peak ow calibration yields a slightly improved representation of the receding part of the ow 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 tranfer 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 conrmation that multiple time scales operating simultaneously lead to a more exible 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 dierent from the classically admitted orders of magnitude.The 80% condence intervals are reported to be X 1 ∈ [100mm, 1200mm], X 2 ∈ [−5mm, +3mm] and X 3 ∈ [20mm, 300mm] (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 ows and peak ows.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 ows around t = 770 d and t ≈ 1000 d.In contrast, the proposed model is better in reproducing the base ow 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 rst 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 dierent 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 congurations than the other two single reservoir models.

Discussion
The application examples in Section 4 tend to indicate that the proposed model can be helpful in the eld of catchment modelling.However, one may be entitled to questioning the physical relevance of a hydrological model with an innite characteristic time.Indeed, the concept of an innite memory leads to wonder (i) whether an innite characteristic time will not lead to a strong initialization bias, making the model worthless in practice and (ii) how realistic an innite characteristic time transfer function is when all hydrosystems have been existing and will exist only for a nite 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.

Sensitivity to initial conditions
For a linear model, the sensitivity of the outowing discharge to the initial condition is by denition 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: with analytical solution The sensitivity of q to the initial condition h 0 is given by Note that, this model being non-linear, the sensitivity of q to the initial condition is dependent on the history of the inows r (τ ) , 0 ≤ τ ≤ 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: 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 ow and peak ow 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.

Physical relevance of an innite memory transfer function
Applying models with innite memory to hydrosystems that have existed and are bound to exist for a nite time may also appear questionable, if not nonsense.Several answers may be given to such an objection.
Firstly, widely accepted non-linear models yield reponses with innite 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 innite characteristic time for B > 1.
Secondly, bearing in mind that the model is to be operated over a nite time interval, the unit response of the model does not strictly need to have an innite characteristic time.The important point is that it behave over the simulation time interval as if its characteristic time were innite.
This remark is illustrated by the following PDF 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 + τ0+T α 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 nite and innite characteristic times.In other words, only using a model over an innite simulation period would allow the innite memory concept to be validated or invalidated.The nite 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 nite characteristic time T c = i θ i /k i .

Relevance of a fractional dierential governing equation
The third issue tackled in this discussion is whether using the asymptotic version ( 14) may turn out more protable and convenient than computing the convolution ( 16), even approximately.Two arguments militate in favour of a negative answer.
Firstly, the fractional dierential equation ( 14) and its equivalent (15) are asymptotic versions of the exact geoverning 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 bu equation ( 15) is physically inacceptable.
Consider a unit, instantaneous recharge in the form of a Dirac pulse, r (t) = δ (t).The resulting unit reponse q (t) is obtained from equation (15): where the extension (41) of Cauchy's formula to non-integers has be 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 throught the outlet at a rate that is initially innite.Besides, the integral of q (t) does not converge at 0 for 0 < α < 1, which means that mass conservation cannot be guaranteed.

Conclusions
In this paper, a transfer function with an innite characteristic time has been proposed for conceptual rainfall-runo modelling.With only two parameters, the function ( 11) proves very exible.
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 dierent catchments: the Dud Koshi catchment in Nepal and the Durzon karst catchment in France.In the long time limit, the model obeys a fractional dierential equation.
The proposed model has ve parameters: the inltration coecient C, the maximum reservoir size h max , the evapotranspiration depth threshold h min and the two parameters (α, τ 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 signicantly 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 inow 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 dierent 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, rst-order dierential equations (19c) proves decisive in terms of computational eciency.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 specic volume and the outowing 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 eects.The shape of the hysteresis paths in the (h, q) space are radically dierent in the Dud Koshi and the Durzon models, a further illustration of the exibility 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 specied explicitly and does not incorporate as much exibility 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 eects induced by the innite 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 owpath 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 identied as highly desirable in both karst and mountainous catchment modelling (Grae 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 eect 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 dierent short-term behaviour, studying theoretically and experimentally the eect of placing the proposed transfer function in cascade with nite or other innite characteristic time transfer functions, investigating whether incorporating a non-linear behaviour as suggested in Section 3.2.2can 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.
The Riemann-Liouville approach (Liouville, 1832;Riemann, 1953) considers the αth-order derivative as the nth-order derivative of the (n − α)th-order integral: while the Caputo derivative (Caputo, 1967(Caputo, , 1969) ) considers the αth-order derivative as the (n − α)thorder integral of the nth-order derivative: (the lower integration bound a is set to zero in this paper).The particular case n = 1 yields with the advantage that this integral converges for α < 1 provided that f is dierentiable over the integration interval.
The Grünwald-Letnikov denition is (Grünwald, 1867) A last possibility is to dene the derivative from the Laplace transform: where • and • denote respectively the Laplace and inverse Laplace transforms of a function.
Fractional derivatives dier from integer-order derivatives by a number of important features: the fractional derivative is a non-local operator.The derivative a D α f (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, dierent denitions may yield dierent 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 specic conditions.The requirements for commutation are not the same for the Riemann-Liouville and the Caputo denition (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 dierent 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 ≤ t ≤ 10τ 0 , while the bottom graphs illustrate the long-term behaviour over 6 orders of magnitude, for τ 0 ≤ t ≤ 10 6 τ 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 τ 0 .Beyond this range, oscillations appear in the interpolated proles (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 τ 0 after calibration is larger than 4 days, which gives 10 5 τ 0 > 10 3 years.This is far longer than any realistic rainfall-runo simulation.It can be concluded that the approximation (17) of the kernel (11) is extremely satisfactory for rainfall-runo modelling purposes.The number of functions needed (hence the number of sub-reservoirs in the model), R = 15, makes the CPU eort moderate.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/τ 0 for the proposed convolution (11).w * = w (t) /w (0) is the convolution kernel normalized by its value at t = 0.              3. The dimensionless plot coordinates are dened as t * = t/τ 0 and w * = w (t) /w (0).w * = w (t) /w (0) is the convolution kernel normalized by its value at t = 0.                3. The dimensionless plot coordinates are dened as t * = t/τ 0 and w * = w (t) /w (0).

Figures
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 inow model used in the application examples is very simple, especially in the way evapotranspiration is considered.It is acknowledged that much more sophisticated models have been proposed in the literature(Fazal et al., 2005;Hartmann et al., 2013a; Juki c and Deni c-Juki c, 2009).

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 ows 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 runo, snow melt and ice melt, with glacier covering about 14% of the catchment area.The rst 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 uxes 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.
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

Figure 2 .
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/τ 0 for the proposed convolution (11).w * = w (t) /w (0) is the convolution kernel normalized by its value at t = 0.

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

Figure 15 .
Figure15.Approximation of the convolution kernel with a set of exponential functions for the parameter sets in Table3.The dimensionless plot coordinates are dened as t * = t/τ 0 and w * = w (t) /w (0).

Figure 2 :
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/τ 0 for the proposed convolution (11).w * = w (t) /w (0) is the convolution kernel normalized by its value at t = 0.

Figure 3 :
Figure 3: Partitioning an innite characteristic time reservoir into a set of linear sub-reservoirs with nite characteristic times.

Figure 10 :
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 specic volumes stored in the upper and lower compartments of the GR4J model.

Figure 17 :
Figure 17: Approximation of the convolution kernel with a set of exponential functions for the parameter sets in Table3.The dimensionless plot coordinates are dened as t * = t/τ 0 and w * = w (t) /w (0).

Table 1 :
Dud Koshi catchment.Simulation results for the single reservoir model.N.A.:Not