Unite Hydrograph Estimation
Nash Model
Cascade of n Linear Reservoirs
The following procedure assumes that basin response can be described using either the Nash Model (cascade of n linear reservoirs of equal storage constant k ) so that the instantaneous unit hydrograph can be writen as:
> u:=(t)->1/k*(t/k)^(n-1)*exp(-t/k)/GAMMA(n);
Linear Channel-Linear Reservoir Model
or using the linear channel-linear reservoir model so that the instantaneous unit hydrograph can be writen as:
> u:=piecewise(t>T,1/k*exp(-(t-T)/k));
> with(student);with(linalg);
> N:=6; M:=4;
Define discharge hydrograph matrix in units of m3/s.
> Q:=matrix(N,1,[30,250,500,400,180,30]);
Define rainfall hyetograph matrix in units of cm/h.
> H1:=matrix(M,1,[0.2,1.5,1.0,0.5]);
Define time interval, DT, in seconds.
> DT:=2700;
Compute area of the basin by determining the area A that would produce an input volume (effective rainfall volume) equal to the observed outflow volume (direct runoff volume). Area is in units of m2.
.
> A:=sum('Q[k,1]','k'=1..N)*DT/(sum('H1[k,1]','k'=1..M)*DT/100/3600);
Transform the rainfall hyetograph matrix from units of intensity (cm/h) to units of flux (m3/s). All future computations are in consistent units in the SI system.
> H:=scalarmul(H1,A/3600/100);
Compute the first and second normalized moments of the input rainfall hyetograph and the output discharge hydrograph. Observe that both input and output are expressed in units of m3/s.
First moments, M1(Q) and M1(H) (in units of time (seconds)),
> M1Q:=DT*(Q[1,1]/2*DT/2+sum((Q[l-1,1]+Q[l,1])/2*(l-0.5)*DT,l=2..N)+Q[N,1]/2*(N+.5)*DT)/(sum('Q[l,1]','l'=1..N)*DT);
> M1H:=DT*(sum(H[l,1]*(l-0.5)*DT,l=1..M))/(sum('H[l,1]','l'=1..M)*DT);
Second moments, M2(Q) and M2(H) (in units of time squared),
> M2Q:=DT*(Q[1,1]/2*((DT/2)^2+DT*DT/12)+sum((Q[l-1,1]+Q[l,1])/2*(((l-0.5)*DT)^2+DT*DT/12),l=2..N)+Q[N,1]/2*(((N+.5)*DT)^2+DT*DT/12))/(sum('Q[l,1]','l'=1..N)*DT);
> M2H:=DT*(sum(H[l,1]*(((l-0.5)*DT)^2+DT*DT/12),l=1..M))/(sum('H[l,1]','l'=1..M)*DT);
Use the method of moments to obtain estimates of the parameters of the Nash Model, n and k (Cascade of n equal linear reservoirs with parameter k ).
> nk:=M1Q-M1H;
> k:=((M2Q-M2H)-2*nk*M1H-nk*nk)/nk;
> n:=nk/k;
Use the method of moments to obtain estimates of the parameters of Linear Channel-Linear Reservoir model, T and k .
> Tpk:=M1Q-M1H;
> kk:=sqrt((M2Q-M2H)-2*Tpk*M1H-Tpk*Tpk);
> T:=Tpk-kk;
In order to use either the Nash model or the Linear Channel-Linear Reservoir model to predict discharges for the given rainfall event, we must obtain the ordinates of the corresponding unit hydrograph. This is accomplished by deriving the corresponding pulse response function, h(t ), for each one of the models. The UH ordinates, U , are then obtained as
U[n-m+1]=h[(n-m+1)dt], where dt is the duration of each rainfall pulse.
First obtain the UH ordinates for the Cascade of Linear Reservoirs model (Nash model). In the analysis that follows, the value of n obtained above representing the number of linear reservoirs in our model has been rounded to the nearest integer to be consistent with the modeling assumption. Integer
> n:=round(n);
> hterm:=(m)->int((t/k)^(n-1)*exp(-t/k)/GAMMA(n)/k, t=(m-1)*deltat...m*deltat)/deltat;
> deltat:=1*DT;
> UCLR:=matrix(N-M+1,1,[evalf(hterm(1)),evalf(hterm(2)),evalf(hterm(3))]);
Second, obtain the UH ordinates for the Linear Channel-Linear Reservoir model. Observe that u(t) for this model is equal to zero for 0<t<T , and that this must be taken into account when evaluating U[n-m+1]=h[(n-m+1)dt].
> h1term:=(T)->int(exp(-(t-T)/kk)/kk, t=T...deltat)/deltat;
> hhterm:=(m,T)->int(exp(-(t-T)/kk)/kk, t=(m-1)*deltat...m*deltat)/deltat;
> ULCLR:=matrix(N-M+1,1,[evalf(h1term(T)),evalf(hhterm(2,T)),evalf(hhterm(3,T))]);
> sum('UCLR[m,1]','m'=1..N-M+1)*DT;
> sum('ULCLR[m,1]','m'=1..N-M+1)*DT;
Define P matrix to be used in the discrete convolution equation.
> P:= matrix(6,3,[[H[1,1]*DT,0,0],[H[2,1]*DT,H[1,1]*DT,0],[H[3,1]*DT,H[2,1]*DT,H[1,1]*DT],[H[4,1]*DT,H[3,1]*DT,H[2,1]*DT],[0,H[4,1]*DT,H[3,1]*DT],[0,0,H[4,1]*DT]]);
Compute prediction of the observed hydrograph using the UH obtained from the Cascade of Linear Reservoirs model.
> QHATCLR:=multiply(P,UCLR);
Compute prediction of the observed hydrograph using the UH obtained from the Linear Channel-Linear Reservoir model.
> QHATLCLR:=multiply(P, ULCLR);
Define the new matrix P to be used in the discrete convolution equation.
> PNEW:=matrix(7,3,[[0.003*A,0,0],[0.007*A,0.003*A,0],[0.015*A,0.007*A,0.003*A],[0.003*A,0.015*A,0.007*A],[0.001*A,0.003*A,0.015*A],[0,0.001*A,0.003*A],[0,0,0.001*A]]);
Compute prediction of the hydrograph that would be observed for the given rainfall hyetograph, using the UH of the Cascade of Linear Reservoirs model.
> QHATLCLR:=multiply(PNEW,ULCLR);
Compute prediction of the hydrograph that would be observed for the given rainfall hyetograph, using the UH of the Linear Channel-Linear Reservoir model.
> QHATCLR:=multiply(PNEW,UCLR);
>