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);

[Maple Math]
[Maple Math]
[Maple Math]

Warning, new definition for norm

Warning, new definition for trace

[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]

> N:=9;M:=4;

[Maple Math]

[Maple Math]

Define time interval duration in seconds.

> DT:=6*3600;

[Maple Math]

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

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

[Maple Math]

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

[Maple Math]

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]]);

[Maple Math]

Obtain the transpose of [P].

> PT:=transpose(P);

[Maple Math]

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

> PTP:=multiply(PT,P);

[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]

Obtain the matrix inverse of [PTP].

> PTPINV:=inverse(PTP);

[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]

Obtain the so-called pseudo-inverse matrix:

> Z:=multiply(PTPINV,PT);

[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]

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

> U:=multiply(Z,Q);

[Maple Math]

Check for volume of unit hydrograph.

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

[Maple Math]

Make predictions of the discharge hydrograph using the estimated Unit Hydrograph

> QHAT:=multiply(P,U);

[Maple Math]

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

[Maple Math]
[Maple Math]

Define arrays of decision variables.

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

[Maple Math]

[Maple Math]

[Maple Math]

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

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

[Maple Math]

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};

[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]
[Maple Math]

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 .

[Maple Math]
[Maple Math]

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);

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

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;

[Maple Math]

[Maple Math]

[Maple Math]

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;

[Maple Math]

[Maple Math]

[Maple Math]

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))]);

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

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))]);

[Maple Math]

[Maple Math]

[Maple Math]

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

[Maple Math]

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

[Maple Math]

>