Journal of Hydrology 12 (1971) 363-376; © North-Holland Publishing Co., Amsterdam N o t to be reproduced by pbotoprint or microfilm without written permission from the publisher
C A L C U L A T I O N OF M E A N AREAL D E P T H OF P R E C I P I T A T I O N
J. E. AKIN Assistant Professor, Department of Engineering Mechanics, The University of Tennessee, Knoxville, Tenn., U.S.A.
Abstract: A new procedure for the calculation of the mean areal depth of precipitation is presented and compared with the commonly used methods. The procedure is conceptually simple and easily programmed for computer application. The analysis of precipitation data for hydrologic studies is, of course, related to the precipitation over a given area, yet the precipitation measurements must be made at specific points. The three c o m m o n methods that have been used to describe the mean areal depth of precipitation in terms of the precipitation measurements at points within the area of interest are: (a) the mean station method, (b) the Thiessen method, and (c) the Isohyetal map method1). An alternate method for the analysis of precipitation data will now be introduced. The development is analogous to some elementary concepts used in finite element analysis techniques2). Consider a region of interest as illustrated in Fig. la. The dashed line could be the boundary of the watershed area, whereas the points indicate the location of the precipitation gages within and near the area. A series of straightlines are arbitrarily choosen to connect the gage points and to form a finite number (m) of triangular sub-areas. These sub-areas are assigned an identifying sequence of integers, beginning with one. (The area numbers are underlined in the figures.) An enlarged view of the nth typical triangular sub-area of Fig. la is shown in Fig. lb. Each point (gage) is also assigned an identifying integer, beginning with one. The integers that define the corners of the typical sub-area (proceeding about the perimeter in a counter-clockwise manner) will be denoted by the letters i, j, and k. The coordinates of the corner points (gages), of the nth sub-area, are (x i, Yi), etc; and the precipitation measurements are Hi, Hi, and Hk. It is assumed that the precipitation over the sub-area varies linearly between the three corner gage points. Thus at any point (x, y) interior to the nth sub-area the precipitation is H n (x, y)
= noel -~- nO~2X -~- n ~ 3 Y
363
(l)
364
J.E.
AKIN
where the e's are constants related to the gage measurements at the corners. Since Eq. (1) is assumed to be valid for all values of x and y associated with the sub-area then it is valid at each corner so that H i = .cq + nO~2Xi -I'- nO~3Yi
(2)
H i = .oq + .o~2xi + nO~3Yi H k = ,cq + no~2x k "~ n~3Yk . 1
',",,
I -H \
\
x
~
4
5 (a)
/7/j
jH]
yil H •r
)
(b)
x
Fig. 1. Notation for triangular sub-areas. Solving the above set of equations for the ~'s yields .oq = [ a i H i + a j H j + a k H k ] / 2 A . .~2 = [ b i l l * + b j H j + b k H k ] / 2 A .
(3)
.ct 3 = [ciHi + c j H j + C k H k ] / Z A .
where a i = xjy k - Xky j b i = y j -- y~ Ci =
X k -- Xj
(4)
365
CALCULATION OF MEAN AREAL DEPTH OF PRECIPITATION
etc. (following a cyclic permutation of i , j , k), and where A. is the sub-area: (5)
A n = [a i + a; + a k ] / 2 .
The differential volume of rainfall at any point within the sub-area is dQ = H (x, y) dA
(6)
so that the total volume of rainfall associated with the sub-area is:
Q. = f f [.o~, + .~zx + .=3y] dx dy
(7)
Q. = A. [.~, + .~2 (x, + x; + xk)13 + .~3 (Y, + Y; + yJ/3].
(8)
which becomes
Therefore, the total volume of rainfall over the triangular mesh is
Q = ~, Qn
(9)
n=l
IO j
I ,8
J
,
,
t3 /
,/\
J
20
.2
t
I
/
I
~ 22
z_~
i'z
. . . .
\/,~
5
I
7
1 3 - P O I N T NUMBER, ETC. 2__Q--AREA NUMBER, ETC. Fig. 2.
--
j4
L. L ~ -
Mesh for Douglas D a m area.
....
/
366
J, E. AKIN
and the total area (of the triangular mesh) is
A=Z n=l
A n
so that the mean areal precipitation (of the triangular mesh) becomes R = Q/A.
By means of this procedure the mean areal precipitation is easily calculated if the gage locations and gage measurements are known. This procedure is easily programmed for computer application. As an example of the proposed technique consider the average precipitaTABLE 1
Node number
x
Coordinates* y
Gage measurement
Gage number
1 2 3 4 5 6 7 8 9 10 1l 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
35.0 43.4 47.3 48.0 45.9 43.7 46.8 49.0 40.4 24.4 15.7 8.2 4.6 0.0 13.8 24.4 32.5 34.7 37.1 29.8 31.0 43.1 34.5 19.0 18.5 25.8 20.6
0.0 3.2 9.7 15.7 22.5 32.8 38.6 48.6 48.4 53.1 45.7 38.7 30.3 20.7 18.7 6.2 4.3 11.3 21.3 28.6 35.8 42.4 45.3 37.9 28.0 22.0 15.0
9.99 15.69 13.56 9.51 18.36 14.52 13.10 10.35 8.83 4.76 4.17 3.81 4.19 4.62 5.45 4.96 19.52 11.4l 4.26 5.27 4.90 10.48 5.21 7.23 6.80 6.64 6.26
663 283 280 278 277 233 235 378 231 294 214 213 212 593 538 249 190 254 267 262 672 234 224 475 255 26l 656
* In miles from an arbitrary reference.
CALCULATION OF MEAN AREAL DEPTH OF PRECIPITATION
367
tion in the Douglas Dam (Tennessee, U.S.A.) drainage area for September 1957. Figure 2 shows the drainage area, the twenty-seven gages, and the thirty-five sub-areas. Using the data tabulated in Tables 1 and 2 one obtains an average areal precipitation o f / I = 7 . 7 l in. By way of comparison the arithmetic mean method gives/t--8.26 in. ; the Thiessen method (see Fig. 3) yields/7--7.91 in. ; and the Isohyetal method yields (see Fig. 4 ) / t - - 7 . 7 6 in. This example illustrates that the proposed method will yield results that are similar to the commonly used methods. TABLE 2 Sub-area
I
J
K
Sub-area
I
J
K
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
1 2 3 3 4 4 5 6 6 6 6 7 8 9 9 10 11 11
2 18 18 4 19 5 6 20 21 22 7 8 9 23 10 11 24 12
17 17 2 18 18 19 19 19 20 21 22 22 22 22 23 23 23 24
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
12 13 13 13 15 15 15 16 16 18 18 19 26 20 25 21 21
13 25 14 15 26 27 16 17 18 26 19 20 20 21 21 23 22
24 24 15 25 25 26 27 18 27 27 26 26 25 25 24 24 23
An advantage of the proposed method is that it allows the user to assign more "weight" to those gages that are most accurate. This is a result of the freedom of choice of the triangular mesh connecting the gage points. For example, in Fig. 2 point number 18 (gage no. 254) is used in the calculations associated with the eight surrounding sub-areas (numbers 2-5 and 26-29), whereas point 17 (gage no. 190) is used in the calculations of only three sub-areas (numbers 1, 2, and 26). However, if the measurements at point 17 are known to be more accurate than those at point 18 then that portion of the mesh can be rearranged as shown in Fig. 5. Then point 18 is associated with four areas (4, 5, 28, 29) and point 17 is associated with seven areas ( 1 4 , 26-28). This change in the mesh results in some minor changes in the input data and yields a new average precipitation of R = 8.01 in.
368
J.E. AKIN
In general the average precipitation calculated will depend on the choice of the mesh. Of course, in the limit as the sub-areas become smaller the linear interpolation would approach the true precipitation distribution. However, if one must choose a mesh for a coarse distribution of gage points, of equal accuracy, it is best to try to have the sub-areas of approximately
\
®
1278
656
254
/ ®255- LOCATION AND NUMBER OF RAINFALL STATION
Fig. 3.
~_~ (~ 249
------./
9 280
J ,,._..~0
Thiessen network.
equal size, and to have approximately the same number of lines radiating from each gage point. Another advantage is that the method is easily automated, and if a gage is added or removed from service this requires only some minor changes in the data preparation. In addition to the above features the output of the proposed method can be used to (automatically) generate a M a x i m u m Average Depth - Area Plot. This is accomplished by calculating the area, A,, and average precipitation,
CALCULATION OF MEAN AREAL DEPTH OF PRECIPITATION
,.'
..
4
i/
/ P'~ ~
/
,,'~ )
(/I/(t,,"
'
L-
\ ",L®
L______J~ / \
CONTOUR INTERVAL= 2"
\\
>~
®
--~/~----~/m b
<-----A .6 I
<~i
-~,~..~ ,
m-'~:o;_~, Fig. 4.
Isohyetal map.
2_.99 14
27 3 16
Fig. 5.
I Change in part of the original mesh.
369
370
J . E . AKIN
17, = Q . / A . , of each triangular sub-area, and ranking (from highest to lowest) the sub-areas in terms of their average precipitation. The points for the plot are then calculated by combining the sub-areas, in the order of their rank, and calculating the average precipitation of each combined area formed in this manner. Figure 6 shows a comparison between such a plot obtained from the Isohyet analysis and two plots obtained from the proposed method (using the two meshes discussed earlier). The plots compare well but the proposed method tends to underestimate the average maximum depth of
19
~ I
I
16.
a_ 14LIJ
-ISOHYETAL METHOD
Z 13
o_ F-
o_E
MESH NO. Z
QELI 9 <~ IE LIJ
MESH NO. I ~
~ ~ ~ . ~ . ~ . ~
8
6 5
I
2
~)0
[
4~
i
,
600
I
I
800
I
I
I000
I
AREA, SQUARE MILES
Fig. 6. Depth-area-duration curve for Sept., 1957. rainfall at the lower end of the area scale. This will be true unless the triangular areas are relatively small in the regions of maximum rainfall. An apparent disadvantage of the proposed method is that in general the perimeter o f the watershed will not coincide with the perimeter of the triangular mesh. This difference in enclosed areas will not be great if, as in previous example, a reasonable number of gage points are located on or
CALCULATION OF MEAN AREAL DEPTH OF PRECIPITATION
371
near the watershed boundary. However, in the Isohyetal and Thiessen methods stations outside the basin are also used to determine the average precipitation depth over the total area of the basin.
Conclusions The proposed method has been shown to yield results that are comparable with the commonly used techniques. Its major advantages are its simplicity and the fact that the procedure is easily automated. However, there are certain limitations associated with the proposed procedure that should also be mentioned. This procedure assumes a linear distribution of precipitation between stations. Thus its accuracy is directly dependent on the accuracy of this assumption. There are many cases where the true distribution of precipitation is not linear. In these cases the accuracy of the calculations would depend on the fineness of the mesh. A coarse linear distribution probably would not adequately model the true distribution. In such cases, one should obtain better accuracy by using the nonlinear distribution model presented in the Appendix. Of course, if the heaviest precipitation occurs at locations without precipitation gages then the average precipitation could be greatly under estimated. Finally, no specific criteria has yet been developed to guide the user in the optimum choice of the mesh used to cover the watershed (and portions of its surroundings).
Acknowledgement The author would like to thank Professor Bert Mullins who performed the calculations given for the mean station method, the Thiessen method, and the Isohyetal method.
References 1) S. S. Butler, Engineering Hydrology (Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1957) 2) O. C. Zienkiewicz and Y. H. Cheung, The Finite Element Method in Structural and Continuum Mechanics (McGraw-Hill Ltd., 1967) 3) A. S. C. E., Hydrology Handbook, Manual of Enginecring Practice, No. 28, (1949)
Appendix 1 QUADRILATERAL SUB-AREAS
The mesh used in the proposed method of analysis could be extended
372
J.E. AKIN
to include qttadrilateral sub-areas like those shown in Fig. 7. This would reduce the amount of input data and should increase the accuracy of the calculations, but it would complicate the computer program. As shown in Fig. 7, each sub-area would have four corner points and thus four gage measurements associated with it. The precipitation at any interior point could be assumed to be (12)
H , ( x , y ) = nO~1 -[-.O~2X -~ nOt3y + nO~4xy.
(Q)
//•Hk
Hj
'L N
I (xl,Yl)
(b)
Fig. 7. Quadrilateral mesh.
This relation is similar to Eq. (1), but it allows a nonlinear variation of precipitation overthe sub-area. The e q u i v a l e n t f o r m o f E q . ( 2 ) i s : H i = n~l + n~zXi + n~3Y~ + n~4Xlyi H i = n~l + n~2Xj + n~3yj + n~4Xjyj H g = n ~ 1 + n~2Xk + n~3yk + n~4Xkyk
Ht = n~l + n~2Xt + ,~3Yt + ,,~4xkYt"
(13)
CALCULATION OF MEAN AREAL DEPTH OF PRECIPITATION
373
The four ~'s can be determined by solving these four equations; a step which is perhaps most easily carried out numerically. The remaining operations would be the same as those in the main text of the paper except that the integration of Eq. (7) would probably be carried out numerically.
C C
MEAN
AREAL
PRECIPITATION
CALCULATIONS
C L
DIMENSION X I [ O O l l Y(LOO), H I I O O I , NGAGE(LOO|, I I ( 2 0 0 l , JJ(2OO), KK(2OO)p AREA(2OO)t HSUB(2OOIt QI2OO)t A I 2 0 0 ) , B(2OO)y C I 2 0 0 ) , TITLE(20) DIMENSION IRANK(20O) READ ( 5 , L O ) TITLE IO FORMAT ( 20A4 ) WRITE ( 6 j t O ) T I T L E READ ( 5 , 1 } NOPTS~ NOAREA I FORMAT ( 215 l DO 2 K = i , NOPTS, 2 READ ( 5 , 3 ) L , X ( L ) ~ Y ( L ) , H ( L ) ~ NGAGE(L) 3 FORMAT ( 1 5 , 3 F 1 0 o 0 , 15 ) DO 4 K = £y NGAREA, 4 READ I 5 , 5 ) L, II(L), JJ(L), KK(LI 5 FORMAT ( ~ I 5 ) QTOTAL = 0 . 0 ATOTAL = 0 . 0 DO 6 L = I , NOAREA, I I = TILL} J = JJILI K : KK(L} A(I) = X(J)*Y(K) - X(K)*Y(J) A(J) = X(K)*Y(I) - XII)*YIK) A(K) = X ( 1 ) * Y ( J ) - X ( J ) * Y I I ) B ( 1 ) = Y ( J ) - YIK B l J ) : Y(K) - Y ( I BIKI : Y(1) - YIJ C(1) = X(K) - X(J C ( J ) = X ( 1 ) - X(K C(K) = X ( J I - X I I AREA(L) = I A ( 1 ) ÷ A I J ) ÷ A I K ) ) I 2 . IF ( A R E A ( L I . L T . O . O ) WRITE (6,II) L I I FORMAT ( 5X,IBH*ERROR* TRIANGLE , [ 5 , 2 0 H HAS A NEGATIVE AREA I Q(L| : ( A ( I I * H I I I ÷ AIJ)*H(JI * A(K)*HIK) )/ 2. 2 * ( B(1)*H(1) + B(J)*H(J) + B(K)*H(K) )*I XIII ÷ XIJ) ÷ XiKI)I6. 3 * ( C(1)*H(1) ÷ CIJ)*HIJ) • CIKI*H(KI )*( Y(1) ÷ Y(J) ÷ YIK))/6. HSUB(LI = Q I L ) / A R E A ( L ) ATOTAL = ATOTAL + AREAIL) 6 QTOTAL = ~TOTAL + G I L l HBAR : QTOTAL / ATOTAL WRITE ( 6 ~ 7 ) N O P T S , NOAREA, HBAR, ATOTAL 7 FORMAT ( I H ,15HNO. OF POINTS = , 1 5 , 3Xt L4HNO. OF AREAS = , 15, 2 3X, 23HAVERAGE PERCIPITATION = , E [ 5 . 8 , 3X, 12HTCIAL AREA = • 3 El5.8 lid WRITEib,8) ( Le X ( L ) e Y I L ) ~ H i l l NGAGE(L), L=I,NOPTS) 8 F O R M A T ( S × , 5 H P O I N T , I g X , [ H X ~ I g X ~ I H Y , I q X , I H H ~ I 3 X , T H G A G ENO / ( I [ O , 2 3E20.8~I20 ) I WRITE 1 6 ~ 9 ) (L,IIiLI,JJ(LI,KK(L),AREA(LI~HSUB(LI,L=I~NOAREA) ~ FORMATIIHI~IOH AREA N O . ~ g X , I H I , g X , I H J , g X , I H K , ~ 5 X , 4 H A R E A ~ I S X , 2 4HHSUB / ( I X , 4 I I O , 2 E I 9 . 8 ) ) CALCULATXCNS FOR MAX DEPTH-AREA-DURATID~ PLOTS JMAX = N~AREA - 1 HMAX = H S U B ( I } 2 3
2 3 4 5 6 7 8 9 lO l[ 12 13 14 15 £6 17 18 Ig
20 2I 22 23 24 25 26 27 28 29 30 3[ 32
33 34 35 36 37 38
39 ~0 4l 42 C 43 4~
Fig. 8a
374
J. E AKIN kS 46 4l 48
[2
49 50
51 52 53 54 55
13
O0
56 57 58 5q 60 61 62 63 b4 65 66 67 68 69
7C 7L
72 73
76 75 76
HMIN = HMAX DO 12 I = I, NOAREA IF ( HSUB(1).GT.HMAX ) HMAX = H S U 8 ( I ] IF ( H S U B ( I ) . L T . H M I N ) HMIN = H S U B ( I } T I = HMAX T2 : HMIN DO 13 I : L, NOAREA IF ( H S U B ( I ) . E Q . T I ) [RANK(I] = [ IF l H S U S I I ] . E Q . T 2 ) IRANKINCARFAI = I O0 1 4 J = 2 , JMAX 15 I = l, NOAREA ( HSU~(I).GT.T2.AND.HSUB(I).LT.FI ] GO TO 15 T2 : HSUBII} IRANKIJ] = [ CONTINUE fI = T 2 T2 : HMIN ATOT = 0 , 0 QTOT = O . O DO 1 7 J = I, NOAREA [ = [RANKIJ] A T O T = ATGT + A R E A ( [ ) QTOT = WTOT + O ( I ) A(J) = ATOT
IF
I6 lb 14
GO TO 16
17 B ( J ) = QTOT/ATOT WRITE ( 6 , I 8 ) 18 FORMAT ( I H I , I I X , Z 6 H P O I N T S FOR AREA-DEPFH PLOT /I ] WRITE ( 6 , 1 9 | ( J, B(J),A(JI,IRA~K(J),J= I,NOAREA ) L9 FORMAT (~X,4HRANK, IOX,HMAX OEPTH,£X,IOHTOTAL AREA, 2 5X,IOHAREA ~OOEO I ( [ t O , Z E I . 8 , [15 ] ) STOP END
Fig. 8b Fig. 8.
PCINI | Z
5 5 8 g lO It [2 1{ 16 15 lh 17 18 19 20 2J 22 23 2W 2b 26 2t
O.35000000E 0.43399qOOE 0o67300000E O.~EOCOOOOE 0.65809990E O,k3609990E C,468GOOOOE 0o40000000F 0.603qO990E 0°24399990E O,[56qqO90E O,8tqq 9906 0°660C0000£ O°OCCCO000E 0.[3800000E 0.2439~990E O.~25GUOCOE I).346q ggOF O.37100000F 0.29800000~ 0.310000COE O.k3|OOOOOE 0,365000COE 0.[9000000£ 0.[85C00C0£ 0.25800000£ (%206000006
X 02 C2 02 02 02 02 02 02 02 02 02 Ol O! ~0 (12 02 02 02 C2 OZ 02 02 02 02 02 02 02
Mean areal precipitation calculations.
0.00000000E O,~1999q OE 0o96999990E O,[56qqOOOE 0.22500000E 0,32800000E 0,~8600000E O.4BbCOCOOE O,63BOOOOOE Oo53IO000CE 0.~569g990E 0°38600990£ 0.30300000E O.206900qOE O.186999qOE O,6[£gg£g0E 0.45000000E O,|t3OOOOOF 0°21300C00E O°Z8600000F 0,~5800000E 0°62399990~ 0,65300C00E 0,5789g990E O.ZSOOOOOOF 0.22000000F 0,]50000006
Y 00 Ol O[ 02 02 02 02 02 02 02 02 02 02 02 02 C[ O[ 02 02 02 02 02 C2 OZ OE 02 02
Fig. 9a
O.OO899990E 0.[5689990F 0.~560000E O.05[O0000E 0.18560000F 0.14520000E O.i3IOOOOOE 0.103500C0~ 0.88299990E 0.47600000E 0.617000006 0o38[00000F 0.4[899990E 0.46[9qOOOF 0.56499990E 0.4960OOOOE 0,t95200006 O.I[~OOqOOE 0.42600000E 0.~27000(;0F O,489990qOE 0.104199906 ~°b2[00000£ 0.72299990E 0,68000000F 0°56600000£ 0.~2600000E
H OI OZ 02 O[ 02 02 02 02 O[ 0[ 0| Ol O[ O[ O| O[ 02 02 O[ O[ O[ 02 0[ O[ 01 O! O[
GAGE NO 663 28~ 280 278 277 233 255 378 2It 29~ 2[6 2[~ 2[2 50~ 5~8 ?6q t90 254 267 262 672 gJ6 224 415 255 261 656
~ o ~ o o ~ o o
o. .o .o .o .o .o .o
~
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
~
.
o
.
.
.
.
.
.
.
.
.
.
.
-
o
°
°
.
~
m
~ .....
o ~ o o o o o o o o o o o o o o o o o o o o o o o o o o
o ~ o o o : , o o o o o o o o o o o ~ o o o o o o o o o o o o o o o o o
~
. . . . . . . . . . .
~
og~ .o .o .o .o .o ~ g .o .o .o .o .o ~ .o o.o . . o.o o. o. ~ o. o. o. o. o. o. o. mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
~'~ooo
o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o . . . . . . . . . . 0 , . ° , . . ° . . . . . ° , 0 , . . 0 . , . , .
.
~
0
.
.
o~oo
mmmmm
ooooo
o o o o o o o o o o o o o o o o o o o ~ o o o o o o o o o o o o o o
m m
mmmm~%mmmmmmmmmmmmmmmmmmmmmmmmmmm~mm~ o o o o o o o ooooo
~"
mmm
~ ooo~og
~
o .....
o o o o ~ o o o o o o o o o o o o o o o
mmmmmmmmmmmmm~
oog
.....
~ ~
~
~
'-,
O ~ ~ '
N
0
~
z
~ o o o o o o o o o o ~ o o o o o o o o o o o o o o o o o o o o o o 0° .... ,,,0 . . . . . . . . . . ..... , . 0 , . , ° ° , ,
N
376
J.E. AKIN
Appendix 2 The computer program shown in Fig. 8 uses the data shown in Tables 1 and 2 and calculates the following quantities: the area and average precipitation for each triangle (Fig. 9b), and the area and average precipitation for the entire triangular mesh. The above quantities are then used to calculate the points for a maximum average depth versus area plot (Fig. 9c). The input format should be changed if the program is applied to a problem where a single mesh geometry is used and various sets of gage measurements are run.