MATLAB Application

Implementing the Horton & Philips Models

Prelude

The primary aim of a MATLAB simulation of the theoretical models governing infiltration and ponding time is obtaining the required results with undue efficiency and accuracy. On a more discreet level, the benefit of attending to the algorithm details while writing the program, outweighs any other advantage one might gain from participating in this MATLAB model application.

Being based on relatively the same algorithm for determining infiltration and ponding time under variable rainfall intensity, the modifications needed to transform the Philips model I was originally assigned to Horton’s were minor, and consequently, this report encompasses both models.

After acquainting the reader with the problem to be tackled and its basic solution mechanism, a copy of each of Philips’ and Horton’s mother program is attached in batch format, in which all subfunctions are clearly highlighted and indexed in detail at the end of the report. Furthermore, to demonstrate the flexible user-MATLAB interface, the reader is then guided explicitly throughout the solution of a sample problem previously solved manually by Horton in class. Finally, the Philips model is used to solve a more advanced problem, and the results are tabulated and plotted, once using MATLAB’s graphing facility and then on a more customized Excel spreadsheet.

(Top)

Algorithm

The problem considered is: given a rainfall hyetograph defined using the corresponding constant time-pulse data representation, and the parameters of Philips’ or Horton’s infiltration equations, determine the ponding time, infiltration rate, cumulative infiltration, and the excess rainfall hyetograph.

Trying to make the program as user-friendly and flexible as possible, minimal effort is required from the user when preparing and typing the input data; for example, only the pulse duration and the final reading are required for the program to automatically provide the user with the corresponding time row-vector, by this also assuring constant time-pulse duration. Although this program does not support variable pulse duration, the user has the option of entering rainfall data in either rate or cumulative format. Moreover, the user is reminded of the units of the variables he/she is prompted to input and warned whenever the number of rainfall entries does not conform to that of the time readings.

To start with, below is a schematic overview of the algorithm employed for both models, Philips’ and Horton’s:

“Applied Hydrology” by Chow, Maidment & Mays

(Top)

Horton Model

Below is a copy of the Horton Model program in batch format:

(Top)

Here the user is prompted for input data entry. Minimal effort in entering the data is required; the time matrix, for example, is automatically generated as soon as the user enters the pulse duration and final time reading. This also secures a constant time interval.

Enter K in t^-1: 1

Enter fo in L/t: 1.2

Enter fc in L/t: 0.2

Enter pulse duration (constant value): 0.5

Enter time for last reading: 1

t =

0 0.5000 1.0000

n =

3

After being given the choice of entering the rainfall date in either cumulative or rate form, the number of entries is checked to match that of the previously entered time data.

Press (1) if rainfall data is cumulative or (2) if it is in rate form: 1

Enter cumulative rainfall data (L) in row vector form: [0 0.4 1.5]

P =

0 0.4000 1.5000

dP =

0 0.4000 1.1000

i =

0.8000 2.2000 0

F =

0 0 0

f =

1.2000 0 0

This being an explicit version of the program I’ve included in the previous page (most “;” have been omitted), the user can trace the loops and conditional statements through those intermediate task executions.

dF =

0 0.4000 0.3754

Fx =

0 0.4000

fx =

0 0.9000

F =

0 0.4000 0

f =

1.2000 0.9000 0

F =

0 0.4000 0.7754

f =

1.2000 0.9000 0.6246

f =

1.2000 0.9000 0.6246

F =

Once all vector entries have been evaluated successively, the sought after results are listed and found to match those obtained in class.

0 0.4000 0.7754

dF =

0 0.4000 0.3754

dPe =

0 0 0.7246

ie =

0 1.4491 0

(Top)

Philips Model

Below is a copy of the Philips Model program in batch format:

(Top)

Based on the parameters and rainfall data given in problem 5.4.6, Philips’ model is used to solve the problem and then have the results tabulated and plotted.

(Top)

n	t	P	dP	i	f	F	dF	dPe	ie
1	0	0	0	25	∞	0	0	0	0
2	1	25	25	45	78.541	25	25	0	0
3	2	70	45	45	45	70	45	0	4.2893
4	3	115	45	25	37.6777	110.7107	40.7107	4.2893	0
5	4	140	25	20	35.2504	135.7107	25	0	0
6	5	160	20	20	33.8316	155.7107	20	0	0
7	6	180	20	0	32.7108	175.7101	20	0	0

(Top)

Subfunctions

Below is the alphabetically sorted list of user-defined functions used in the Horton & Philips mother programs.

function Ft=cumulative(k,s,t,F,f)

Ft=F+k*t-s^2/(2*(f-k))+s*(t+s^2/(4*(f-k)^2))^0.5;

q CUMULATIVE:

function Pd=depth(P)

n=length(P);

Pd(1)=0;

for j=2:1:n

Pd(j)=P(j)-P(j-1);

end

q DEPTH:

function ie=eintensity(dt,Pe)

n=length(Pe);

ie(n)=0;

for j=1:1:n-1

ie(j)=Pe(j+1)/dt;

end

q EINTENSITY:

function Fp=Fponding(k,s,i)

Fp=((s^2)*(i-k/2))/(2*(i-k)^2);

q FPONDING:

function Ft=hcumulative(k,F,fc,dt,ft)

Ft=F+fc*dt+(ft-fc)*(1-exp(-1*k*dt))/k;

q HCUMULATIVE:

function Fp=hFponding(k,fo,i,fc)

Fp=1/k*(fo-i+fc*log((fo-fc)/(i-fc)));

q HFPONDING:

function ft=hrate(k,f,Ft,F,fc,dt)

ft=f-k*(Ft-F-fc*dt);

q HRATE:

function i=intensity(t,P)

n=length(P);

i(n)=0;

for j=1:1:n-1

i(j)=(P(j+1)-P(j))/t;

end

q INTENSITY:

function ft=rate(k,s,F)

ft=k+s*k/((s^2+4*k*F)^0.5-s);

q RATE:

function tp=tponding(Fp,F,i)

tp=(Fp-F)/i;

q TPONDING:

(Top)