A Flood-Frequency 

Relation Based on 

Regional Record Maxima

GEOLOGICAL SURVEY PROFESSIONAL PAPER 434-F

Prepared in cooperation with the 

U.S. Atomic Energy Commission
A Flood-Frequency 

Relation Based on 

Regional Record Maxima

By P. H. CARRIGAN, JR.

STATISTICAL STUDIES IN HYDROLOGY

GEOLOGICAL SURVEY PROFESSIONAL PAPER 434-F

Prepared in cooperation with the 

U.S. Atomic Energy Commission

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1971
UNITED STATES DEPARTMENT OF THE INTERIOR 

ROGERS C. B. MORTON, Secretary

GEOLOGICAL SURVEY 

W. A. Radlinski, Acting Director

Library of Congress catalog-card No. 75-611783

For sale by the Superintendent of Documents, U.S. Government Printing Office 

Washington, D.C. 20402 - Price 40 cents (paper cover)
CONTENTS

Abstract._ _________________________________________

Introduction _______________________________________

Acknowledgments _-_--__----_______________________

Conover and Benson method of determining probabilities.

Method for determining probabilities from dependent

records ______-__---_-_--___--___-____________

Simulation model.______________________________

Flow chart for simulation,_______________________

Influence of interrecord correlation on exceedance probability. ___-__----_-------_---__--_____-_____

Equicorrelation matrix __________________________

Verification that p, is distribution free_________

Page 

Fl

1

2

2

3

3

3

3

4

4

Page

Influence of interrecord correlation on exceedance probability C ontinued 

Equicorrelation matrix Continued

Variation of p* with k and n_________________ F5

General correlation matrix.______________________ 5

Application of the method in the Big Lost River basin,

Idaho _____________-____---____-_-___-_--_-_--_._ 7

Discussion _________________________________________ 7

Conclusions._______________-_-__--__-__-_-_____-___ 9

References.________----____--______-___---_-__--___ 9

Computer program___--_-___---__--_-_-__--_-__--___ 11

ILLUSTRATIONS

FIGURE 1. Flow chart for program to simulate {pi(k); i=l, 2, ... n}_______________-____---___--_-___---_--------

2-8. Graphs:

2. Variation of {?,-(&); i= 1, 2, 3} with p for k= 10._______________-_______----_____--_-_---__-----

3. Variation of Pi(k) and p2 (k) with p for bivariate gamma and normal distributions of annual extremes; 

fc=10_.. ---__-______________________-_________-________-___--_______-_---_-------------_

4. Variation of pi with k for n 3 and p = 0, 0.6, and 1_____________________________________________

5. Variation of Pi(k) with n for k = 10 and p=0, 0.6, and l_________________--___----__---_____---_-

6. Effect of mean square deviation _ 32 on [pi'(k); i=l, 2, 3} for k= 10 and p = 0.5__-______---_________

7. Upper limit of p/ for -fm«__--_-__-______--____--______-______--__-_-----__---_-_----___---_-

8. Flood-frequency relation for Big Lost River below Mackay Reservoir_____________________________

Page 

F4

TABLES

TABLE 1. Name of station and period of record for floods in central Idaho having identical distributions.______________

2. List of annual flood records for analysis of flood-event frequency in the Big Lost River basin, Idaho, for the 

concurrent period 1922-67-_-_________________________________-_-____-_-_---_-_------------.------

3. Plotting positions for the relation between the rare-flood magnitude and rare-flood frequency at Big Lost River 

below Mackay Reservoir.______________________________________________________---__-------_-----

Page 

F7

STATISTICAL STUDIES IN HYDROLOGY

A FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD

MAXIMA

By P. H. CAEEIGAN, JR.

ABSTRACT

Flood intensity varies randomly in both time and space over a 

region. In flood-frequency analysis the regional variation is commonly related, through regression analysis, to variations in basin 

and climatic characteristics. Thus, the sample size is reduced to 

one time sample applicable to the entire region. This type of 

analysis virtually limits the recurrence-interval estimate for the 

maximum annual flood to the length of time for which records are 

available.

Techniques of analysis introduced in this paper increase the 

sample size by taking account of the random variations of flood 

intensity in both time and space. The magnitude of increase in 

the sample size is limited by the statistical dependence between 

annual flood records for the region; that is, the less the correlation 

between records the greater the effective sample size. As a result 

of increasing the sample size, the recurrence interval for the maximum event is increased; the increase may conceivably be as great 

as the number of station-years of record, if the records are 

independent.

The recurrence interval for the maximum annual flood in the 

region is estimated by a computer-simulation model. In order to 

apply the model, the assumption is made that the sets of concurrent records for a region are reducible to identical distributions bv 

use of an appropriate scaling factor.

The dependence among records is preserved by including in 

the model the sample correlation coefficients for all pairs of 

record sets. By means of a principal-component transform, the 

matrix of correlation coefficients is used in connection with a 

normal random-number generator to generate synthetic sets of 

normally distributed flood intensities equal in number and duration to observed sets of records. The maximum event from each 

synthetic record is found, and the maxima are ranked in descending order. The exceedance probabilities (reciprocals of recurrence 

intervals) of the ordered maxima are repeatedly determined to 

estimate the average normal exceedance probability for each 

order. The scheme of computing exceedance probabilities of known 

distribution for ordered maxima selected from synthetic records 

of known distribution provides an estimate of probability that 

is independent of the statistical distribution of the observed 

floods.

The construction of a flood-frequency relation for any station 

in the region can be accomplished by graphically relating the 

observed ordered maxima, reduced to a common base, to the

corresponding estimated recurrence interval. The procedure of 

constructing a flood-frequency relation is illustrated by applying 

it to floods in the Big Lost River basin near the National Reactor 

Testing Station in southeastern Idaho.

INTRODUCTION

Current practices of regional flood-frequency analysis 

define regions in which the records of annual hydrologic 

extremes have, for the most part, identical distributions. 

The method of analysis is that of multiple regression, 

whereby floods of specific recurrence intervals or the 

statistical parameters of flood distributions are related 

to basin and climatic characteristics. The method 

composites the experience over a region but, except by 

extension based on the assumed distribution of events, 

it produces results that are limited virtually to the 

length of the period of sampling. Thus, no advantage is 

taken of the facts that spatial sampling to some degree 

is equivalent to time sampling and that some of the 

experience within the region may furnish information 

about the expectancy of events for a period much 

longer than the actual period of record. Work by 

Conover and Benson (1963) suggests that this information can be realized by appropriate treatment of 

combined records from several streamflow stations in 

the defined area. In the Conover and Benson method of 

analysis, the process is assumed to be ergodic, so that 

spatial sampling is equivalent to time sampling of 

extreme events. The result is to increase the size of the 

sample to nk events in a region (fc=years of record, 

n number of sets of record), provided the records are 

independently and identically distributed.

If the Conover and Benson method of regional floodfrequency analysis were applied to records of annual 

extremes, the records could be reduced to identical 

distributions by normalizing them to the magnitudes

Fl
F2 STATISTICAL STUDIES IN HYDROLOGY

of the mean annual events at each streamflow station. 

However, the records for a region usually are not independent but are correlated. Hence, another method of 

analysis must be developed for treating the combined 

information from several identically distributed sets of 

records. This report presents a procedure to estimate 

exceedance probabilities of floods from records that are 

dependent and identically distributed. The maximum 

recurrence interval for the region, estimated by this 

procedure, lies between fc-fl and/(p)fc-fl where /(p) 

is a function of the correlation between record sets 

(l</(p)<r&) for the range in correlation coefficients

The scheme of Conover and Benson for treating the 

combined information from independent and identically 

distributed records is to compute expected exceedance 

probabilities (independent of the distribution) for the 

ordered record maxima. This scheme is the key to 

developing regional frequency relations for dependent 

records.

Once a region having identically distributed sets of 

concurrent annual flood records is defined, the regional 

flood-frequency relation is developed as follows:

1. The records are reduced to a common base.

2. The maximum events from each reduced record are 

ordered.

3. The correlation matrix for the records is computed.

4. Computer routines developed in this study are used 

to provide, through simulation alogorithms that 

utilize the correlation matrix, estimates of exceedance probabilities associated with ordered maxima 

from n, number of fc-year records.

5. The frequency relation is defined graphically by 

relating exceedance probabilities to the ordered 

maxima. Reversing the transformation process that 

reduces records to an identical distribution produces regionalized frequency relations for annual 

floods at individual stations.

After briefly illustrating the Conover and Benson 

method of obtaining exceedance probabilities for independent flood records, this report introduces the 

mathematical and computer-simulation foundations 

for estimating exceedance probabilities for dependent 

flood records. Then follow discussions on the application 

of the procedure to simple models which help to describe 

the influence of the number of records, correlation coefficients, and record length on exceedance probabilities; 

the development of regional flood-frequency relations 

in the Big Lost River basin, Idaho; and the extension 

of the procedure to incomplete and to serially correlated 

records. The complete computer program is presented 

at the end of the report.

ACKNOWLEDGMENTS

This study of flood-frequency relations was part of 

a research project on the statistical occurrence of extreme hydrologic events in relation to the design and 

selection of sites for nuclear reactor facilities. The work 

was done on behalf of the Division of Reactor Development and Technology, U.S. Atomic Energy Commission.

Early considerations for estimating the probability 

of occurrence of rare events favored analytical approaches. However, after discussions with W. J. 

Conover of Kansas State University, it was decided 

to discard methods relying on analytical approaches. 

U.S. Geological Survey personnel made the following 

contributions to the study: M. A. Benson recommended 

use of Monte Carlo techniques; discussions with N. C. 

Matalas resulted in much of the structuring in the 

Monte Carlo model; E. J. Gilroy drew attention to the 

use of a moving-average scheme to simulate multivariate distributions; and M. S. Hellman kindly supplied an improved computer algorithm for finding the 

eigenvectors of a correlation matrix.

CONOVER AND BENSON METHOD OF 

DETERMINING PROBABILITIES

Assume that n concurrent sets of k annual extremes 

are obtained from n identical and independent distributions. The maximum events, one from each record, 

are ordered from highest to lowest. The probability 

that another event j exceeds these ordered maxima

{/,; i=l, 2, . . .n} is

4=4 ______r&! 

-Z_/ ji

771=0 / \IT-r

for i=l, 2, . . .n (Conover and Benson, 1963, p. E159). 

Example. The following three sets of four random 

events {x«; i=l, 2, 3 and j=\t 2, 3, 4} are from 

independent and identical distributions:

fn=3, fc=4]

1 3 69

2 38 24

3 17 61

4 32 30

3

48

60

83
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA F3

The ordered record maxima {/*; i=l, 2, 3} are, respectively, 83, 69, and 38. The probability that another 

event / exceeds the second highest maxima is

3!

771=0

j=Q

3! 3!

314(3+1) 2,4(3+!) (3 + 1-

METHOD FOR DETERMINING PROBABILITIES 

FROM DEPENDENT RECORDS

The probability of another random event/exceeding 

one of the ordered record maxima jt where the records 

are from dependent distributions cannot be analytically 

determined. The probability can be determined, however, through use of Monte Carlo simulation techniques. 

In this application of simulation for dependent records, 

it is assumed that the probability of another random 

event exceeding an ordered record maximum is independent of the identical distributions for which the 

records of random events are a sample. Such an assumption has been verified by Conover and Benson (1963) 

and proven by Conover (1965) for independent random 

events.

The simulation algorithm, which incorporates a 

normal random-number component having zero mean 

and unit variance, generates n sets of k events that 

simulate the dependence among n records of k hydrologic extremes. The maximum event yt from each set 

of normally distributed variables {xn, xi2, . . . xik ; 

i=l, 2, . . . n} is found, and the maxima are ordered 

from highest to lowest, so that 2/i>2/2> . . .^>y n- All yt 

are from normal distributions having zero mean and 

unit variance, so their exceedance probabilities are 

given by

1 C' <° ~^2

P[y>yi]=j= J e~*~dy,toTi=l,2,...n.

By repeated generation of this simulated set of records, 

estimates of the average exceedance probabilities for 

the ordered maxima are obtained.

SIMULATION MODEL

The simulated sets of records are generated by the 

normal multivariate model

where X=matrix of n groups of k events,

=nXk matrix of independent normal random 

numbers with zero mean and unit variance, and

B=E\=riXn principal-component transform 

matrix for the records of annual extremes, 

in which

E=eigenvector matrix associated with the correlation matrix R obtained from n records 

of k hydrologic events, and 

A=diagonal matrix whose n diagonal elements 

are the square roots of the eigenvalues 

forE".

(See Morrison, 1967, p. 221-247.) This model is based 

on techniques of synthetic hydrology developed by 

Matalas (1967, p. 940).

After generation of X, the n maxima are selected and 

ordered and the exceedance probabilities are computed. 

The algorithm for computing normal exceedance 

probabilities is that of Zelen and Severs (1965, p. 932, 

eq. 26.2.19).

FLOW CHART FOR SIMULATION

The average exceedance probability

- .>yn ]

1,000

was estimated through use of a digital computer. The 

main program flow chart is shown in figure 1. Routines 

for computer input-output listings of matrices, matrix 

operations, and generation of random numbers are 

adaptations of subroutines developed by IBM (1968) 

for use in their 360-series computer systems. Formal 

declarations of most subroutines have been eliminated 

from the program to increase efficiency. A complete 

listing of the Fortran language program for the computations is included in the computer program at the 

end of this report. The parameters n and k may each be 

as large as 100 using the present configuration of the 

IBM 360/65 computing system (380,000 bytes available 

storage) operated by the U.S. Geological Survey.

INFLUENCE OF INTERRECORD CORRELATION ON 

EXCEEDANCE PROBABILITY

An n-record correlation matrix contains n(n 1)/2 

correlation coefficients. If r^>3, problems in management of matrix elements begin to obscure the study of 

the influence of interrecord correlation on exceedance 

probabilities. Hence, most simulations of exceedance 

probabilities were derived from the simple, yet informative, 3-record correlation matrices.
F4 STATISTICAL STUDIES IN HYDROLOGY

c Start Read n, 

k, and R

 > Print title 

n, k, and R

 > 

Call Eigen 

subroutine 

to compute 

E and eigenvalues

 > Compute 

A

-0

Compute cumulative 

sum of P[y>y. ],2p.

FIGURE 1. Flow chart for program to simulate (p,-(/c); i=l, 2,.. .n}.

EQUICORRELATION MATRIX

Using the equicorrelation matrix R (all correlation 

coefficients equal)

p P P~\ pip 

LP P IJ

for 0>p>l, the expected probabilities \pi\ i l, 2, 3} 

were simulated for Ar=10, 20, ... 50.

As shown in figure 2, pt for any order i changes 

monotonically, as had been anticipated very gradually 

from the theoretical probability at p=0 (three independent records) to about p=0.8, and then rapidly 

to the common theoretical probability at p l (three 

identical records). Comparison of the simulated probability to the theoretical probability at p=0 and at 

P=l, shown in figure 2, indicates that 1,000 iterations is 

sufficient for estimates of p t. Similar monotonic variations in {p^ i=l, . . . n} between p=0 and p=l occurred for n=2 and n=5 where Ar=10, 20, ... 50.

VERIFICATION THAT A IS DISTRIBUTION FREE

The probabilities p t for n=2, k=lO, and p=0, 0.1, 

0.2, ... 1.0 were simulated by the previously de0.18

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

FIGURE 2. Variation of {pi(/c); i=l, 2,3} with p for *=10.

scribed normal variate generating model. The data for 

these simulation runs are plotted with their related 

curves in figure 3.
0.10

0.08

0.06

0.04

0.02

0

FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA

0.10

F5

TYPE OF VARIATE 

 Norma I 

x Gamma

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

P

FIGURE 3. Variation of pi(k) and p^(k) with p for bivariate gamma and normal 

distributions of annual extremes; fc=10.

Another model for simulating a multivariate distribution \Zi_j] i=l, 2, . . . n and j=l, 2, ... k} is to 

let

in which u(l-i) d+i. .,= exponentially distributed independent random numbers having unit 

mean and unit variance, and 

2 i( ^=random numbers following a marginal gamma distribution having 

mean ra and variance ra, for 

which the correlation coefficient

s

P(Zi,

m d

Using this simulation model the mean exceedance 

probabilities {p { ; i=l, 2} for the ordered maxima from 

ztij (for ra=5, n=2, and &=10) were estimated for 

P=0.2, 0.4, 0.6, and 0.8 (d=4, 3, 2, 1); the coefficient 

of skew for zitj was 0.89 (the mean coefficient for a 

50-record sample of annual floods was 0.9).

Results for this bivariate gamma method of estimating pi are shown in figure 3. Within the experimental 

error experienced, the estimates appear to confirm the 

assumption that pt is distribution free.

VARIATION OF p t WITH k AND n

The probability pt varies smoothly and monotonically 

with either k or n. (See figs. 4 and 5 for examples.) The 

variation with k of p t for 0<p<l and n constant is 

similar to the variation of the theoretical probability 

(p=0 and 1) with k. The relation between p t and n for 

constant k shifts from being highly curvilenar for p 

near zero to being nearly horizontal for p near 1.

0.08

0.06

0.04

0.02

p=0.6

p=0

10 20 30 

k

40 50

FIGURE 4. Variation of pi with k for n=3 and p=0, 0.6, and 1. 

GENERAL CORRELATION MATRIX

The relative insensitivity of p^ to changes in p (fig. 2) 

suggests the use of a simple parameter

n-l

Pn =

n(n S P«.

This parameter weights the central tendency of correlation coefficients in the general correlation matrix 

(pij are not all equal) and replaces p in estimating 

Pi (the prime referring to a general matrix) from relations such as those in figure 2. Results of tests showed 

that all the variation in p tf could not be accounted 

for by variations in ~p n . A parameter that described the 

dispersion of the coefficients in R about p», such as a 

mean deviation

n-l

"" 71(71-1) tifr

or a mean square deviation

o n 1 n «52 ^ ^~

UM 

\PiJ~Pn\

seemed necessary to explain variations in pt' more 

fully.

Experiments using third-order correlation matrices 

were performed to define the variation of \pi\ ^ = 1, 2, 

3 j with increases in a3 or 5Jj for constant p3 (see fig. 6). 

Under these conditions, ptr was found to increase
F6 STATISTICAL STUDIES IN HYDROLOGY

monotonically for i=l and i=2 and to decrease 

monotonically for i=3.

Under the constraints that R be positive, semidefinite 

(see Morrison, 1967, p. 60), ~p n be constant, and 0<py 

<1, a maximum a n exists. For p3 >%, a3 is maximized 

if one of the p«=l, say pi 2 =l, and the other coefficients

o.io

0.08

0.06

0.04

0.02

FIGURE 5. Variation of pi(k) with n for fc=10 and p=0, 0.6, and 1.

0.12

0.08

0.04

;=2

X Defined in figure 3 

for i = 1,2

0.04 0.08 

5?

0.12 0.16

FIGURE 6. Effect of mean square deviation Stf on !p,'(fc); i=l, 2, 3) for fc=10 and 

p=0.5.

are as follows:

_ 3p3 Pis P23 

hence, Pi3=p23=p for n 2. For ps<O£, a3 is maximized 

if pi2=3p3 and Pis=p23=0. The maximum is

These conditions for which a3 is a maximum obtain 

the maximum for d%. The upper limit for p/ at a 

maximum dl is shown in figure 7.

The variation of ptf with 532 , shown in figure 6, 

brings out a point about the maxima. For i=l, p\ 

increases to the limit pl for n 2; for i 3, p%' decreases 

to the limit p2 for n=2; but p2 f does not reach either 

limit. If Pi2=l, then only two ordered maxima occur, 

one in records 1 and 2 and one in record 3. One of the 

maxima may be first order and the other may be third 

order; the second-order maximum does not really 

exist. It would appear in figure 6 that p2 ' is approaching 

the limit for p=l (^=0.0909 for fc=10). Results 

from two other sets of tests for A:=50 and n=3 and 

for fc=10 and n=5 indicate that the variations of 

pf with Pn and 5 32 are much like those shown in figures 

6 and 7.

0.10

0.09

0.08

or

o

a-

 Computed limit

X Limit based on relation 

P! =f(p)forn =

0.07

0.06

0.05

0.04

0.03

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

FIGURE 7. Upper limit of pi' for fi|maz. Arrow indicates shift in plotting position.
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA F7

An attempt to relate ptf fully to a practical range in 

pn , d nz, and k appears to be economically prohibitive 

because of computer cost. Rather, these results suggest 

that p^ should be estimated for each specific correlation 

matrix and period of record.

APPLICATION OF THE METHOD IN THE BIG 

LOST RIVER BASIN, IDAHO

Frequencies of floods on the Big Lost River near 

the National Reactor Testing Station in southeastern 

Idaho are required in connection with an investigation 

for the Atomic Energy Commission of the potential 

hazard from extreme floods. The first step of the 

flood-frequency analysis is to correlate annual flood 

records in this part of Idaho with distributions that 

are identical save for some scaling factor.

From a report on regional flood-frequency relations 

(Thomas and others, 1963), 48 candidate records were 

selected for analysis. The coefficient of variation, square 

of the coefficient of skew, kurtosis, and the standard 

error of these moments (see Kendall and Stuart, 

1963, p. 228-233) were computed for each candidate 

record. The records of annual floods were assumed to 

have distributions identical with the selected base 

record (Big Lost River below Mackay Reservoir) 

if they satisfied the following criteria:

1. The confidence limits of each moment for the record, 

defined by its standard error, were within the 

limits of the same moment for the base record.

2. The median test (Ostle, 1963, p. 473) indicated the 

acceptability of the hypothesis of identical 

distribution with the base record at the 5-percent 

level of significance.

These two criteria were met by 15 records (see table 1). 

The next step in the analysis is to reduce the selected 

records to identical distributions. This is accomplished 

by normalizing them to the computed mean annual 

flood.

TABLE 1. Name of station and period of record for floods in 

central Idaho having identical distributions

Name Period of record

Beaver Creek at Spencer________________________ 1950-63

Lime Creek near Bennett__________________________ 1946-56

Challis Creek near Challis_________________________ 1944-63

Combined discharge, Big Wood River and Slough at

Hailey___________________________ ______________ 1936-57

Big Lost River below Mackay Reservoir near Mackay_ 1919-67

Surface inflow to Mackay Reservoir near Mackay____ 1919-67

Big Lost River at Ho well Ranch near Chilly.. _______ 1920-67

Salmon River at Salmon_________________________ 1920-67

Valley Creek at Stanley....________________________ 1921-67

Salmon River below Yankee Fork near Clayton______ 1922-67

Salmon River near Challis______-____._____________ 1929-67

Big Lost River at Wild Horse near Chilly_________ 1944-67

South Fork Boise River at Anderson Ranch Dam__ _ 1946-67

Big Wood River near Ketchum___________________ 1948-67

Birch Creek near Reno____________________________ 1950-63

At this point in the analysis, the objective is to select 

those of the 15 concurrent records whose combination 

maximizes nk to realize the dual advantage of obtaining 

large recurrence intervals and relating outstanding 

events in the region to these intervals. The records for 

which the station years are maximized are listed in 

table 2.

The correlation matrix R for the six station records

123 456 

1 0. 379 0. 418 0. 239 0. 256 0.415' 

1 . 941 . 711 . 728 . 802 

1 . 677 . 711 . 823 

1 . 971 . 914 

1 . 934 

1

and a k value of 46 were inputs to the computer program for which the maxima [p'i\ i=l, ... 6} were 

computed.

The ordered reduced variate maxima, p\, and the 

associated flood magnitude in the Big Lost River 

below Mackay Reservoir are listed in table 3. The 

flood magnitude is equal to the product of the ordered 

reduced variate maximum and the mean annual 

flood at this station. The recurrence interval, T, is

The rare-flood reduced variates are plotted against 

their recurrence intervals in figure 8. The annual 

floods, in reduced variate form, are also plotted in this 

figure to define the relation between flood magnitude 

and flood frequency for the Big Lost River below 

Mackay Reservoir.

DISCUSSION

The procedure for determining the probability of 

occurrence of rare flood events may be extended to 

other kinds of dependent records, only a few alterations being required in the computer program. If the 

records follow a first-order Markov process as well as 

being cross correlated, then the simulation model is

X&i =AXj + Btj+i ,

where Xj=riX 1 matrix of the jih events (.7=1, 2, ... 

k 1) inn records,

~

BBT=R0 -R1R0~1R1 T,

B=E\, where E is the eigenvector matrix

and X is the diagonal matrix whose

elements are the square roots of eigenvalues,

.R0 =lag-zero correlation matrix, 

JRj=lag-one correlation matrix, and 

 j=nXl matrix of random normal numbers

with zero mean and unit variance. 

(See Matalas, 1967, for more detail.)
F8 STATISTICAL STUDIES IN HYDROLOGY

TABLE 2. List of annual flood records for analysis of flood-event frequency in the Big Lost River basin, Idaho,

for the concurrent period 1922-67

Station

1

2 

3 

4 

5 

6

Coefficient of 

variation

Value Standard 

error

Big Lost River below Mackey Reservoir near

Big Lost River at Howell Ranch near Chilly--.. 

Surface inflow to Mackay Reservoir near Mackay .

Salmon River below Yankee Fork near Clayton.

0.417 

.394 

.491 

.474 

.354 

.380

0.033 

.035 

.052 

.048 

.036 

.035

Square of skew 

coefficient

Kurtosis Maximum flood

Value Standard Value Standard Cubic feet Reduced 

error error per second variate

0.397 

.170 

.0493 

1.16 

.217 

.153

0.289 

.165 

.089 

.71 

-.317 

.144

2.57 

2.49 

2.36 

3.46 

2.80 

2.64

0.49 

.43 

2.33 

1.05 

.63 

.34

2,640 

4,420 

3,240 

2,000 

10,300 

16, 500

1. 

2. 

2. 

2. 

1. 

1.

40 

13 

23 

07 

99 

90

TABLE 3. Plotting positions for the relation between the rareflood magnitude and rare-flood frequency at Big Lost River below 

Mackay Reservoir th

£k"V

Exceedance Recurrence interval Reduced variate Flood J? 

Order probability (years) maximum magnitude ol 

(cfs) ei

1 0. 00662 151. 1 2. 23 4, 190 ^ 

2 .0115 86.9 2.13 4,000 ai 

3 .0159 62. 7 2. 07 3, 890 ,/ 

4 . 0218 45. 8 1. 99 3, 750 tl] 

5 . 0298 33. 5 1. 90 3, 580 

6 . 0432 23. 2 1. 40 2, 640 j

10,000 

8000 

6000

4000 

Q z. 

o 

o

LJ 

</)

(X 2000

LJ 

CL

1- 

UJ 

LJ 

LJ_

2 1000

en

0 800

z

- 600

LJ

0 

g 400 

Q

200 

inn

1 1

x^JC-X^ ^e****^ 

^*^

s$^X

1 I

1 1 1 1

X X.

"X

1 1 1 1

^"^

Through special investigations for historical floods, 

e period of record for the maximum event may be 

tended from k to k+l years, where I is the period 

extension. The estimation of p tf for the (k+l)- 

rent record may be undertaken provided the correla3n matrix denned by fc-event records is considered 

)plicable to the longer period and the n maxima for 

e (k+l) -event records are denned. 

A somewhat analogous condition to the problem 

analyzing historical floods occurs when some records

i i

 

*

i i

i i i i

«^^^^*~

* r=kfi' r

.r = -V f 

p;

I I I I

_ ^-^^"<

i = rank o1 

jr ordered

^- 

observ 

maxim

. 

3d flood 

a

10 20 50 

RECURRENCE INTERVAL r, IN YEARS

100 150 200 300

FIGURE 8. Flood-frequency relation for Big Lost River below Mackay Reservoir based on the probability of occurrence of rare events; w=6, fc=46, and

p6= 0.661.
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA F9

covering a period k have missing events. Assuming 

known record maxima, some uncertainty about the 

correlation matrix is introduced for records of this kind. 

If the matrix is positive semidefinite, a test of the sensitivity of PI to changes in individual correlation coefficients would seem appropriate. The variation in the 

lengths of record may be so large that the correlation 

matrix is inconsistant. Some suggestions for developing 

consistent matrices were given by Fiering (1968).

The objective in the analytical procedure for developing a regional flood-frequency relation is to maximize the number of station years (increasing the 

sample size) which tends to minimize p % '. Because 

Pi varies almost directly with l/nk, the objective can 

be attained by seeking a large number of short records 

or a few long records. The objection to the use of 

short records is that the sample may be nonrepresentative in comparison with norms established by much 

longer sequences of events. Unless statistical comparisons between long and short records indicate nonhomogeneity, no reason exists for rejecting short 

records. Sampling errors will be larger for short records, 

but this disadvantage may be outweighed by the substantially larger values of nk that will be provided by 

the shorter than the longer records.

CONCLUSIONS

The technique of regional flood-frequency analysis 

that has been introduced takes advantage of the random variation of flood intensities in both time and space. 

Previous methods of analysis are virtually limited to 

the consideration of variation in time. A result of applying this new technique is to increase the recurrence interval for the maximum flood event; because of correlation between records, it is unlikely that this recurrence 

interval would approach the number of station-years of 

record available.

Using Monte Carlo simulation techniques, the exceedance probability for rare flood events may be estimated by the simultaneous consideration of concurrent 

records of annual extremes which are dependent and 

identically distributed. The method for obtaining rareevent probabilities is to generate sets of normal multivariates equal in number and duration to the observed 

sets of records through the simulation model

to order the simulated record maxima and to compute 

the normal exceedance probabilities of the ordered 

maxima. This scheme of generating synthetic variates 

preserves the observed dependence in records.

In applying the procedure of estimating exceedance 

probabilities to regional flood-frequency problems it is 

assumed that a means of defining the region of identical 

distributions is at hand, that the estimate of exceedance 

probability is independent of the underlying statistical 

distribution of observed floods, and that the sample 

correlation matrix is a reasonable estimate of the true 

matrix. Indications are that the probability estimates 

are independent of the distribution of observed floods.

If the correlation coefficients for the records are all 

equal, the rare-event probabilities vary in a predictable 

manner with changes in numbers of records, period of 

record, and magnitude of the coefficient. If the correlation matrix is a mixture of coefficients, variations in the 

probabilities are a function not only of an equivalent 

mean coefficient, the number of records, and the period 

of record, but also of the deviation of the coefficients 

about the mean coefficient. Even though the probabilities associated with a general correlation matrix 

also appear to vary predictably with these variables, 

any general solution would not seem economically 

justified because of excessive computer costs (about 

Ink seconds per run). It would seem to be preferable in 

the interest of economy to estimate the rare-event exceedance probabilities for each set of data under 

analysis.

A set of criteria for selecting records having identical 

distributions has been identified in this report. Other 

sets of criteria for making this selection are possible, 

and the most definitive must be sought. The simulation 

program is equally applicable to many kinds of rare 

hydrologic and hydrometeorologic events such as precipitation, floods, and temperatures.

REFERENCES

Conover, W. J., 1965, A k-sample model in order statistics: 

Annals Math. Statistics, v. 36, no. 4, p. 1223-1235.

Conover, W. J., and Benson, M. A., 1963, Long-term flood 

frequencies based on extremes of short-term records, in 

Geological Survey research, 1962: U.S. Geol. Survey Prof. 

Paper 450-E, p. E159-E160.

Fiering, Myron B., 1968, Schemes for handling inconsistent 

matrices: Water Resources Research, v. 4, no. 2, p. 291-297.

International Business Machines, 1968, System/360 scientific 

subroutine package (360-CM-03X) version III, programmer's manual: White Plains, N.Y., International 

Business Machines, 454 p.

Kendall, M. G., and Stuart, Alan, 1963, Distribution theory 

Volume 1 of The advanced theory of statistics [2d ed.]: 

New York, Hafner, 433 p.

Matalas, N.C., 1967, Mathematical assessment of synthetic 

hydrology: Water Resources Research, v. 3, no. 4, p. 937- 

945.

421-897 O - 71 - 2
F10 STATISTICAL STUDIES IN HYDROLOGY

Morrison, D. F., 1967, Multivariate statistical methods: New

York, McGraw-Hill, 338 p. 

Ostle, Bernard, 1963, Statistics in research [2d ed.]: Ames, Iowa,

Iowa, Iowa State Univ. Press, 585 p. 

Thomas, C. A., Broom, H. C., and Cummans, J. E., 1963,

Magnitude and frequency of floods In the United States 

Part 13, Snake River basin: U.S. Geol. Survey WaterSupply Paper 1688, 250 p., 6 pis.

Zelen, Marvin, and Severs, M. A., 1965, Probability functions, 

in Abramowitz, Milton, and Stegren, I. A., eds., Handbook 

of mathematical functions: Natl. Bur. Standards Appl. 

Math. Ser. 55, p. 925-995.
COMPUTER PROGRAM
*********************

 

**  PROGRAM RAREVENT **

 

** ********************

* DISCLAIMER *

* *

ALTHOUGH PROGRAM HAS BEEN TESTED BY ITS CONTRIBUTOR, NO 

WARRANTY, EXPRESSED OR IMPLIED, IS MADE BY THE CONTRIBUTOR OR 

THE GOVERNMENT, AS TO THE ACCURACY AND FUNCTIONING OF THE 

PROGRAM AND RELATED PROGRAM MATERIAL, NOR SHALL THE FACT OF 

DISTRIBUTION CONSTITUTE ANY SUCH WARRANTY, AND NO 

RESPONSIBILITY IS ASSUMED BY THE CONTRIBUTOR OR THE 

GOVERNMENT, IN CONNECTION THEREWITH.

PURPOSE

COMPUTES ESTIMATE OF EXPECTED EXCEEDANCE PROBABILITY OF 

ORDERED COLUMN MAXIMA FROM N CROSS-CORRELATED, CONCURRENT 

RECORDS OF K EVENTS, WHERE THE RECORDS ARE SIMULATED

DESCRIPTION OF 

INPUT 

N

VARIABLES

- NUMBER OF RECORDS. MUST BE 20 OR LESS. 

CAPACITY OF THE PROGRAM CAN BE INCREASED FROM N = 20 

AND K = 100 BY CHANGING DIMENSION STATEMENTS. 

ASSUMING K = 100, N CAN BE AS GREAT AS 100 WITH 

PRESENT CONFIGURATION OF THE U.S. GEOLOGICAL 

SURVEY COMPUTER SYSTEM.

K - NUMBER OF EVENTS IN A RECORD. MUST BE 100 OR LESS.

IX - SEED FOR NORMAL RANDOM NUMBER GENERATOR. FORMAT 

IS 14.

RO - CORRELATION MATRIX {READ IN BY SUBROUTINE MTRXIN)

IS - CODE FOR STORAGE MODE 

0 - GENERAL

1 - SYMMETRIC

2 - DIAGONAL 

FORMATS (PER DATA SET)

N,K 

RO

OUTPUT 

IER

- 7F10.0, LISTING BY ROWS, STARTING WITH 

DIAGONAL ELEMENT AND CONTINUING WITH 

ELEMENTS TO ITS RIGHT.

- II (INDICATES END OF DATA SET)

ERROR CODE FOR INPUT DATA CARDS OF CORRELATION

MATRIX

0 - NO ERROR

1 - INCORRECT NUMBER OF CARDS 

N - NUMBER OF RECORDS 

K - NUMBER OF EVENTS IN RECORD 

RO - CORRELATION MATRIX 

Bl - DIAGONAL MATRIX CONTAINING SOUAREROOTS OF

EIGENVALUES

R - EIGENVECTOR MATRIX 

NSD - LEAST NUMBER OF SIGNIFICANT

CHARACTERISTIC EQUATION. IF

CORRELATION MATRIX 

P(J)- EXPECTED EXCEEDANCE PROBABILITY FOR JTH-ORDFR

COLUMN MAXIMUM

DIGITS IN ELEMENT OF 

0, CHECK STRUCTURE OF

F12
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA F13

SUBROUTINES AND FUNCTION SUBPROGRAMS 

MTRXIN 

MTRXOT 

EGN 

GMPRD 

MAX 

WIN 

LISD

METHOD

REFFR TO P.H.CARRIGAN, JR., 1971, A FLOOD FREQUENCY RELATION 

BASED ON REGIONAL RECORD MAXIMA, U. S.GEOL.SURVEY PROF. 

PAPER 434-F

REMARKS

EXECUTION TIME ABOUT 2NK SECONDS PER DATA SET. FORT PLUS

LKEO TIME LESS THAN 15 SECONDS

MAXIMUM LINES OF OUTPUT = 800 + 205*(NO. OF DATA SETS)

DIMENSION X(20,100),Y(20),P(20),JUMK(20) #1 

DOUBLE PR EC IS I ON 01 ,02 ,03 ,D4 ,D5 ,06 ,U ,RO ( 210 ) ,R (400 ) , B ( 400 ) , Bit 400 )#1 

DOUBLE PRECISION XI!20) , E(20 ) ,RBC(210) #1 

DATA Dl,02,03/0.049867347,0.0211410061,0.0032776263/ #1 

04,05,06/0.0000380036,0.0000488906,0.0000053830/ #1

#1

#1

 EMPIRICAL EXCEEDANCE PROBABILITIES'/I3X,'FOR ORDER#1 

'-RECORD MAXIMA'/18X,'WITH',14,' EVENTS PER RECORD'#1

#1

#1 

//) #1

#1 

ORDER SIGNIFICANT DIGIT#1

DATA

101 FORMAT(14)

102 FORMAT(2I4) 

201 FORMAT('1',12X 

1ED SET OF ' , 12

202 FORMAT! 20X

203 FORMAT!'O'/' 

204 FORMAT! 0«/' '

205 FORMAT!'0',10X

'CORRELATION MATRIX'//) 

20X,'SQRT EIGENVALUE MATRIX 

20X,'EIGENVECTOR MATRIX'//) 

 THE LEAST NUMBER OF HIGH

IS !NSD) IN A CHFCK OF CHARACTERISTIC EQUATION IS«, I 3,'.'///) #1

206 FORMAT! «0« ,9X, 'ORDER',8X,8!12,10X)) #1

207 FORMAT(4X,'PROBABILITY 1 ,5X,8(F9.6,3X) ) #1

208 FORMAT!'0',10X,'INCORRECT NUMBER OF DATA CARDS FOR CORRELATION MAT#1 

1RIX'//) #1

209 FORMAT!«0',10X,'MATRIX OF CORRELATION COEFFICIENTS DOES NOT TEST T#1 

10 BE POSITIVE SFMIOEFINITE.(DETERMINANT LESS THAN ZERO.)'/ #1 

2' ',10X,'THIS MAY BE A RESULT OF INACCURACIES IN THE COMPUTER ALGO#1 

3RITHMS.'/' I ,10X,«***THIS PROGRAM WILL NOT OPERATE ON THIS DATA SE#1 

4T BEYOND STATEMENT 70***')' #1

READ!5,101 ) IX 

1 READ(5,102,END=500)N,K

READ IN RO

CALL MTRXIN(RO,N, IER ) 

WRITE(6,201)N,K

CHECK FOR ERROR IN INPUT

IF!IER)60,61,60

60 WRITF(6,208) 

GO TO 1

61 WRITF(6,202)

PRINT CORRELATION MATRIX

CALL MTRXOT(RO,N,1)

COMPUTE EIGENVECTORS AND EIGENVALUES OF RO

NSD=0

CALL EGN(RO,R,RBC,N,NSO)

CHECK FOR POSITIVE SEMIOEFINITE CORRELATION MATRIX

#1

#1

#1

#1

#1

#1

#1

#1

#1

#1

0000

0010

0020

0030

0040

0050

0060

0070

,0080

0090

,0100

0110

,0120

0130

,0140

0150

,0160

,0170

,0180

,0190

,0200

,0210

,0220

,0230

0240

,0250

,0260

,0270

,0280

,0290

,0300

,0310

#1.0320

,0330

,0340
F14 STATISTICAL STUDIES IN HYDROLOGY

DO 70 1 = 1, N

IF(RO( IJ) .GE.0.0) GO TO 70 

WRITE(6,209) 

GO TO 71 

70 CONTINUE

COMPUTE SOUAREROOTS OF EIGENVALUES AND SET OFF-DIAGONAL ELEMENTS 

OF RO TO ZERO. PUT ELEMENTS IN GENERAL STORAGE MODE.

DO 2 1=1, N 

DO 2 J=1,N 

IF(I-J)5,6,6

5 IJ=I+< J*J-J)/2 

GO TO 3

6 IJ=J+( 1*1-1 )/2

IFU.NE.J) GO TO 3 

R0( IJ)=DSORT(RO(IJ ) ) 

GO TO 4

3 R0( IJ)=0.0

4 K1 = N*( J-l ) + I 

2 BKK1)=RO( IJ )

PRINT SORT OF EIGENVALUES AND EIGENVECTOR MATRICES

WRITE(6,203)

CALL MTRXOT(B1,N,0)

WRITE(6,204)

CALL MTRXOT(R,N,0)

PRINT SIGNIFICANT DIGITS MESSAGE

IF(NSD.EO.O) GO TO 73

72 WRITF( 6,205)NSD

COMPUTE PRINCIPAL COMPONENT MATRIX B

73 CALL GMPRD(R,B1 ,B,N,N,N ) 

INITIALIZE P(J) TO ZERO

DO 7 J = 1,N 

7 P(J)=0.0

#1.0350

#1.0360

#1.0370

#1.0380

#1.0390

#1.0400

#1.0410

#1.0420

#1.0430

#1.0440

#1.0450

#1 .0460

#1.0470

#1.0480

#1.0490

#1 .0500

#1.0510

#1.0520

#1.0530

#1.0540

#1.0550

#1.0560

#1.0570

#1.0580

#1.0590

#1.0600

#1.0610

IN THE 1000 ITERATION LOOP ENDING IN STATEMENT 10, AN N BY K 

RECORD OF X(J,KK) IS SIMULATED (USING MATRIX B AND A NORMAL RANDOM 

NUMBER E), N COLUMN MAXIMA Y(J) ARE SELECTED FROM X(J,KK), 

Y(J) ARE PLACED IN DESCENDING ORDER ( J=l , 2 , . . . , N ) , NORMAL 

EXCEEDANCE PROBABILITY C FOR EACH Y(J) IS FOUND, AND THE 1000 

ITEM SUM OF C FOR EACH J, P(J), IS COMPUTED.

DO 10 M=l,1000 #1.0620 

GENERATE X(J,KK) ARRAY

DO 11 KK=1,K #1.0630 

DO 15 1 = 1, N #1.0640

SUBROUTINE GAUSS IS IMBEDDED IN MAIN PROGRAM BY FOLLOWING 

10 STATEMENTS

AA=0.0 #1.0650

DO 12 I 1 = 1,12 #1.0660

IY=IX*65539 #1.0670

IF{ IY)13,14,14 #1.0680

13 IY=IY+1+2147483647 #1.0690

14 YFL=IY #1.0700

YFL=YFL*0.4656613E-9 #1.0710

IX=IY #1.0720

12 AA=AA+YFL #1.0730

15 E( I )=AA-6.0 #1.0740
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA

FIND K VECTOR PRODUCTS OF XI - X(N,K)

CALL GMPRD(B,E,X1,N,N,1) #1.0750

DO 11 1=1,N #1.0760

J=I #1.0770

11 X(J,KK)=X1(I) #1.0780

FIND COLUMN MAXIMA

DO 20 J=1,N #1.0790

D=-1000.0 #1.0800

DO 21 KK=1,K #1.0810

CM = AMAX1(X(J,KK ),D) #1.0820

21 D=CM #1.0830

20 Y(J)=D #1.0840

ORDER COLUMN MAXIMA

NN=M-1 #1.0850

30 DO 31 J=1,NN #1.0860

IF((Y(J)-Y(J+1>).GE.0.0) GO TO 31 #1.0870

A=Y(J+1) #1.0880

Y(J+1)=Y(J) #1.0890

Y(J)=A #1.0900

31 CONTINUE #1.0910

NN=NN-1 #1.0920

IF(NM.GT.O) GO TO 30 #1.0930

COMPUTE NORMAL EXCEEDANCE PROBABILITY C OF EACH COLUMN MAXIMUM

Y(J). (REF. - EON. 26.2.19 NBS HNDBK (1965))

DO 10 J=1,N #1.0940

IF(Y(J) .GE.0.0) GO TO 40 #1.0950

U=-Y(J) #1.0960

GO TO 41 #1.0970

40 U=Y(J) #1.0980

41 IF(U.GT.5.0) GO TO 42 #1.0990

IF(U.LT.0.0001) GO TO 43 #1.1000 

G=l .0 + D1*U+D2*U#U + D3#U*U#IJ+D4#U##4 + D5*U#*5 + D6*U**6 #1 . 1010

C=0.5*G**(-16) #1.1020

GO TO 44 #1.1030

42 C=0.0 #1.1040

GO TO 44 #1.1050

43 C = 0.5 #1.1060

44 IF(Y(J) .GE.0.0) GO TO 10 #1.1070

C=C-1.0 #1.1080

SUMMATION OF C = P(J)

10 P(J)=P(J)+C #1.1090 

COMPUTE MEAN EXCEEDANCE PROBABILITY P(J) FOR JTH-ORDER MAXIMUM

DO 45 J=1,N #1.1100

45 P(J)=P(J)/1000.0 #1.1110

PRINT P(J)

JS=1 #1.1120

IF(N.LE.B) GO TO 50 #1.1130

JE = 8 #1..1140

GO TO 51 #1.1150

50 JE=N #1.1160

51 DO 52 J=JS,JE #1.1170

52 JUNK(J)=J #1.1180

WRITE(6,206)(JUNK(J),J=JS,JE) #1.1190

WRITE(6,207)(P(J),J=JS,JE) #1.1200

JS=JS+8 #1.1210

IF(JS.GT.N) GO TO 53 #1.1220

JE=JE+8 #1.1230

F15
F16 STATISTICAL STUDIES IN HYDROLOGY

IF(LE.LT.N) GO TO 51 #1.1240

JE = N #1.1250

GO TH 51 #1.1260

53 CONTINUE #11.1270

71 GO TO 1 #1.1280

500 STOP #1.1290

END #1.1300

SUBROUTINE MTRX IN ( A ,NORD , I ER ) #1.1310

PURPOSE

READS DATA ELEMENTS OF SYMMETRIC MATRIX FROM CARDS

DESCRIPTION OF VARIABLES 

A - INPUT MATRIX

NORD - ORDER OF INPUT MATRIX. MUST BE 20 OR LESS. 

IER - ERROR IN DATA CARDS CODE 

IER = 0 NO ERROR

= 1 INCORRECT NUMBER OF CARDS

LISTING OF DATA CARDS

FORMAT IS 7F10.0. COLUMNS 71-80 FOR IDENTIFICATION. DATA 

ELEMENTS PUNCHED BY ROW STARTING WITH DIAGONAL ELEMENT AND 

PROCEEDING TO THE RIGHT IN THAT ROW, WITH LISTING CONTINUED 

TO NEW CARDS IF NECESSARY. LAST CARD FOR LAST ROW IN MATRIX 

FOLLOWED BY A CARD WITH 9 PUNCHED IN COLUMN 1.

REFERENCE

ADAPTATION OF SUBROUTINE MATIN. REFER TO INTERNATIONAL 

BUSINESS MACHINES, 1968, SYSTEM/360 SCIENTIFIC SUBROUTINE 

PACKAGE (360-CM-03X) VERSION III, PROGRAMMER'S MANUAL. 

WHITE PLAINS, NEW YORK, INTERNATIONAL BUSINESS MACHINES, 

P.453

DOUBLE PRECISION A(1),CARD(8) #1.1320

1 FORMAT(7F10.0) #1.1330

IER=0 #1.1340

ICOLT=NORD #1.1350

IROCR=1 #1.1360

COMPUTE NUMBFR OF CARDS IN ROW

2 IRCDS = ( ICOLT-1)/7+l #1.1370 

SET UP LOOP FOR NUMBER OF CARDS IN ROW

DO 3 K=l , IRCDS #1.1380

READ<5,1)(CARD(I),1=1,7) #1.1390

L=0 #1.1400

COMPUTE COLUMN NUMBER FOR FIRST FIELD IN CURRENT CARD

JS=(K-1)*7+NORD-ICOLT+1 #1.1410 

JE=JS+6 #1.1420

SET UP LOOP FOR DATA ELEMENTS WITHIN CARD

DO 4 J=JS,JE #1.1430

IF(J-NORD)5,5,3 #1.1440

5 IF( IROCR-J)6,7,7 #1.1450

6 IJ=IROCR+(J*J-J )/2 #1.1460

GO TO 8 #1 .1470

7 IJ=J+(IROCR*IROCR-IROCR)/2 #1.1480

8 L=L+1 #1.1490

4 A( IJ)=CARD(L) #1.1500

3 CONTINUE #1.1510

IROCR=IROCR+1 #1.1520

IF(NORD-IROCR)9,10,10 #1.1530

10 ICOLT=ICOLT-1 #1.1540

GO TO 2 #1.1550
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA F17

9 READ(5,1)CARD(1) #1.1560

IF(CARD(l)-9.09)11,12,11 #1.1570

11 IER=1 #1.1580

12 CONTINUE #1.1590

RETURN #1.1600

END #1.1610

SUBROUTINE MTRXOT(A,NORD,IS> #1.1620

PURPOSE

PRINTOUT LISTING OF ANY SIZED MATRIX

DESCRIPTION OF VARIABLES 

A - OUTPUT MATRIX 

NORD - ORDER OF MATRIX 

IS - STORAGE MODE OF A

IS = 0 GENERAL

IS = 1 SYMMETRIC

IS = 2 DIAGONAL

SUBROUTINES REQUIRED 

LOG

REFERENCE

ADAPTATION OF SUBROUTINE MXOUT. REFER TO INTERNATIONAL 

BUSINESS MACHINES, 1968, SYSTEM/360 SCIENTIFIC SUBROUTINE 

PACKAGE (360-CM-03X) VERSION III, PROGRAMMER'S MANUAL. 

WHITE PLAINS, NEW YORK, INTERNATIONAL BUSINESS MACHINES, 

P.454 

DOUBLE PRECISION A(1),B(8) #1.1630

201 FORMAT(20X,7(3X t I3,10X)) #1.1640

202 FORMAT( '0' ,1IX,13,7(D16 .6)) #1.1650 

J=l #1.1660

1 LS=1 #1.1670

2 JNT=J+6 #1.1680 

IF(JNT-NORD)4,4 t 3 #1.1690

3 JIMT = NORD #1.1700

NUMBER COLUMNS

4 WRITE<6,201)(JCUR , JCUR=J,JNT) #1.1710

LE=LS+57 #1.1720

DO 5 L=LS,LE #1.1730

PRINT LINE

DO 6 K=l,7 #1.1740

KK=K #1.1750

JT=J+K-1 #1.1760 

CALL LOG (L,JT,IJ,NORD t NORD,IS) #1.1770

B(K)=0.0 #1.1780

IF(IJ)8,8,7 #1.1790

7 B(K)=A(IJ) #1.1800

IF LAST COLUMN, WRITE ROW'S NUMBER AND ELEMENTS

8 IF(JT-NORD)6,9,9 #1.1810

6 CONTINUE #1.1820

9 WRITE(6,202)L, (B(JW) ,JW=1,KK) #1.1830

IF END OF ROWS, CHECK FOR MORE COLUMNS

IF(NORD-L)10,10,5 #1.1840 

5 CONTINUE #1.1850

CHECK FOR MORE OUTPUT

LS=LS+58 #1.1860 

GO TO 2 #1.1870

#1.1880
F18 STATISTICAL STUDIES IN HYDROLOGY

C END DF COLUMNS, RETURN

#1.1890

10 IF(JT-NORDJll,12,12 #1.1900

11 J=JT+1 #1.1910 

GO TO 1 #1.1920

12 RETURN #1.1930 

END #1.1940 

SUBROUTINE EGN(A,R,B,N,NSD) #1.0000 

C PURPOSE

C COMPUTE EIGENVALUES, EIGENVECTORS, AND LEAST NUMBER OF 

C SIGNIFICANT DIGITS IN CHARACTERISTIC EQUATION OF A RFAL 

C SYMMETRIC MATRIX 

C DESCRIPTION OF VARIABLES

C A MATRIX FOR WHICH EIGENVECTORS AND EIGENVALUES ARE 

C DESIRED. MATRIX DESTROYED IN COMPUTATION AND RESULTANT 

C EIGENVALUES STORED IN DESCENDING ORDER IN DIAGONAL. 

C R RESULTANT EIGENVECTOR MATRIX 

C B - ORIGINAL MATRIX TRANSFERRED TO THIS AREA 

C N ORDER OF MATRICES

C NSD - IF EIGENVALUES ARE TO BE COMPUTED, NSD MUST BE ZERO AT 

C INPUT, RESULTANT IS LEAST NUMBER OF SIGNIFICANT DIGITS 

C IN CHARACTERISTIC EQUATION STORED IN NSD. IF NSD = 1, 

C EIGENVALUES AND LEAST NUMBER OF SIGNIFICANT DIGITS 

C ARE NOT COMPUTED. 

C

C FUNCTION SUBPROGRAMS REQUIRED 

C LISD 

C MAX 

C MIN 

C REMARKS

C THIS EXTENSION TO DOUBLE PRECISION FORM OF SUBROUTINE 

C EIGEN DEVELOPED BY MARSHALL HELLMAN, COMPUTER CENTER 

C DIVISION, U.S. GEOLOGICAL SURVEY. REFER TO INTERNATIONAL 

C BUSINESS MACHINES, 1968, SYSTEM/360 SCIENTIFIC SUBROUTINE 

C PACKAGE (360-CM-03X) VERSION III, PROGRAMMER'S MANUAL. 

C WHITE PLAINS, NEW YORK, INTERNATIONAL BUSINESS MACHINES, 

C P . 1 64 

C

DIMENSION Ad), R(l), B(l) #1.0010

DOUBLE PREC IS ION A,R,ANORM,ANRMX,THR,X,Y,SI NX,SINX2,COSX,COSX2,SIN#l.0020

1CS,RANGE,B,W,W1,Z,Q #1.0030

201 FORMAT( '0' ,10X,'NUMBER OF SIGNIFICANT DIGITS IN ELEMENT ',13,' OF #1.0040

1CHARACTERISTIC EOUATION /1IX,'MATRIX IS 0. MAGNITUDES OF ELEMENTS #1.0050

2COMPARED ARE ,D24.16/1IX,'AND ', D24.16,'. IF MAGNITUDE OF AN ELE#1.0060

3MENT IS NOT CLOSE '/11X,' TO ZERO, IT WOULD BE ADVISABLE TO STUDY #1.0070

^STRUCTURE OF CORRELATION MATRIX'//) #1.0080

MSDR=15 #1.0090

MSD=15 #1.0100

MV=NSD #1.0110

RANGF=1.0D-14 #1.0120

C

C COPY A INTO B 

C

NX2=N#(N+1)/2 #1.0130 

DO 1 I=1,NX2 #1.0140

1 B( I )=A(I) #1.0150 

C

C GENERATE IDENTITY MATRIX 

C

IF (MV-1) 2,5,2 #1.0160

2 IQ=-N #1.0170 

DO 4 J=1,N #1.0180 

IQ=IQ+N #1.0190 

DO 4 1 = 1,N #1.0200 

IJ=IO+I #1.0210 

R( IJ)=O.ODO #1 .0220 

IF (I-J) 4,3,4 #1.0230

3 R(IJ)=1.0DO #1.0240

4 CONTINUE #1.0250
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA 

COMPUTE INITIAL AND FINAL NORMS, ANORM AND ANORMX

5 ANORM=O.ODO #1.0260

DO 7 1 = 1,N #1.0270

DO 7 J=I,N #1.0280

IF (I-J) 6,7,6 #1.0290

6 IA=I+(J*J-J)/2 #1.0300 

ANORM=ANORM+A(IA)*A(IA) #1.0310

7 CONTINUE #1.0320

, IF (ANORM) 33,33,8 #1.0330

8 ANORM=DSQRT(2*ANORM) #1.0340 

ANRMX=ANORM#RANGE/DFLOAT(N) #1 .0350

INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR

IND=0 #1.0360

THR=ANORM #1.0370

9 THR=THR/DFLOAT(N) #1.0380

10 L=l #1.0390

11 M=L+1 #1.0400

COMPUTE SIN AND COS

12 MG=(M*M-M)/2 #1.0410

LQ=(L#L-L)/2 #1.0420

LM=L+MQ #1.0430 

IF (DABS(A(LM) )-THR) 26,13,13 #1.0440

13 IND=1 #1.0450

LL=L+LQ #1.0460

MM=M+MQ #1.0470

Z=A(LL) #1.0480

Q=A(MM) #1.0490

X=(Z-0)/2 #1.0500 

Y = -A( LM)/DSQRT(A(LM)#A(LM)+X#X) #1 .0510

IF (X) 14,15,15 #1.0520

14 Y=-Y #1.0530

15 SINX=Y/DSQRT(2*(1+(DSQRT(1-Y#Y) ))) #1.0540

SINX2=SINX*SINX #1.0550

COSX2=1-SINX2 #1.0560

COSX=DSORT(COSX2) #1.0570

SINCS=SINX*COSX #1.0580

ROTATE L AND M COLUMNS

ILQ=N*(L-1) #1.0590

IMQ=N*(M-1) #1.0600

DO 25 1=1,N #1.0610

IO=(I*I-I)/2 #1.0620

IF (I-L) 16,23,16 #1.0630

16 IF U-M) 17,23,18 #1.0640

17 IM=I+MO #1.0650

GO TO 19 #1.0660

18 IM=M+IO #1.0670

19 IF (I-L) 20,21,21 #1.0680

20 IL=I+LO #1.0690

GO TO 22 #1.0700

21 IL=L+IO #1.0710

22 X = A( IL)*COSX #1.0720

W=A(IM)*SINX #1.0730

Y=A(IL)*SINX #1.0740

W1=-A(IM)*COSX #1.0750

A( IM) = Y-W1 #1.0760

A(IL)=X-W #1.0770

23 IF (MV-1) 24,25,24 #1.0780

24 ILR=ILQ+I #1.0790

IMR=IMO+I #1.0800

Z=R(ILR)*COSX #1.0810

0=R( IMR)*SINX #1.0820

X=Z-0 #1.0830

Z=R(ILR)*SINX #1.0840

F19
F20 STATISTICAL STUDIES IN HYDROLOGY

Q=-R( IMR)#COSX

R(IMR)=Z-(0

R(ILR)=X

25 CONTINUE

X=2*A( LM)*SINCS

Z=A(LL)*CnsX2

0=-A( MM)*SINX2

Y=Z-0

Z=A(LL)*SINX2

Q=-A(MM) *COSX2

W=Z-0

W1=-W

A(LM)=0

A(LL)=Y-X

A(MM)=W+X

TESTS FOR COMPLETION 

TEST FOR M = LAST COLUMN

26 IF (M-N) 27,28,27

27 M=M+1

GO TO 12

TEST FOR L = SECOND FROM LAST COLUMN

28 IF (L-(N-l)) 29,30,29

29 L=L+1

GO TO 11

30 IF (IND-1) 32,31,32

31 IND=0

GO TO 10

COMPARE THRESHOLD WITH FINAL NORM

32 IF (THR-ANRMX) 33,33,9

SORT EIGENVALUES AND EIGENVECTORS

33 IQ=-N

DO 37 1=1 ,N

H)=IO+N

LL=I+(I*I-I)/2

JO=N*( 1-2)

DO 37 J=I ,N

JQ=JO+N

MM=J+( J*J-J)/2

IF ( A( LL)-A(MM) ) 34,37,37

34 X=A( LL)

A(LL)=A(MM)

A(MM)=X

IF (MV-1) 35,37,35 

35D036K=1,N

ILR=IO+K

IMR=JO+K

X=R( ILR)

R( ILR)=R( IMR)

36 R( IMR)=X

37 CONTINUE

DETERMINE LEAST NUMBER OF SIGNIFICANT DIGITS

IF (MV-1) 38,48,38

38 CONTINUE

COMPUTE MATRIX ELEMENTS OF PRODUCT OF ORIGINAL MATRIX AND CHARACTERISTIC VECTORS, X, AMD OF EIGENVALUES AND CHARACTERISTIC 

VECTORS, Z.

#1.0850

#1.0860

#1.0870

#1.0880

#1.0890

#1.0900

#1.0910

#1.0920

#1.0930

#1.0940

#1.0950

#1.0960

#1.0970

#1.0980

#1.0990

#1.1000

#1.1010

#1.1020

#1.1030

#1.1040

#1.1050

#1.1060

#1.1070

#1.1080

#1.1090

#1.1100

#1.1110

#1.1120

#1.1130

#1 .11^0

#1.1150

#1.1160

#1.1170

#1.1180

#1.1190

#1.1200

#1.1210

#1.1220

#1.1230

#1.1240

#1.1250

#1.1260

#1.1270

#1 .1280

#1.1290

#1.1300

#1.1310
FLOOD-FREQUENCY RELATION BASED ON REGIONAL RECORD MAXIMA F21

39

41

42

49

43

44

45

46

47

48

DO 45 12 = 1 ,N

Y = A( !2*( 12 + 1 )/2)

I1=N*(12-1)

DO 43 L=l f N

X = 0

DO 42 M = l ,N

IF (M-L) 39,40,40

I=L*(L-l)/2+M

GO TO 41

I=M*(M-l)/2+L

K=I 1 + M

X=X+B(I)*R(K)

Z=Y*R( Il + L)

DETERMINE LEAST NUMBER OF SIGNIFICANT DIGITS IN X OR Y FOR EACH 

ELEMFNT f)F MATRIX V

NSD=LISD(X,Z,15)

IF(NSD)43,49,43

I1L=I1+L

WRITE(6,201 ) I1L,X,Z 

MSD=MIN(NSD,MSD) 

MSDR=MIM <MSD,MSDR)

K1=I2*( I2 + D/2 + 1

IF ( I2-N) 44,46,46

A(K1)=MSD

MSD=15

DO 47 1=1,N

K=I*(1+1)/2

L = K+1

B(K)=A(L)

B(NX2)=MSD

NSD=MSDR 

CONTINUE

RETURN 

END

FUNCTION LISD(X,Y,N) 

PURPOSE

TO FIND SMALLER NUMBER OF SIGNIFICANT DIGITS IN TWO NUMBERS

X AND Y, IF NUMBER IS LESS THAN N

REMARKS

PROGRAMMED BY MARSHALL HELLMAN, COMPUTER CENTER DIVISION, 

U.S. GEOLOGICAL SURVEY

PRECISION X,Y,EPS,T

#1.1320

#1.1330

#1.1340

#1.1350

#1 .1360

#1.1370

#1 .1380

#1.1390

#1.1400

#1.1410

#1 .1420

#1.1430

#1.1440

#1 .1450

#1 .1460

#1.1470

#1.1480

#1.1490

#1 .1500

#1.1510

#1 .1520

#1.1530

#1.1540

#1.1550

#1 .1560

#1.1570

#1 .1580

#1.1590

#1 .1600

#1.1610

#1 .1620

#1.1630

#1.1640

DOUBLE

LISD=0

IF(DABS(Y).GE.DABS(X))

T = X

IF(T.EO.O.O) GO TO 2

GO TO 3

T = Y

IF(T.EO.O.O) GO TO 2

EPS = DLOG10(UAB5 (T) )

IX=EPS

IF(EPS.GE.0.0) GO TO 4

IF(EPS.GE.IX) GO TO 4

IX=IX-1

IF(Y.NE.X) GO TO 5

LISD=N

GO TO 6

T=DLOG10(DABS(Y-X))

IT = T

IF(T.GE.O.O) GO TO 7

IF(T.GE.IT) GO TO 7

IT=IT-1

LISD=IX-IT-1

LISD=MAX(LISn,0)

RETURN

END

GO TO 1

#1.1650

#1.1660

#1.1670

#1.1680

#1.1690

#1 .1700

#1.1710

#1 .1720

#1.1730

#1.1740

#1.1750

#1 .1760

#1.1770

#1 .1780

#1.1790

#1 .1800

#1.1810

#1 .1820

#1 .1830

#1 .1840

#1.1850

#1 .1860

#1.1870

#1 .1880

#1.1890
F22 STATISTICAL STUDIES IN HYDROLOGY

FUNCTION MAX(M,N) #1.1900 

PURPOSE

TO DETERMINE IF INTEGER M IS LARGER THAN SELECTED INTEGER N

RFMARKS

PROGRAMMED BY MARSHALL HELLMAN, COMPUTER CENTER DIVISION, 

U.S. GEOLOGICAL SURVEY

IF(M.GT.N) GO TO 1 #1.1910

MAX=N #1.1920

RETURN #1.1930

MAX=M #1.1940

RETURN #1.1950

END #1.1960

FUNCTION MINIM,N) #1.1970 

PURPOSE

TO DETERMINE IF INTEGER M IS SMALLER THAN SELECTED INTEGER N

RFMARKS

PROGRAMMED BY MARSHALL HELLMAN, COMPUTER CENTER DIVISION, 

U.S. GEOLOGICAL SURVEY

IF(K.LE.N) GO TO 1 #1.1980

MIN=N #1.1990

RETURN #1.2000

MIN=M #1.2010

RETURN #1.2020

END #1.2030

SUBROUTINE GMPRD,( A , B ,R , N ,M , L ) #1.2040

PURPOSE

PRODUCT OF TWO GFMERAL-MODF-STORAGF MATRICES 

DESCRIPTION OF VARIABLES

A - FIRST INPUT MATRIX, N ROWS, M COLUMNS

B - SECOND INPUT MATRIX, M ROWS, L COLUMNS

R - OUTPUT MATRIX, N ROWS, L COLUMNS

REFERENCE

INTERNATIONAL BUSINESS MACHINES, 1968, SYSTEM/360 SCIENTIFIC 

SUBROUTINE PACKAGE (360-CM-03X) VERSION III, PROGRAMMER'S 

MANUAL. WHITE PLAINS, NEW YORK, INTERNATIONAL BUSINESS 

MACHINES, P.99

DOUBLE PRECISION A{1 ) ,B(1 ) ,R(1 ) #1.2050

IR=0 #1.2060

IK=-M #1.2070

DO 1 K=1,L #1.2080

IK=IK+M #1.2090

DO 1 J=1,N #1.2100

IR=IR+1 #1.2110

JI=J-N #1.2120

IB=IK #1.2130

RUR)=0.0 #1.2140

DO 1 1 = 1, M #1.2150

JI=JI+N #1.2160

IB=IB+1 #1.2170 

R( IR)=R(IR)+A(JI)*B(IB) #1.2180

RETURN #1.2190 

END

U. S. GOVERNMENT PRINTING OFFICE : 1971 O - 421-897