CE522 Engineering Hydrology

Jorge A. Ramirez

Continuous Muskingum Routing Model

> restart;with(student);with(inttrans);with(DEtools);

Continuity Equation:

> diff(S(t),t)=In(t)-Q(t);

Storage - Discharge Relationship:

> S(t)=k*[x*In(t)+(1-x)*Q(t)];

Differentiating the S vs. Q relationship and substituting into the continuity equation obtain an ordinary linear differential. This ODE characterizes the spatially-lumped dynamics of flow in a channel reach whose storage discharge relationship is given by the S vs. Q equation above:

> ODE:=diff(k*(1-x)*Q(t)+k*x*In(t),t)=In(t)-Q(t);

Use the Laplace transform technique to solve the resulting ODE. This technique consists in taking the Laplace transform of the ODE in order to convert it into an algebraic equation.

> laplace(ODE,t,s);

Observe that this is a simple algebraic equation in which the unknown quantity is the laplace transform of Q(t). Solving for laplace(Q(t)) leads to:

laplace(Q(t),t,s):= (k*(1-x)*Q(0)+k*x*In(0))/(k*(1-x)*s+1)+laplace(In(t),t,s)/(k*(1-x)*s+1)-k*x*s*laplace(In(t),t,s)/(k*(1-x)*s+1);

Finally, taking the inverse laplace transform leads to the desired solution for Q(t):

> QQ:=invlaplace(laplace(Q(t),t,s),s,t);

which can be re-written as:

> Q:=t->(k*(1-x)*Q(0)+k*x*In(0))*exp(-t/k/(1-x))/k/(1-x)-k*x*In(t)/(k*(1-x))+int(In(tau)*exp(-(t-tau)/k/(1-x))/k/(1-x)/(1-x),tau = 0 .. t);

which is the solution presented in class.

In order to obtain the impulse response function u(t) for this system, that is, the instantaneous unit hydrograph, we let the initial conditions be equal to zero, and the forcing function be equal to an instantaneous impulse at time t=0 -a Dirac delta function. This leads to the following solution:

> Q(0):=0;In(0):=0;In:=t->Dirac(t);

> u:=t->-x*Dirac(t)/(1-x)+exp(-t/k/(1-x))/(k*(1-x)*(1-x));

The parameters of the Muskingum Routing equation are k and x . The parameter k represents an average travel time of a perturbation through the river reach; and x represents the relative importance of the input in defining the S vs. Q relationship. These parameters can be estimated using the Method of Moments as presented in class. In order to do so we need to obtain the First and Second moments of the IUH, as follows:

First Moment of u(t) about the origin:

> M1:=Int(t*u(t),t=0..infinity);

Second Moment of u(t) about the origin:

> M2:=Int(t^2*u(t),t=0..infinity);

For positive values of k and for 0< x <1, the ratio 1/ k /(1- x ) > 0. Thus,

> assume(1/k/(1-x)>0);

Finally, evaluating the above integrals leads to:

First Moment:

> M1:=simplify(int(t*u(t),t=0..infinity));

and Second Moment:

> M2:=simplify(int(t^2*u(t),t=0..infinity));