Flood Frequency Analysis

Flood Frequency Analysis

Log-Pearson Type III, Expected Moments, Low-Flood Test and Confidence Limits

Scope

The tool estimates the flood frequency curve of a gaged site from its record of annual peak flows, by the two federal guidelines of the United States: Bulletin 17C (U.S. Geological Survey Techniques and Methods 4-B5, 2018, version 1.1 of May 2019) and its predecessor Bulletin 17B (1982). Both fit the Pearson Type III distribution to the base-10 logarithms of the peaks by the method of moments; they differ in how they treat low floods, censored peaks and historical information, and in the confidence limits. Every relation on the page comes from the two bulletins, which are U.S. government works, and from the public U.S. Geological Survey R package MGBT for the p-value of the low-flood test. The work of Cohn, Stedinger, Hirsch and Hosking that the bulletins rest on is cited, not reproduced.

Entering a record

Each line of the peaks box is one water year: the year and the peak flow. Bulletin 17C describes what is known about a year by a flow interval, and the box accepts the same forms: year lower upper for a flood known within a range, year 0 base or year <base for a peak that stayed below a gage base (a crest-stage gage that was not reached), year >level for a flood known only to have exceeded a level, and 0 for a year without flow. The perception thresholds box is optional; each line gives a range of years and the smallest peak that would have been recorded in them. A range that holds a peak in fewer than half of its years is taken as a historical or broken-record period: its years without a peak are known to lie below the threshold, and that knowledge lengthens the record. A range over gaged years is a gage base. The flows are plain numbers in the unit chosen beside the box; only the check flow goes through the unit dictionary.

The distribution and the frequency factor

With X = log₁₀ Q, mean M, standard deviation S and skew G, the flow of annual exceedance probability p is Q_p = 10^{M + K S}. K is the standardized Pearson III quantile. For G ≠ 0 the distribution is a shifted gamma distribution with shape α = 4/G², scale β = S G/2 and location τ = M − 2S/G (Bulletin 17C equations 8 to 10), so K = (y − α)/√α with y the quantile of the incomplete gamma function at probability 1 − p for positive skew, and the mirror image for negative skew. The page solves this exactly and reproduces the K table of Bulletin 17B appendix 3 to its five decimals. The Wilson–Hilferty approximation K ≈ (2/G)[(1 + G z/6 − G²/36)³ − 1] is listed beside it; between annual exceedance probabilities of 99 and 0.2 percent it stays within 0.01 of the exact value for skews up to ±0.5, within 0.035 up to ±1 and within 0.10 up to ±2.

Bulletin 17C: low floods

Small annual peaks can bend the fitted curve at the flood end although they say nothing about large floods. Bulletin 17C screens them with the Multiple Grubbs–Beck test. The logarithms of the gaged peaks are sorted; for the k-th smallest the statistic ω_k = (X_[k] − μ_k)/σ_k is formed with the mean and standard deviation of the observations larger than it, and p(k; n) is the probability of a statistic that small for the k-th smallest of n normal variates. The page computes p(k; n) by the method of Cohn and others (2013) as coded in the MGBT package: an integral over the quantile of the k-th order statistic of a noncentral t probability whose parameters come from the moments of a normal sample truncated below. A sweep from the median outward stops at the first k with p below 0.005, a sweep from the smallest inward stops at the first k with p of 0.10 or more, and the larger k counts. The threshold is the next larger peak. Zero flows and peaks below a gage base count as the smallest observations; peaks of a historical period stay out of the test.

A threshold entered by hand replaces the test. The bulletin allows this where the test does not serve the data, provided the reason is documented. One of the published examples needs it: for Orestimba Creek (1932 to 2013) Bulletin 17C reports 30 low floods below 782 ft³/s, while this page finds p(38; 82) = 0.0049, just under 0.005, and therefore 38 low floods below 1,130 ft³/s. A Monte Carlo evaluation of the same probability gives 0.0043, so the difference lies in the numerical evaluation inside PeakFQ 7.1 at a borderline value, not in the rule. The example record loads with 782 entered, which reproduces the bulletin's moments and flows.

Bulletin 17C: the Expected Moments Algorithm

Peaks below the low-flood threshold are recoded as lying between zero and the threshold, and the lower perception threshold of every year is raised to it. The moments then come from an iteration (equations 7-1 to 7-9). A known logarithm contributes X, (X − M)² and (X − M)³ as it is. A censored year contributes the expected values of the same three quantities over its interval, computed from the Pearson III distribution of the current fit. The sums over the known values carry the bias factors n/(n − 1) and n²/((n − 1)(n − 2)), n the number of known values; the sums are divided by the total number of years. The new moments define a new distribution and the step is repeated until the moments change by less than 10⁻¹⁰. Without censored years the result is the ordinary sample mean, standard deviation and skew. The skew is not allowed below −1.41.

For skews up to ±1 the expectations are evaluated by Gauss–Legendre quadrature of the standardized density, which stays accurate as the skew approaches zero, where the incomplete gamma form of equation 7-9 loses digits; beyond ±1 the incomplete gamma form is used, because the density is no longer smooth at its bound.

Weighted skew

The station skew of a short record is uncertain, so both bulletins weight it with a regional skew G_r inversely to the mean square errors: G_w = (v_r G_s + v_s G_r)/(v_r + v_s). The page takes the mean square error of the station skew from Bulletin 17B equation 6 with N the number of years in the analysis. Bulletin 17C says the value is computed by the Expected Moments Algorithm; the equation 6 value with the full record length reproduces the mean square errors that PeakFQ prints on the figures of the bulletin's examples (0.101, 0.132 and 0.186), and it is what the page uses. With a weighted or regional skew the mean and standard deviation are refitted with the skew held fixed. The regional skew and its mean square error come from a regional study; the map of Bulletin 17B (Plate I, mean square error 0.302) is superseded in most states.

Confidence limits

For Bulletin 17C the page computes a first-order covariance matrix of the three moments. The estimating equations of the algorithm are linearized about the fitted distribution, each year being observed only when its flow exceeds its perception threshold; this is the construction the bulletin outlines in equations 7-22 to 7-30. The variance of log Q follows from the derivatives of the quantile with respect to the moments (equation 7-37), a regional skew adds its own mean square error, and the limits are log Q ± z √Var (equation 7-31).

These limits are not those of PeakFQ. The bulletin replaces the simple interval by an adjusted one (equations 7-33 and 7-34) whose terms are given only by reference to Cohn and others (2001) and to an unpublished revision of 2015, and they are not built here. Compared with the tables of the bulletin's examples, the variances computed here agree within about 20 percent for annual exceedance probabilities of 4 percent and larger and reach up to twice the printed value at 0.2 percent; the printed limits are not symmetric in log Q and are wider on the upper side. Where limits must agree with the bulletin, use PeakFQ or HEC-SSP.

For Bulletin 17B the limits are those of its appendix 9, K^U,L = [K ± √(K² − a b)]/a with a = 1 − z²/(2(N − 1)) and b = K² − z²/N. They treat the skew as known and are narrower at small probabilities.

Bulletin 17B procedure

The second method follows the flow chart of Bulletin 17B appendix 12. The statistics of the nonzero systematic peaks are computed; a Grubbs–Beck test at the 10 percent level gives the thresholds X̄ ± K_N S, with K_N from appendix 4 (the table is carried for N from 10 to 149). High outliers are tested first when the station skew exceeds +0.4 and low outliers first when it is below −0.4. Low outliers are deleted. When a historical period is described in the thresholds box, high outliers and historical peaks are the largest events of that period and the remaining peaks are weighted by W = (H − Z)/(N + L) (appendix 6, equations 6-1 to 6-4); without one, high outliers stay in the record. When zero flows, below-base peaks or low outliers have been removed, the conditional probability adjustment of appendix 5 multiplies the exceedance probabilities by the share of years above the truncation level and replaces the statistics by synthetic ones from the adjusted flows at 0.01, 0.10 and 0.50 (equations 5-3 to 5-5). The bulletin reads these three flows from a plot; the page computes them, so the synthetic skew can differ from the printed examples in the second decimal. The bulletin's examples also round the skew to the nearest tenth to use the K table; the page does not.

Plotting positions and the probability plot

Without thresholds the i-th largest of n peaks is plotted at p = (i − a)/(n + 1 − 2a), a = 0 (Weibull), 0.40 (Cunnane), 0.44 (Gringorten) or 0.50 (Hazen). With perception thresholds the Bulletin 17C method uses the threshold-exceedance positions of Hirsch and Stedinger (appendix 5, equations 5-1 to 5-5): the floods above each threshold share the exceedance probability of that threshold, estimated from the years in which such a flood would have been seen, and the floods below share the rest. The Bulletin 17B method uses the weighted order numbers of its equations 6-6 to 6-8. The plot's horizontal axis is the standard normal variate of the non-exceedance probability, which is what probability paper does.

Comparison fits and risk

The log-normal, Gumbel and GEV curves are for comparison. The Gumbel and GEV parameters come from the first three L-moments of the flows by Hosking's formulas (scale λ₂/ln 2 for the Gumbel; for the GEV the shape k = 7.8590 c + 2.9554 c² with c = 2/(3 + τ₃) − ln 2/ln 3). They use the known nonzero peaks only. The risk that a flow of annual exceedance probability p is exceeded at least once in n years is 1 − (1 − p)ⁿ. The expected probability adjustment of Bulletin 17B appendix 11 is not applied.

What the page reproduces

Bulletin 17C appendix 10: the moments and the eight flows of examples 2 (Orestimba Creek, with the threshold entered), 3 (Back Creek), 5 (Bear Creek) and 6 (Santa Cruz River) to every printed digit, the low-flood counts of examples 3, 5 and 6 with the thresholds of examples 3 and 6 (2,000 and 380 ft³/s; the bulletin prints none for example 5), and the moments of example 1 (Moose River). For example 1 the flows printed in table 10-5 are reproduced with the station skew 0.397; the text of the bulletin says they were computed with the weighted skew 0.421, which gives flows about half a percent higher at the 1 percent flood (4,985 against 4,957 ft³/s), as the default case of this page shows. Examples 4 and 7 (paleoflood records run with PeakfqSA) were not checked. Bulletin 17B appendix 12: the statistics, outlier thresholds, weights, frequency factors, flows and confidence limits of examples 1 to 4, with the differences described above where the bulletin reads a plot.

Limits

The page analyses one station. It does not estimate a regional skew, extend a record from a nearby station, weight independent estimates, separate mixed populations or adjust for regulation, urbanization or trends; Bulletin 17C discusses each. At least five known nonzero peaks are required, and the bulletin advises caution below ten. Plain flows are not converted when the unit system is switched.