nonlinearity and nonstationarity

in RAINFALL-RUNOFF RELATIONS FOR KARSTIC SPRINGS

 

David LABAT ^(1,2), Rachid ABABOU ⁽¹⁾, Alain MANGIN ⁽²⁾

⁽¹⁾ Institut de Mécanique des Fluides de Toulouse, ⁽²⁾ Laboratoire Souterrain de Moulis

 

Corresponding authors

D. LABAT, R. ABABOU, Institut de Mécanique des Fluides de Toulouse

Allée du Professeur Camille Soula, 31400 Toulouse, FRANCE

Tel : +33 (0)5 61 28 58 73 ; Fax : +33 (0)5 61 28 58 99

e-mail : labat@imft.fr , ababou@imft.fr

 

 

Abstract

Karstic aquifers are distributed all around the world (USA, China, Europe, ...). They provide large underground reserves of water of high quality but are quite difficult to exploit since their behaviour is not well-understood (floods and low-flow). A temporal and/or spatial model would have applications not only in water management (pollution and pumping) but also in flood prevention.

A karstic watershed is a three-dimensional basin which involves both surface and subsurface flows including infiltration, surface runoff and subsurface processes through fissures, underground streams and reservoirs. Due to the complexity of the physical structure of these systems, we focus on the relations between two temporal signals : an input constituted by the basin-averaged rainfall rate and an output constituted by runoff at the supposed main outlet of the basin. We will study more carefully two karstic basins located in the Pyrénées mountains in France (Aliou and Baget) for which we have long records of daily and semi-hourly rainfall-runoff data.

The proposed input-output causal relations are based on kernels of increasing orders, which relate input to output via nonlinear Volterra convolution expansions. Truncated to the first order, this yields the classical linear convolution model. Higher order terms correspond to nonlinear input/input products. However, we restrict our approach to order two. In the present work, we analyse the limitations of linear models and the performance of the nonlinear ones using a stochastic interpretation of the signals. The latter are able to take into account in a satisfying statistical way the complexity of the temporal behaviour of karstic basins.

Finally, we show that we can also take in account nonstationarity thanks to other methods : the continuous wavelet transform and the orthogonal multiresolution wavelet analysis. The mathematical aspects of those techniques are briefly exposed, and they are both applied to springflow pumping and intermittent springflow data from two other karstic springs (Larzac plateau, Aveyron, France) in order to filter and to identify specific sub-processes.

 

Keywords: Karstic systems, rainfall-runoff, floods, convolution, nonlinear Volterra expansions, continuous wavelet transform, orthogonal multiresolution analysis.

 

karst physical structure

Karstic systems are highly fissured calcareous formations characterised by great heterogeneities at many different scales of observation. The karstification process (dissolution of limestone under the action of meteoritic waters) allows the development of network of conduits and reservoirs. Karstic aquifers can be decomposed in three main parts corresponding to different flow regimes (Mangin, 1975) (Fig. 1) :

Ä     the non karstic impluvium comprises all poorly permeable surfaces which allows a superficial runoff and constitutes the first degree of organisation of the drainage system. It may contain superficial water storage and could be sensitive to evapotranspiration.

Ä     few meters below the surface, a more permeable zone called "unsaturated karst" is characterised by two kinds of infiltration. Vertical fissures and sinkholes allow rapid, vertical single phase flow during intense recharge whereas in the fissure system, water flows in a slower two-phase flow regime (air/water).

Ä     the "saturated karst" constitutes the main storage of water, with a highly organised and hierachised drainage system composed of a main drainage conduit and of annex systems (secondary drains and reservoirs).

 

 

Figure 1 : Karst physical structure (Mangin, 1975)

 

Therefore, karstic media appear as very different from classical porous media. The main difference is the dissociation of transmissive and capacitive functions : transmissive functions are ensured by conduits and drains whereas capacitive functions are ensured by reservoirs of annex systems. Moreover, the physical description of karstic systems leads us to emphasize two kinds of discontinuities : a physical heterogeneity with natural wells, holes and fissures of all sizes, and a dynamic heterogeneity due to the different time responses of parts of the system.

 

study area

We study here two karstic watersheds located in the Pyrénées (Ariège, France) : Aliou and Baget. Those two basins are geographically and geomorphologically similar. Their median altitude is about 940 m, their surface is about 13 km², their specific runoff is both 36 l/s/km² and their mean daily runoff about 0.45 m³/s. We have at our disposal long and ininterrupted chronological series of daily rainfall rate and runoff on the period from 1970 to nowadays. Because of their geographic proximity, rainfall rate is considered as identical. Concerning Aliou basin, we have rainfall-runoff data at a sampling rate of 30 minutes on a three years period. This finer resolution is due to the very short time response of this watershed known for its unexpected swells of runoff : discharge rate can increase from 0.1 m³/s to nearly 30 m³/s in less than 10 hours.

Finally, we also use two other karstic springs located in Larzac (Esperelle spring and Cernon spring) for wavelet analysis focusing on non-stationarity.

 

Linear and nonlinear rainfall-runoff convolution models

Rainfall rates and runoff are considered as two random processes autocorrelated and intercorrelated. We make the hypothesis of a stationary convolution model between those temporal processes. This model is based on Volterra kernels noted h_n(t₁,t₂, .. ,t_n) fullfilling the causality principle. Rainfall rates (P(t)) and runoff (Q(t)) are then related by the general convolution equation :

(1)

This study is restricted to first and second order expansion and (1) reduces to:

1^st order : (2)

2^nd order : (3)

Different methods to determine Volterra kernels are exposed below.

 

Linear model identification

Concerning the first order linear model, we use a statistical method based on orthogonality principle (Papoulis, 1964). The optimal causal kernel h^*(t) which minimises the mean square error between real and simulated runoff is given by another convolution relation (4) linking R_PQ and R_PP (respectively rainfall/runoff and rainfall/rainfall correlations) :

(4)

The discretisation of the integral formulation (4) leads to a linear system solution followed by an optimisation of the memory length of the kernel.

 

nonlinear model identification

In order to identify first and second order kernels corresponding to the nonlinear models, we use two methods.

Ä Volterra/Meixner : this first method (Xia, 1991) is based on the projection of the relation (3) on Meixner orthogonal polynomials M_p(i)=F_p(i-1) where the discrete-time functions F_p(i) are defined by the recursive relation :

(5)

with and , i=1,2, ...

We note {c_p} and {d_q1,q2} the coefficients of the decomposition of the Volterra kernels on the Meixner orthogonal polynomial series, e.g :

and (6)

Using the symmetric property of the second order kernel (Amorocho et al, 1971), the projection of (3) leads to a linear system in term of the coefficients {c_p} and {d_q1,q2} :

(7)

with (m : memory length of Volterra kernels)

In the literature, other polynomial series have also been employed such as Laguerre polynomials (Wiener, 1958), Tchebytchev polynomials (Papazafiriou, 1976), Hermite polynomials (Hino, 1979) or Stiefel-Forsythe polynomials (Papazafiriou, 1976).

Ä Volterra/Nash : this second method (Napiorkowski, 1986) constitutes a more physical approach based on cascade model of N nonlinear reservoirs. The water-stock/discharge (S/Q) relation of each nonlinear reservoir i is of the form :

and , i=1, ... ,N , (8)

where b appears as a nonlinearity coefficient. The first and second order Volterra kernels are obtained by matching the perturbed solution of model (8) to the Volterra model (2)-(3) whence:

(9)

The last step is the optimisation of a, b and N (number of reservoirs), minimising the mean square error between observed and simulated runoff (Diskin, 1984). Note that the two approaches above (Volterra/Meixner and Volterra/Nash) are similar, but the kernels and the solution methods differ.

 

application to pyrenean springs : results and discussions

We apply linear and nonlinear models (Volterra/Meixner) to identify the daily rainfall-runoff relationship of two pyrenean karstic basins. We use a calibration/validation approach. The calibration step consists in identifying the different kernels and/or to optimise them on a given period. The validation step consists in generating a runoff sequence to analyse the goodness-of-fit of the model on a distinct period. The calibration period extends from 25/04/1970 to 01/01/1985 and the validation period extends from 01/01/1985 to 25/02/1997. For the Volterra/Meixner method, we compare the quality of the model using the Nash coefficient F=1-Var(e)/Var(Q) where e is the error Q_DATA-Q_MODEL. Note that represents the normalised r.m.s. error (root mean square error).

Table I sums up the results. The use of the nonlinear model increases significantly the accuracy of the fit for both calibration and validation periods for the two springs. Figure 2 shows the different kernels obtained for Aliou spring at the daily sampling rate with linear and nonlinear models.

Finally, the Volterra/Nash model (Fig. 3) appears quite good at predicting floods at finer sampling rate. The nonlinear part of the model was decisive in achieving this goodness-of-fit.