A robust rainfall-runoff transfer model

WATER RESOURCES RESEARCH, VOL. 37, NO. 12, PAGES 3207-3216, DECEMBER

A robust

rainfall-runoff

transfer

2001

model

G. Capkun Departmentof Mathematicsand Departmentof Rural Engineering(IATE/HYDRAM) SwissFederal Institute of Technology,Lausanne,Switzerland

A. C. Davison Departmentof Mathematics,SwissFederal Instituteof Technology,Lausanne,Switzerland

A. Musy Departmentof Rural Engineering(IATE/HYDRAM), SwissFederal Instituteof Technology Lausanne, Switzerland

Abstract. This paper describesa simpleand flexiblenew approachto rainfall-runoff transfermodeling,whichjointly modelsthe mean and the variancestructureof runoff. Its meanis taken to be a linear autoregressive combinationof presentand previousrainfall and previousrunoff, while its variancealso dependson rainfall history.Inference for our model is performedusingclassicallikelihoodmethods,and alsoby the more robust techniqueof quasi-likelihood,presupposing no particulardistributionfor runoff. The model hasbeen developedon the Biorde catchmentin Switzerlandand calibratedand validatedfor 12 smallSwisscatchments,with promisingresults. ence. Delleur [1991] givesan overviewof so-calledsystemtheoretictransferfunctionmodels,whichattemptto establish Hydrologicalmodelingsplitsbroadlyinto physicallybased a causallinkagebetweenseveralphenomenawithout describand empiricalor system-based approaches. In the first the aim ing the physicalprocesses involved.The mostcommonmodels, is to understandunderlyingphenomena,with the goal of ex- suchas deterministicautoregressive movingaverage,applied plaininghowthingshappenandwhyprocesses behaveasthey unit hydrographtransferfunction, and classicalARMA moddo. The physicalreality underlyinghydrologicalsystemsis so els, are specialcasesof a generalsingleinput-singleoutput complex,however,that suchmodelsare difficultto validateon model of the form the scantdata usuallyavailable,and this compromises their •(b) C(b) practicalusefulness. The systems approachis concernedwith A(D)Yt= b*F-•-•Xt + D(D) et, the operationof the systemunder considerationrather than with its natureor the physicallawsgoverningit and is intended where X,, Y,, and e, are the input, output, and white noise to bypassphysicalmodelingby establishing a relationshipthat sequences at time t;, is a delaytime betweeninput andoutput; canbe usedto reconstructpastevents,to predictfuture ones, andA, B, C, D, and F are polynomialsin the backwardshift or to generatedata like thosealreadyobserved.The choiceof operatorb. In all thesemodels,X,, Y,, and e, are stationary approachdependson the requirementsof the user,the prob- processes. This approachhas the major drawbackthat stanlem, the availabilityof data, and computationalfacilitiesand dard linear time series models are intended to model variables economicconsiderations. Often, in practice,rainfallandrunoff lying on the real line; in fitting routinesthe white noise is measurements are the only availabledata, and a systemsap- generallytaken to be Gaussian,leadingto symmetricdistribuproachmustbe chosen.The goalof thispaperis to proposea tions for Yt. Rainfall and runoff are nonnegativeand highly however,sotransformations mustbe appliedtoX, methodologyfor solvingthe "inverseproblem,"i.e., the iden- asymmetric, tificationof the systemfunctionfrom giveninput and output. and Y,. This makesthe modelparametersdifficultto interpret that nonIn the termsusedby Singh[1988]and Chow [1988]our pro- evenif their numberis relativelysmalland suggests posalis a symbolic,mathematical,theoretical,nonlinear,time- Gaussianmodels may be more useful. Hipel and McLeod [1994]givereferencesto non-Gaussian modelingof time series varying,lumpedstochasticmodel. in hydrology. During the last3 decades,severaltypesof stochastic models In thispaperwepropose a simpleandflexiblenewapproach for hydrologicaldata have been proposed[Salasand Smith, 1981]. Delleur [1991] usesa conceptualwatershedmodel to to rainfall-runofftransfermodeling,whichneedsonly rainfall and runoff series. The mean of current runoff is treated as an suggestthat the groundwaterstorageand runoff processes autoregression on presentand previousrainfall and previous belongto the classof autoregressive movingaverage(ARMA) runoff, while its variancealsodependson rainfall history.Our processes. The order of the ARMA processfitted dependson model, which is linear in its mean structure but not in the usual the type of processanalyzedand is in practicechosenby stahydrological sense[Singh,1988;Singh,1995],is detailedin tisticalgoodness-of-fit criteria,judgment,or personalprefer1.

Introduction

section 2. Inference is described in section 3, where we show

Copyright2001 by the American GeophysicalUnion. Paper number2001WR000295. 0043-1397/01/2001WR000295 $09.00

howit canbe fitted asa generalizedlinearmodel(GLM) with runoff treatedas a gammaresponse,usingclassicallikelihood theoryand martingalelimit resultsasa basisfor inference.We 3207

3208

CAPKUN

ET AL.: RAINFALL-RUNOFF

TRANSFER

MODEL

then put our approachinto the more generalcontextof quasi- expectall the/3i to be positive,thoughwe do not imposethis likelihood theory and estimatingequations,assumingno par- when fitting the model. ticular distributionfor runoff, and showthat resultsobtained Primary interest attachesto how the mean responseis afusingthesetwo techniquesare very similarfor our model. In fectedby covariates, but we havealsofoundit usefulto model section 4 we outline its fit to data from the Biorde catchment changesin crt2. In the usuallinearmodel,thereis no relation small Swiss

between themeanandvariance of theresponse, andthisisalso

catchments, with promisingresults.The final sectioncontainsa

the casefor conventionalARMA models.It is rarely true for positivedata,however,for whichthe varianceis oftenproportional to the mean squared,as for the exponential,gamma, lognormal,and Weibull distributions.In our caseit is usefulto accommodate the possibilitythat the currentvarianceof runoff dependson recent rainfall historyand to supposethat crt2 __

in Switzerland

and then extend

it to 12 further

brief discfission.

2.

Model

Let Xt and Yt representhourly rainfall and runoff at times wheretheprecision parameter vt = t = 1, ..., T. These measurementsare nonnegative,depen- #,t2/vt, dependson the indicator dent, and positively skewed time series.Although rainfall showsno stronglysuggestivestochasticstructure,the runoff seriesis not time-reversible.Our aim is to find an appropriate , 1, otherwise,

itt• ={0,Xt_ 1..... Xt_ k'--0,

form of the transfer function that transforms rainfall into run-

off. When studyingthe rainfall-runoffrelationship,the usual practiceis to replacemeasuredrainfall with effectiverainfall, i.e., the portion of rainfall that contributesto direct runoff, includingsurfaceandpromptsubsurface runoff.Variousmethods for obtainingeffectiverainfall have been studied[Park et al., 1999; Musy, 1998], but to incorporatea processwhose model cannot be validated could introduce seriouserrors, and

exp(a1 + ot2It,•: )

of rainfall in the precedingk days.This impliesthat the con-

ditionalcoefficientof variationof runoff,var(YtlDr)l/2/ E(YtIDt) = v;-1/2, depends onlyon recentrainfall. In the data analysisbelow,we saythat the periodpreceding

timet is rainless if xt_i < 0.1 mmfor allj = 1, ..., k or if xt - 1 ..... xt - • = 0 for ANETZ stationstreatingallxt with less than 0.1 mm of rain as zero. Otherwisethe period is

we prefer to use observedrainfall. As usualwith stochastic models,the key idea is to treat the availabledata as the outcome of some random experiment.As our goal is to model runoff as a functionof observedrainfall,we supposethat the runoffseries{ Yt) is randomandtreat the rainfallseries{Xt }

period it showsa slowessentiallydeterministicdeclinetoward zero. ANETZ is a Swissautomaticmeasuringnetworkof climate variablesrecordedevery 10 min. Our parametervector

as known.

0-- (/30,ooo,/3k,•/1, øøø,•/1,Otl,Or2) T contains termsforthe

considered to be rainy,with It,•: = l. Duringa rainyperiod, runoff increasesand can be rather variable, but in a rainless

dependenceof E(Yt) on previousrunoff and rainfall and for at time t, i.e., dependenceof its varianceon the presenceof recentrainfall. Our proposalis a type of generalizedlinear model (GLM) D t-- (Xt,..., Xt-k,... ;Yt-1,..., Yt-t.... ). [McCullaghand Nelder,1989].The linear predictorr/t is equal to the mean, so the link functionis the identity,r/t = tzt. The This is the informationavailableto the modelerfor prediction variance crt2 depends on tzt throughthe variancefunction of the runoff Yt at time t. In practice,we use a smallsubsetof V(p,t) = tzt 2, givingcrt2 = V(ld, t)• t = Id, t2/l,•t , whereqbt= 1,•1 Dt, typicallytakingk and l to be smallintegersand ignoring is the dispersionparameterof the GLM. The complicationis the earlier observationsXt_k_ 1, Xt-•,-2, ... and Yt-t-1, that our model is a nonstationarytime series,and so standard Yt-t-2,.... The assumption that Yt depends only on resultson inferencefor GLMs do not apply.We discussinferYt- 1, ..., Yt-t amountsto the/th orderMarkov propertyfor ence in the next section. { Yt ). We supposethat conditionalon Dt, Yt comesfrom some Let D t denote the presentand past rainfall and past runoff

densityh(Ytldt; 0), whereYt anddt represent valuesof the .randomvariablesYt and D t and the vectorparameter0 sum3. Statistical Inference marizestheir dependencestructure. We specifythe modelin termsof the dependence of the first 3.1. Parameter Estimation two momentsof Yt on covariates,the mean and the variance Likelihoodfunctionsare centralto parametricstatisticalinbeing definedas ference.They are basedon probabilitymodelsfor the observed #,t-- E(YtIDt), o't 2'- var(YtIDt). response,in our case runoff, regarded as functions of the parametersof the underlyingdistribution.The densityof the The mean runoff at time t is supposedto be a linear combi- time seriesof runoffy 1, -.., Yr conditionalon the time series nation of presentand previousrainfall and previousrunoff, of rainfallxl, ..-, x r can be written as k

l

Id't--E [3iXt-i '-[-E •j Yt-j i=0

f(yi ....

, yrxl, ßßß, xr; O)

t=r+l,...,T, T

j=l

wherer = max(k, l ) andthe valuesof k andI are determined by a model choiceproceduredetailed in section3.3. The di-

= f(Yl..... yrX•,... , Xr;O) I-I f(Ytldt;0),

(1)

t=r+l

mensionless parameters Ti measurethe association between where r = max(k, l). The joint distributionof Y1, ..., Yr expectedcurrentandpastrunoff,while the parameters/3ihave the dimensions of runoff/rainfall

and measure the effect of

rainfall on expectedcurrentrunoff, assumedto be linear. We

cannotbe determinedwithout additionalassumptions. A simple alternative is to base the inference only on the second component of (1) which is the conditional likelihood of

CAPKUN ET AL.: RAINFALL-RUNOFF

Yr+l,''',

TRANSFER MODEL

3209

estimation Yr givenY1, ..., Yr. The corresponding condi- by solving(3). This is knownas quasi-likelihood [Wedderburn, 1974];it is closelyrelatedto the theoryof estimatingfunctions (seeDesmond [1991]for furtherreferences).

tional log likelihoodcanbe written T

More detailsare givenin section3.4.

t=r+l

3.2.

Parameter Uncertainty

whereet(O) = logf(Ytldt;O) is theloglikelihood contribu-

The standardapproachto likelihoodconfidenceintervals tion from the tth runoff conditionedon previousrunoff and usesthe fact that under regularityconditionsthe maximum currentand previousrainfalldt. As r is generallytiny com- likelihood estimator • hasapproximately a normal distribution

paredwithT, the information andefficiency lossdueto ignor- withmean0 andvariance-covariance matrixJ(0 ) - 1, i.e., ingthe firstcomponent on the rightof (1) is verysmall. •}'