CE522 Engineering Hydrology

Jorge A. Ramirez

Unit Hydrograph Estimation

Least Squares Procedure

The following procedure produces estimates of the unit hydrograph (UH) by minimizing the sum of the squares of the estimation errors. Estimation errors are defined as the difference between the observation and the prediction, (Q observed - Q predicted ).

The discrete convolution equation has been rewritten in matrix form as:

[Q]=[P][U]

where [Q] is the column vector containing the N direct runoff observations; [P] is the precipitation matrix containing the volumes of M rainfall pulses; and [U] is the column vector containing the unit hydrograph (UH) ordinates. In what follows, N is the number of direct runoff ordinates, M is the number of effective rainfall pulses, and DT is the duration of the observation intervals. Thus, the instantaneous direct runoff is observed at every DT, and the rainfall pulses are volumes of rainfall of constant intensity over DT.

The solution presented below corresponds to the following effective rainfall hyetograph:

Time(h) Average Intensity (m3/s)

0 - 6 100

6 - 12 300

12 - 18 200

18 - 24 100

The corresponding direct runoff discharge hydrograph is:

Time(h) Discharge (m3/s)

6 10

12 70

18 165

24 180

30 142

36 79

42 38

48 13

54 3

Initialize linear algebra library.

> with(student); with(linalg);

Warning, new definition for norm

Warning, new definition for trace

> N:=9;M:=4;

Define time interval duration in seconds.

> DT:=6*3600;

Define the direct runoff matrix [Q obs ] of dimensions (Nx1).

> Q:=matrix(N,1,[10,70,165,180,142,79,38,13,3]);

> H:=matrix(M,1,[100,300,200,100]);

Construct matrix [P] of dimensions Nx(N-M+1).

> P:= matrix(N,6,[[100*DT,0,0,0,0,0],[300*DT,100*DT,0,0,0,0],[200*DT,300*DT,100*DT,0,0,0],[100*DT,200*DT,300*DT,100*DT,0,0],[0,100*DT,200*DT,300*DT,100*DT,0],[0,0,100*DT,200*DT,300*DT,100*DT],[0,0,0,100*DT,200*DT,300*DT],[0,0,0,0,100*DT,200*DT],[0,0,0,0,0,100*DT]]);

Obtain the transpose of [P].

> PT:=transpose(P);

Multiply the [transpose of P] and [P].

> PTP:=multiply(PT,P);

Obtain the matrix inverse of [PTP].

> PTPINV:=inverse(PTP);

Obtain the so-called pseudo-inverse matrix:

> Z:=multiply(PTPINV,PT);

Obtain the Unit Hydrograph matrix [U] of dimensions (N-M+1)x1.

> U:=multiply(Z,Q);

Check for volume of unit hydrograph.

> UHvolume:=sum('U[k,1]','k'=1..N-M+1)*DT;

Make predictions of the discharge hydrograph using the estimated Unit Hydrograph

> QHAT:=multiply(P,U);

Linear Algebra Procedure

The following procedure produces estimates of the unit hydrograph (UH) by minimizing the sum of the absolute value of the estimation errors. Estimation errors are defined as the difference between the observation and the prediction, (Q observed - Q predicted ).

The discrete convolution equation can be rewritten in matrix form as:

[Q]=[P][U]

where [Q] is the column vector containing the N direct runoff observations; [P] is the precipitation matrix containing the volumes of M rainfall pulses; and [U] is the column vector containing the unit hydrograph (UH) ordinates. In what follows, N is the number of direct runoff ordinates, M is the number of effective rainfall pulses, and DT is the duration of the observation intervals. Thus, the instantaneous direct runoff is observed at every DT, and the rainfall pulses are volumes of rainfall of constant intensity over DT.

This problem reduces to a constrained minimization problem in which the function to be minimized and the constraints are linear. This type of problems can be solved using linear programming methods. In the illustration below the Simplex algorithm is used.

Initialize the simplex algorithm library.

> with(simplex);

Warning, new definition for basis

Warning, new definition for maximize

Warning, new definition for minimize

Warning, new definition for pivot

Define arrays of decision variables.

> u:=vector(6);theta:=vector(9);beta:=vector(9);

Define the objective function as S ( q n + b n ).

> obj:=sum(theta[k]+beta[k],k=1..9);

Define the set of constraints as [P][u]+[ q ]-[ b ]=[Q obs ].

> consts:={P[1,1]*u[1]+theta[1]-beta[1]=Q[1,1], P[2,1]*u[1]+P[1,1]*u[2]+theta[2]-beta[2]=Q[2,1], P[3,1]*u[1]+P[2,1]*u[2]+P[1,1]*u[3]+theta[3]-beta[3]=Q[3,1], P[4,1]*u[1]+P[3,1]*u[2]+P[2,1]*u[3]+P[1,1]*u[4]+theta[4]-beta[4]=Q[4,1], P[4,1]*u[2]+P[3,1]*u[3]+P[2,1]*u[4]+P[1,1]*u[5]+theta[5]-beta[5]=Q[5,1], P[4,1]*u[3]+P[3,1]*u[4]+P[2,1]*u[5]+P[1,1]*u[6]+theta[6]-beta[6]=Q[6,1], P[4,1]*u[4]+P[3,1]*u[5]+P[2,1]*u[6]+theta[7]-beta[7]=Q[7,1], P[4,1]*u[5]+P[3,1]*u[6]+theta[8]-beta[8]=Q[8,1], P[4,1]*u[6]+theta[9]-beta[9]=Q[9,1], DT*sum(u[l],l=1..6)=1};

Carry out the minimization problem and obtain the UH ordinates [u].

> minimize(obj, consts, NONNEGATIVE);

The unit hydrograph ordinates u n are to be compared to those from the Least Squares procedure, U n .

Nash Model

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 homogeneous units (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;

Linear Channel-Linear Reservoir Model

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].

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 equal to and representing the number of linear reservoirs in our model has been rounded to the nearest integer to be consistent with the modeling assumption.

> n:=3;

> 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(6,1,[evalf(hterm(1)),evalf(hterm(2)),evalf(hterm(3)),evalf(hterm(4)),evalf(hterm(5)),evalf(hterm(6))]);

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(6,1,[evalf(h1term(T)),evalf(hhterm(2,T)),evalf(hhterm(3,T)),evalf(hhterm(4,T)),evalf(hhterm(5,T)),evalf(hhterm(6,T))]);

> sum('UCLR[m,1]','m'=1..N-M+1)*DT;

> sum('ULCLR[m,1]','m'=1..N-M+1)*DT;

>