Frequency Analysis

October 13, 2017 | Pengarang: hajihuzefa | Kategori: Estimator, Bias Of An Estimator, Mean Squared Error, Normal Distribution, Skewness
Share Benamkan


Penerangan Ringkas

Frequency Analysis of Extreme Events...

Penerangan

~

.

~) 1!;,, f';'

~,.

CHAPTER FREQUENCY

18 ANALYSIS

EXTREME

OF

EVENTS

if" ~i, t';: ?,r; ~ .';

~~ i;

Jery R. Stedinger Schoolof (:'ivil and Envir~nme.ntal Engineering

ri!J; "', i; ~w ~:, ~ ~ \;;; ~" ~, I!& ~

~l; ~ ~; " ~1 ~~ ~(1 ~'" ~~, I~i:' ~'c ~. ~:j ~

Cornell

University

Ithaca,

New

York

M.

Vogel

Richard

Department of Civil Engineering Tufts University . Medford,

Massachusetts

Ef " F

f I G . OU OU aeorglou Department of Civil and Mineral Engineering University of Minnesota Minneapolis, Minnesota I

I: Ir ;,:' ~ B~' ~; ~l 18. 1 INTRODUCTION TO FREQUENCY ~[ ANAL YSIS ~' ~~ !1 Extreme rainfall events and the resulting floods can take thousands oflives and cause ~ billions of dollars in damage. Rood plain management and designs for flood control ~ .'orks, reservoirs, bridges, and other investigations ne~d to reflect the likelihood or ~f probability of such events. Hydrologic studies also need to address the impact of il~ unusually low stream flows and pollutant loadings because of their effects on water ~ ('C quality and water supplies. 1-', ,~ ~, TIlt Basic Problem. Frequency analysis is an information problem: if one had a ~ sufficiently long record of flood flows, rainfall, low flows, or pollutant loadings, then ~ .frequency distribution for a site could be precisely detennined, so long as change ~f over time due to urbanization or natural processes did not alter the relationships of ~~ concern. In most situations, available data are insufficient to precisely define the risk ~~ or large floods, rainfall, pollutant loadings, or low flows. This forces hydrologists to i.~ U5Cpractical knowledge of the processes involved, and efficient and robust statistical i~):i 1"

--~

18.2

CHAPTEREIGHTEEN

techniques, to develop the best estimates of risk that theycan."s These techniquesare generally restricted, with I a to 100 sample observations to estjmate events exceeded with a chance of at least I in lOO, corresponding to exceedance probabilities of I percent or more. In some cases,they are used to estimate the rainfall exceededwith a chance of I in 1000, and even the flood flows for spillway design exceeded with a chance of I in 10,000. The hydrologist should be aware that in practice the true probability distributions of the phenomena in question are not known. Even if they were, their functional representation would likely have too many parameters to be of much practical use. The practical issue is how to select a reasonable and simple distribution to describe the phenomenon of interest, to estimate that distribution's parameters, and thus to obtain risk estimates of satisfactory accuracy for the problem at hand. Common Problems. The hydrologic problems addressedby this chapter primarily deal with the magnitudes of a single variable. Examples include annual minimum 7-day-averagelow flows, annual maximum flood peaks, or 24-h maximum precipitation depths. These annual maxima and minima for successiveyears can generally be considered to be independent and identically distributed, making the required frequency analyses straightforward. In other instances the risk may be attributable to more than one factor. flood risk at a site may be due to different kinds of events which occur in different seasons,or due to risk from several sources of flooding or coincident events, such as both local tributary floods and large regional floods which result in backwater flooding from a reservoir or major river. When the magnitudes of different factors are independent,a mixture model can be used to estimate the combined risk (see Sec. 18.6.2). In other instances, it may be necessaryor advantageous to consider all events that exceeda specified threshold because it makes a larger data set available, or because of the economic consequencesof every event; such partia/ duration series are discussedin Sec. 18.6.1.

18.1.1

Probability

Concepts

Let the upper caseletter X denote a random variab/e, and the lower case letter x a possible value of x. For a random variable X, its cumulative distribution function (cdf), denoted F x -0.5)

whenK> 0,X< ( ,+~); Weibull

Jlx- ~ + O.5772a

x- ~

--a-

[

(-a-x-~ )]

-exp

(

-exp

fX ( ,+~)

Pareto

fxO

.ux=a r ( I +i)

FX-0.33) .Here y -In

r( l+i)J}

a

.Ux='+-

a

1-

O"k=p2

"x-2

fx 0, where X(I) and X(n)are, respectively, the largestand smallest observed values; Xmedian is the sample medium equal to X(k+I) for odd sample sizesn = 2k + 1, and 1(X(k)+ X(k+I» for even n = 2k. [When X(I) + X(n)-2xmedian< 0, the formula provides an upper bound so that In( ~ -x) would be normally distributed.] Given ~, one can estimate)1yand cr}by using the sample mean and variance of Yi = In (xi -~), or wi = log (xi -~). The quantile-lower-bound estimator's per. fonnance with the resultant sample estimators of)1y and 0-} is better than method-of. moments estimators and competitive with maximum likelihood estimators.61.126 For the two-parameter and three-parameterlognormal distribution, the secondL moment is A.2= exp ~)1y+ -1 ) erf~ ~ ) = 2 exp ~)1y + -1 ) L (f>~-Ji) -~J

(18.2.12)

The following polynomial approximates, withiaO.0005 forIT31-1(I -qi). Figure 18.3.1 illustrate: use of lognormal paper with Blom's plotting positions. For the Gumbel distribution, G-l(1

-qi)

Thus a plot of X(i) versus G-I(I Gumbel variate

= ~ -a

-qi) Yi = -In

In [-In

(1-

qi)]

(18.3.8:

is identical to a plot of X(i) versus the reduceG . [-In

(I -qi)]

(18.3.9)

It is easy to construct probability paper for the Gumbel distribution by plotting X(i)as a function of Yi ; the horizontal axis can show the actual values of y or, equivalently, the associated qi, as in Fig. 18.3.1 for the lognormal distribution. Special probability papers are not available for the Pearson type 3 or log Pearson type 3 distributions because the frequency factors depend on the skew coefficient. However, for a given value for the coefficient of skewness y one can plot the observation X(i) for a P3 distribution, or log (X(i» for the LP3 distribution, versus the frequency factors Kp(Y) defined in Eq. (18.2.29) with Pi = 1- qi. This should yield a straight line except for sampling error if the correct skew coefficient is employed. Alternatively for the P3 or LP3 distributions, normal or lognormal probability paper is often used to compare the X(i) and a fitted P3 distribution, which plots as a curved

TABLE 18.3.2

Generation of Probability Plots for Different Distributions

Normal probability paper. Plot x(;) versus Zp,given in Eq. ( 18,2.3), where p; = 1 -q;. Blom's formula (a = 318) provides quantile-unbiased plotting positions. Lognormal probability paper. Plot ordered logarithms log [x(;J versus Zp. Blom's formula (a = 318) provides quantile-unbiased plotting positions. Exponential probability paper. Plot ordered observations x(;) versus f. -In (q;)IP or just -In (qJ. Gringorten's plotting positions (a = 0.44) work well. Gumbel and Weibull probability paper. For Gumbel distribution plot ordered observations x(;) versus f.-a In [-In (I-q;)] or just Yj=-ln [-In (l-q;)]. Gringorten's plotting positions (a = 0.44) were developed for this distribution. For Weibull distribution plot In [x(jJ versus In (a) + In [-In (qJ]lk or just In [-In (qj)]. (See Ref. 154.) GEV distribution. Plot ordered observationsx(i) versusf. + (a/K) { I -[-In ( I -qj)]K, or just (IlK) {1 -[-In (1 -qj)]K}. Alternatively employ Gumbel probability paper on which GEV will be curved. Cunnane's plotting positions (a = 0.4) are reasonable.s2 Pearson type 3 probability paper. Plot ordered observations x(j) versus Kpf(Y)'where pj = I qj. Blom's formula (a = 318) is quantile-unbiased for normal distribution and makessense for small Y. Or employ normal probability paper. (See Ref. 158.) Log Pearson type 3 probability paper. Plot ordered logarithms log [-\"(;Jversus Kp,(Y)where p; = 1 -qj. Blom's formula (a = 318) makes sensefor small Y.Or employ lognormal probability paper, (See Ref. 158.) Uniform probability paper. Plot x(j) versus I -q;, where qj are the Weibull plotting positions (a = 0). (See Ref. 154.)

18.26

CHAPTEREIGHTEEN

mal distribution G-I(I

-qJ

= ,u + a -1(1- qJ

(18.3.7

Thus, except for intercept and slope, a plot of the observations X(i) versus G-I [ 1 -qi is visually identical to a plot of X(i) versus -1(1 -qi). The values of qi are after printed along the abscissa or horizontal axis. Lognormal paper is obtained by using ( log scale to plot the ordered logarithms log (x(j» versus a normal-probability scale which is equivalent to plotting log (x(;y verus -1(I -qj). Figure 18.3.1 illustrate: use of lognormal paper with Blom's plotting positions. For the Gumbel distribution, G-l(1

-qj)

Thus a plot of X(i) versus G-I(I Gumbel variale

= -1

relatively

given

and

in

of

ranges.

By

complete-sample

the

to

the

Eq,

samplc

filling

in

estimators

insensitive

lognormal

is

estimators

medians

using

(18.6.14)

for

standard

as

then

made

a

be

in

a

(18.6.13)

constructed,

calculated,

procedure

observations

can

Ifregression

fori=r+

and

[%(r)]

'Ti

with

by

observation

a

is

To

is

be

be

~tatistics

well-descnbed

log

largest

to

as

-q;)

observations,

shown

IOm+szp

-r/n).

sample

can

small

=

IOY(i)

-1(1

...~

r

quantile

observations

complete

variance

missing

of

m+

Once

IS

positions.

the

(1

missing

%(i)=

where

~

the

plotting

of

probability

one

[%(1)]

where

Xp

tics,

data

log

...,

estimator

been

various

water-quallty

values

i

has

and.estimati.ng

When

available

-1[1

and

distribution

~ata.60

dlstnbutlon,

upon

regression

of

the

statistics

assumption

that

thc

distribution.60

ANAL

YSIS

OF

FLOODS

,

Lognormal,

Pearson

able

choices

18.2.

However,

select

for

a

an

suggested

in

for

individual

and

generalized

flood

distribution

a

site.

Sec.

flood

discusses

using

it

and

to

value

the

18.3.3,

section

for

and

extreme

flows

region

This

adopted

Kingdom,

3,

describing

as

procedures

United

type

for

is

to

the

sources

flow.ftequency

use

of

historical

flood

the

to

estimat~

data

United

flow

Sec.

experience

parameters

flow

in

flood

in

regional

of

of

reason.

described

use

number

analysis

the

are

methods

advisable

reduce

describes

distributions

fitting

and

States

particular

and

the,

information.

,

;,

c

{,

;ci

:i

18.7.1

Selection of Data and Sources

I\ convenient way to find information on U ni ted Stateswater data is through the U. S. ~ational Water Data Exchange (NA WDEX) assistance centers. [For information ;ontact NAWDEX, U.S. Geological Survey (USGS), 421 National Center, Reston, Va.22092; tel. 703-648-6848.] Records are also published in annual U .S. Geological ~urveywater data reports. Computerized records are stored in the National Water DataStorage and Retrieval System (W A TSTORE). Many of these records (climate jata, daily and annual maximum stream flow, water-quality parameters) have been puton compact disc read-only memories (CD-ROMs) sold by Earthlnfo Inc. (5541 CentralAve., Boulder, Colo. 80301; tel. 303-938-1788; fax 303-938-8183) so that the jata can be accesseddirectly with personal computers. The W A TSTORE peak-flow recordscontain annual maximum instantaneous flood-peak discharge and stages, anddates of occurrence as well as associated partial duration series for many sites. USGSoffices also publish sets of regression relationships (often termed state equalions) for predicting flood and low-flow quantiles at ungauged sites in the United States.

18.7.2

Bulletin 1 7B Frequency Analysis

Recommendedprocedures for flood-frequency analysesby U.S. federal agenciesare jescribedin Bulletin 17B.79Bulletin no.15, " A Uniform Technique for Determining Rood How Frequencies," released in December 1967, recommended the log-Pear)()ntype 3 distribution for use by U .S. federal agencies. The original Bulletin 17, releasedin March 1976, extended Bulletin 15 and recommended the log-Pearson type3 distribution with a regional estimator of the log-space skew coefficient. Bulletin 17A followed in 1977. Bulletin 17B was issued in September 1981 with correclionsin March 1982. ThomasJ43describes the development of these procedures. The Bulletin 17 procedures were essentially finalized in the mid-1970s, so they did notbenefit from subsequent advances in multisite regionalization techniques. Stud-iesin the 1980sdemonstrated that use of reasonable index flood procedures should providesubstantially better flood quantile estimates, with perhaps half the standard ~rror.9J.112.162 Bulletin 17 procedures are much less dependent Oh regional multisite analyses than are index flood estimators, and Bulletin 17 is firmly established in the UnitedStates,Australia, and other countries. However, App. 8 of the bulletin does describea procedure for weighting the bulletin's at-site estimator and a regional regressionestimator of the logarithms of a flood quantile by the available record lengthand the effective record length, respectively. The resulting weighted estimator reflectsa different approach to combining regional and at-site information than that employedby index flood procedures. Bulletin 17B recommends special procedures for zero flows, low outliers, historic peaks,regional information, confidence intervals, and expected probabilities for cstimatedquantiles. This section describes only major features of Bulletin 17B. The full Bulletin 17B procedure is described in that publication and is implemented in the HECWRCcomputer program discussed in Sec. 18.11. The bulletin describes procedures for computing flood flow frequency curves usingannual flood serieswith at least 10 years of data. The recommended technique fitsa Pearson type 3 distribution to the common base 10 logarithms of the peak discharges. The flood flow Q associated with cumulative probability p is then log(Qp)=X+KpS

(18.7.1)

where X and s are the sample mean and standard deviation of the base 10 logarithms, and Kp is a frequency factor which depends on the skew coefficient and selected exceedance probability; see Eq. { 18.2.28) and discussion of the log-Pearson type 3 distribution in Sec. 18.2.3. The mean, standard deviation, and skew coefficient of station data should be computed by Eq. { t 8.1.8), where Xi are the base 10 logarithms of the annual peak flows. Section 18.1.3 discussesadvantages and disadvantagesof logarithmic transformations. The following sections discuss three major features of the bulletin: generalized skew coefficients, outliers, and the conditional probability adjustment. Expected probability adjustments are also discussed. Confidence intervals for Pearson distri. butions with known and generaiized skew coefficient estimators are discussedin Sec. 18.4.3. Use ofhistoricai information is discussedin Sec. 18.7.4, mixed populations in Sec. 18.6.2, and record augmentation in Sec. 18.5.3. Generalized Skew Coefficient. Because of the variability of at-site sample skew coeffIcients in small samples, the bulletin recommends weighting the station skew coefficient with a generalized coefficient of skewness,which is a regional estimateof the log-space skewness. In the absence of detailed studies, the generalized skew coefficient G9 for sites in the U nited States can be read from Plate I in the bulletin. Assuming that the generalized skew coefficient is unbiased and independent of the station skew coefficient, the mean square error {MSE) of the weighted estimateis minimized by weighting the station and generalized skew coefficients inversely proportional to their individual mean square errors: G = Gs/MSE(Gs) + Gg/MSE(Gg) w I/MSE(Gs) + I/MSE{Gg)

{1872) ..

Here Gwis the weighted skew coefficient, Gs is the station skew coefficient, and G, is the generalized regional estimate of the skew coefficient; MSE( ) is the mean square error of the indicated variable. When generalized regional skew coefficients are read from its Plate I, Bulletin 17 recommends using MSE( Gg) = 0.302. From Monte Carlo experiments,JS9the bulletin recommends that MSE{G,) be estimated using the bulletin's Table I, or an expression equivalent to 10

a+b

;

,

MSE(Gs)

=

~

(18.7.3),t;

';ii

~q

i;;~

where

:j;(

ii

a

=

-0.33

+

0.081Gsi

if

IGsi$

0.90

!\

=

-0.52

+

0.301Gsi

if

IGsi>

=

0.94-

if

IGsi~

1.50

=

0.55

if

IGsi>

1.50

;\!

0.90

';~

";J!!

b

MSE{Gs)

Gs

is

in

Eq.

estimate

essentially

5/n

(18.7.3)

when

of

MSE(

development

G

for

small

Gsand

estimating

s)

of

0.261Gsi

(Ref.

10

..s;

MSE{Gs)

138).

to

McCuen98

skew-coefficient

and

maps,

n

s

50.

avoid

\

should

be

correlation

Tasker

and

Gg

1~

used

in

between

and

Stedingerl38

regression

place

of

Gsand

the

discuss

estimators

the

of

G

g

and

MSE(Gg).

Outliers.

from

Bulletin

the

trend

17B

of

the

defines

remaining

outliers

to

data.

"

In

be

experimental

"Data

points

which

statistics

depart

an

significanlly

outlier

:

is

often

a

;~

;'

-i" ~ ;?"

,",W 'i~ ~li " J;~

rogue observation which may result from unusual conditions or observational or recording error; such observations are often discarded. In this application low outliers are generally valid observations, but becauseBulletin 17 uses the logarithms of the observed flood peaks to fit a two-parameter distribution with a generalized skew coefficient, one or more unusual low-flow values can distort the entire fitted frequency distribution. Thus detection of such values is important and fitted distributions should be compared graphically with the data to check for problems. The thresholds used to define high and low outliers in log space are X :t KnS

(18.7.4)

whereX and S are the mean and standard deviations of the logarithms of the flood peaks,excluding outliers previously detected, and Kn is ~ critical value for sample size n. For normal data the largest observation will exceedX + KnS with a probability of only 10 percent; thus Eq. ( 18.7.4) is a one-sided.outlier test with a 10 percent significancelevel. Values of Kn are tabulated in Bulletin 17B; for 5 s n s 150, Kn can be computed by using the common base 10 logarithm of the sample size Kn = -0.9043 + 3.345 ~

-0.40461og

(n)

(18.7.5)

Flood peaks identified as low outliers are deleted from the record and a conditional probability adjustment is recommended. High outliers are retained unless historical infoffilation is identified showing that such floods are the largest in an extended period. Conditional Probability Adjustment. A conditional probability procedure is recommendedfor frequency analysis at sites whose record ofannual peaks is truncated bythe omission of peaks below a minimum recording threshold,.years with zero flow, or low outliers. The bulletin does not recommend this procedure when more than 25 percentof the record is below the truncation level. Section 18.6.3 discusses other methods. Let G(x) be the Pearson type 3 (P3) distribution fit to the r logarithms of the annual maximum floods that exceeded the truncation level and are included in the record,after deletions ofzero, low outliers, and other events. If the original record spannedn years (n > r), then an estimator of the probability the truncation level is exceededis r

qe = -( ~

18.7.6)

Aood flows exceededwith a probability q :S qein any year are obtained by solving q = qe[1 -G(x)]

(18.7.7)

to obtain

G(x)=l-!L=

I-q{~\ qe

(18.7.8) \

r J

Bulletin 17 usesEq. ( 18.7.8) to calculate the logarithms of flood flows which will beexceededwith probabilities of q = 0.50, 0.10, and 0.0 I. These three values are usedto define a new Pearson type 3 distribution for the logarithms of the flood flows whichreflects the unconditional frequency of above threshold values. The new Pearsontype 3 distribution is defined by its mean Ma, variance S~, and skew coefficient

~; ~

G Q'

which

are

calculated

as

"1~ ,':~

Ga

=

-2.50

+

3.12

log

(QO.99/Qo.90)

'::~

log

(Qo.90/QO.SO)

:;; ,;?)

Sa

=

log

(Qo.99/QO.SO) Ko.99

Ma

for

log-space

tion

skew

obtained

are

likely

to has

which

is

thus

a

A

requested

short

to

flood

flow

exceeded

cal

It

is

q.

used

to

for

3 distribu.

to

describe

near

fit

level

issue flow that

An

is

what

a

exceeded one

the

the

and

threshold low

outliers,

the

flood

question

is

estimators

should

probability

in

q

quantile what

I/T

and

which

have valid

will

be

a

18.2

be yield

\\iU

the

statisti.

il-q

thatm

XI-q:

relatively

argument

exceeded

small

suggests

with

(18.7.10)

variance

that

probability

when

both

X

two

criteria

records

and

they

estimated

q,

wants

so

lead

for

the

will

For

samples

bility

of

be

a

error.

value

However,

which

in

an

the

future

=

random

almost

(18.7.11)

variables.

If

same

design

the

estimates

q

because

of

the

one

had

effect

a

very

long

value

XI-q.

of

uncertainty

the

record,

With

shon in

the

App.

of 2

the

exceeded.

verage

II

in

size

16, on

the

17B

p

=

estimates

of

average.

the

its

the

App.

Ref.

99

percentile notes can

=

0.0

will that

be

20) =

x

provides +

forrnu.

Zps

of

the

lOOp

is

iO.99

17B II

ip

formula

for

Bulletin in

also

estimator

0.99

probability

equations

19 (see

almost-unbiased

For

exceedance

percent

Bulletin

I

be

for

used

to

~

I

+

~

)

exceeded

with

lognormal make

or

an

( 18.7.12)

a

proba.

log-PeaooD

expected

probabil-

adjustment. U

the

nfortunately

, while

expected

would

generally

the on for are

the located

estimated

the human related

expected

T-year

fitted

frequency and to

probability

probability

increase

activities

because

the

exceedance

economic

issues

square to

that

XI-q)

as to

different

that

be

distributions,

lated

to

probabilities

A

based

viewed lead

samples,

percentile

ity

are would

mean

i1-q

parameters.8.113.121

Fornormal las

il-q

or

one

P(X>

these

usinl

x,-qwhich

should

Sec.

provide =

E[i1-q]=Xl-q

equally

thc

threshold

above

zeros

hydrologist

with

wants

unresolved Most

of

be

type

quantiles G(x)

truncation

flood

agreed

estimators

Pearson

not

Fitted

distribution

the

ofestimatorsil-q.

unbiased

The

x.

the

probability

characteristics

P3

2.5.

should

Qo.so.

fundamental

estimate

records.

with

almost

for

+

( 18.7.9)

the

than

bound

and

median

if

less

Probability.

when

Eq.

the

poor

bound

Sa

-2.0

in

below

lower

-Ko.so

between

particularly

lower

-Ko.so

(Qo.so)

moments

flows

be a

Expected

be

the

offlood

values

log

coefficients

with

frequency

=

(18.7.9)"!;

economic Ba

y esian

of

bias at

in

correction computed

estimated

fixed

flood

a

damages

locations is

a

activities

in

(random)

distribution, at

fixed

eliminate event,

basin.4.121

the flood

computed expected levels.

the the

calculated a

level

whereas

can T -year

for

bias

dwellings

This

paradox by

the

damages Expected

in

corrections and ari~

hydrologist are

calcu.

probability

~

inference.121

:~~ ;~ ji~

"~~ .;i

, :

'

~'

~;"' 1\, L 1B.7.3 British Frequency Analysis Procedures !;, ~ TheFlood Studies ReportlOScontains hydrological, meteorological, and flood rout; ing studies for the British Isles. The report concluded that the GEV distribution t. providedthe best description of British and Irish annual maximum flood distribu; tionsand was recommended for general use in those countries. The three parameter PJand LP3 distribution also described the data well (Ref. 105, pp. 241,242). A key recommendation was use of an index flood procedure. The graphically derivednormalized regional flood distributions were summarized by dimensionless GEVdistributions called growth curves. Reevaluation of the original method showed thatthe L-moment index flood procedure is to be preferred.69(See Sec. 18.5.1.) The rtport distinguishes between sites with less than 10 years of record, those with 10 to 25yearsofrecord, and those with more than 25 years of record (Ref. 105, pp. 14 and 243): Siteswith n < 10 Years. The report recommends a regional growth curve with Q oblalnedfrom catchment characteristics (seeSec. 18.5.2), or at-site data extended if possibleby correlation with data at other sites (see Sec. 18.5.3). Siteswith 10 s n s 25 Years. Use either the index flood procedure with at-site data to estimate Q or, if the return period T < 2n, the Gumbel distribution. Siteswith n > 25 Years. Use either an index flood estimator with at-site data to estimateQ or, if the return period T < 2n, GEV distribution (see Sec. 18.2.2). For T > 500. Use Q with a special country-wide growth curve.

1B.7.4

Historical

Flood Information

Available at-site systematic gauged records are the traditional and most obvious sourceof infonnation on the frequency of floods, but they are of limited length. Another source of potentially valuable at-site information is historical and paleoOoodrecords. Historical information includes written and other records of large floodsleft by human observers: newspaper accounts, letters, and flood markers. The .1tffi1 paleoflood information describes the many botanical and geophysical sources of information on large floods which are not limited to the locations ofpast human observationsor recording devices.6.2s.134 Botanical data can consist of the systematic interpretation of tipped trees, scars, and abnormal tree rings along a water course providing a history of the frequency with which one or more thresholds were exceeded.77,78 Recent advances in physical pal eoflood reconstruction have focused on theuseof slack-water deposits and scour lines, as indicators of paleoflood stages,and theabsenceoflarge flows that would have left such evidence; such physical evidence or flood stagealong a water course has been used with radiocarbon and other dating techniquesto achieve a relatively accurate and complete catalog of paleofloods in favorablesettings with stable channels.86 Characterof I nformation. Different processescan generate historical and physical paJeofloodrecords. A flood leaving a high-water mark, or known to be the largest Ooodof record from written accounts, is the largest flood to have occurred in some IXriod of time which generaUy extends back beyond the date at which that flood , occurred.66 In other cases,several floods may be recorded (or none at all), because i!:ctheyexceedsome perception level defined by the location of dwellings and economic .

~iif,w 1 ,:.

~

1 check data for outliers and consistency. Interpretation of data is needed to nt for liquid precipitation versus snow equivalent and observation time differStations submitting data to NCDC are expected to operate standard equipand follow standard procedures with observations taken at standard times.169 infall frequency analysis is usually based on annual maximum series or partial on series at one site (at-site analysis) or several sites (regional analysis). Since 11data are usually published for fixed time intervals, e.g., clock hours, they yield the true maximum amounts for the indicated durations. For example, mnual maximum 24-h rainfalls for the United States are on the average 13 nt greater than annual maximum daily values corresponding to a fixed 24-h 1.63Adjustment factors are usually employed with the results of a frequency sis ofannual maximum series. Such factors depend on the number of observaI reporting times within the duration of interest. (See Ref. 172, p. 5-36). lother source of data which has been used to derive estimates of the probable mum precipitation, and to a lesser extent for rainfall frequency analysis, is the Army Corps of Engineers catalog of extreme storms. The data collection and :ssing were a joint effort of the U.S. Army Corps of Engineers and the U.S. her Bureau. Currently, a total of 563 storms, most of which occurred between and 1940, have been completed and published in Ref. 146; see also Refs. 104 123. There are problems associated with the use of this catalog for frequency 'sis. Jt may be incomplete because the criteria used for including a storm in the :Jgare not well-defined and have changed. Also, the accuracy in the estimation e storm depths varies.

~.2

Frequency

Analysis

Studies

Rainfall Frequency Atlas.63 known as TP-40, provides an extended rainfall lency study for the United States from approximately 4000 stations. The Gumistribution (Sec. 18.2.2; see also Ref. 172) was used to produce the point precipin frequency maps of durations ranging from 30 min to 24 h and exceedance labilities from 10 to 1 percent. The report also contains diagrams for making I ipitation estimates for other durations and exceedance probabilities. The U .S. .ther Bureau, in a publication called TP-49,149 published rainfall maps for duras of 2 to 10 days. Isohyetal maps (which partially supersede TP-40) for durations to 60 min are found in Ref. 46, known as HYDRO-35, and for 6 to 24 h for the em United States in NOM Atlas 2.100Examples of frequency maps can be found 'hap.3. ~or a site for which rainfall data are available, a frequency analysis can be perned. Common distributions for rainfall frequency analysis are the Gumbel, logrson type 3, and GEY distributions with K < 0, which is the standard distribution j in the British Isles.'os Ytapspresented in TP-40 and subsequent publications have been produced by rpolation and smoothing of at -site freq uency analysis results. Regional freq uency lysis, which uses data from many sites, can reduce uncertainties in estimators of "ernequantiles (Refs. 15 and 161; see Sec. 18.5.1 ). Regional analysis requires ction of reasonably homogeneous regions. Schaeferl2° found that rainfall data in shington State have CY and y which systematically vary with mean areal precipion. He used mean areal precipitation as an explanatory variable to develop a ional analysis methodology for a heterogeneous region, thereby eliminating Indary problems that would be introduced if subregions were defined. Models of daily precipitation series (as opposed to annual maxima) are con-

structed for purposes of simulating some hydrologic systems. As Chap. 3 discusses, models of daily series need to describe the persistence of wet-day and dry-day sequences. The mixed exponent distribution, and the Weibull distribution with k= 0.7 to 0.8, have been found to be good models of daily precipitation depths on rainy days, though an exponential distribution has often been used.122,173

,18.8.3

Intensity-Duration-Frequency

Curves

Rainfall intensity-duration-frequency (IDF) curves allow calculation of the average design rainfall intensity for a given exceedanceprobability over a range of durations. IDF curves are available for several U .8. cities; two are shown in Fig. 18.8.1.148When an lDF curve is not available, or a longer data base is available than in TP-25 or TP-40, a hydrologist may need to perfonu the frequency analyses necessaryto construct an IDF curve (see p. 456 in Ref. 20). IDF curves can be described mathematically to facilitate calculations. For example, one can use i=~

(18.8.1)

where i is the design rainfall intensity (inches per hour), t is the duration (minutes), c is a coefficient which depends on the exceedanceprobability, and e and fare coefficients which vary with location. 170For a given return period, the three constants can be estimated to reproduce i for three different t's spanning a range of interest. For example, for a I in 10 year event, values for Los Angeles are c = 20.3, e = 0.63, and f= 2.06, while for St. Louis c= 104.7, e= 0.89, andf=.9.44. More recently, generalized intensity-duration-frequency relationships for the United States have been constructed by Chenl9 using three depths: the 10-year I-h rainfall (RIO),the 10-year 24-h rainfall (R!~), and the 100-year I-h rainfall (RlOO)from TP-40. These depths describe the geographic pattern of rainfall in tenus of the depth-duration ratio (RT / Rr4) for any return period T, and the depth-frequency ratio (RlOOjRlO)for any duration t. Chen's general rainfall IDF relation for the rainfall depth RT in inches for any duration t (in minutes) and any return period T (in years) IS alRJO [xRT=

1)log(TpjIO) + I] (-k) (t+bl)Cl

(18.8.2)

where x = (R loo/ R 10),Tp is the return period for the partial duration series(equal to the reciprocal of the average number of exceedancesper year), and al , hl , and c1are coefficients obtained from Fig. 18.8.2 as functions of(RIO/R!~) with the assumption that this ratio does not vary significantly with T. Chen uses Tp = -I jlri ( I -I j 7) to relate Tp to the return period T for the annual maximum senes (see Sec. 18.6.1); for T > 10 there is little difference between the two return periods. The coefficients obtained from Fig. 18.8.2 are intended for use with TP-40 rainfall quantiles. For many design problems, the time distribution of precipitation (hyetograph) is needed. In the design of a drainage system the time of occurrence of the maximum rainfall intensity in relation to the beginnihg of the stonu may be important. Design hyetographs can be developed from IDF curves following available procedures (see Chap.3).

ployed to estimate the joint probability distribution of selected storm characteristics (such as maximum storm-center depth, storm shape parameters, storm orientation, storm depth spatial variability, etc. ). Then, the probability distribution of the position of the storm centers within the transposition area is determined. This distribution is generally not uniform because the likelihood of storms of a given magnitude will vary across a region.45 Integration over the distribution of storm characteristics and storm locations allows calculation of annual exceedance probabilities for various catchment depths. An advantage of the SST method is that it explicitly considers the morphology of the storms including the spatial distribution of storm depth and its relation to the size and shape of the catchment of interest-

,'. ,i;" 'r :i c

t 8.9

FREQUENCY

ANAL

YSJS

OF

LOW

FLOWS

Low-flow quantiles are used in water-quality management applications including waste-Ioad allocations and discharge permits, and in siting treatment plants and sanitary landfills. Low-flow statistics are also used in water-supply planning to deter\~ mine allowable water transfers and withdrawals. Other applications of low-fl.ow tc frequency analysis include detennination of ~inimum downstream r~~e~se requlre!;{~!~ments from hydropower, water-supply, coolIng plants, and other facultIes. ~" ~ !~;i ~"; "I:' '; ~,

Annual- Event- Based Low- Flow Statistics. Sources of stream flow and low-flow data are discussed in Sec. 18.7. 1. The most widely used low-flow index in the U nited States is the one in 10-year 7-day-average low flow, denoted Q7.0.l0 (Ref. 117). In general, Qd,p is the annual minimum d-day consecutive average discharge not ex-

i: ~); ~ ~.~ ~{; ~p; ~) ld~\ i~ ~r:

ceeded with probability p. Prior to performing low-flow frequency analyses, an effort should be made to "deregulate" low flow series to obtain "natural" streamflows. This includes accounting for the impact of large withdrawals and diversions including water- and wastewater-treatment facilities, as well as urbanization, lake regulation, and other factors. Since low flows result primarily from groundwater inflow to the stream channel, substantial year-to-year carryover in groundwater storage can cause seq~enc~s of annual mi~imu.m low flows to be correlat~d from one year ~o the next (see FIg, 9 ~n Ref. 11 ? or ~Ig. 2 In Ref. 157). Low-flow ~enes should be subjected to trend analysIs so that IdentIfied trends can be reflected In frequency analyses.

18.9.1

Selection

of Data

and

Sources

~;r[ ~f:: The Flow-Duration Curve. Flow-duration curves are an alternative to analysis of fb:;: annual minimum d-day averages (see Sec. 8.6.1 ). A flow-duration curve is the empir~!; ical cumulative distribution function of all the daily ( or weekly) stream flow recorded ~1;}cat a site, A flow-duration curve describes the fraction of the time over the entire ~;:ffi record that different daily flow levels were exceeded. Flow-duration curves are often ~~t u~d in hydrol?gic studies for run-o!-river hydropower, water supply, irrigation pl~nW;iJ mng and design, and water-quallty management.J6,40.4J,5J.l21 The flow-duratIon -;'~;- curves should not be interpreted 1 on an annual event basis, as is Qd,p, because the .);¥{ Dow-durat~on curve pro.vi?es o.nly the fraction of the time that a ~tr~am.flow level was 1~~1exceeded; It does not dIstInguish between regular seasonal vanatlon In flow levels, """ ~; and random variations from s~asonal averages. 1 1~"'c .~)';

~I;

18.9.2 Zeros

Frequency Analysis Methods for Low Flows and Treatment

of

Estimation of Qd,pfrom stream flow records is generally done by fitting a probability distribution to the annual minimum d-day-averagelow-flow series.The literature on low-flow frequency analysis is relatively sparse.The extreme value type III or Weibull distribution (see Sec. 18.2.2) is a theoretically plausible distribution for low flows. Studies in Canada and the eastern United States have recommended the three. parameter Weibull, the log-Pearson type 3, and the two-parameter and three. parameter lognormal distributions basedon apparent goodness-of-fit.24.IJ9.ls4 Fitting methods for complete samples are described in Sec. 18.2. Low-flow series often contain years with zero values. In some arid areas,zero flows are recorded more often than nonzero flows. Stream flows recorded as zero imply either that the stream was completely dry or that the actual stream flow was below a recording limit. At most U.S. Geological Survey gauges there is a lower stream flow level (0.05 ftJ/s) below which any measurement is reported asa zero.This implies that low-flow series are censored data sets, discussed in Sec. 18.6.3. Zero values should not simply be ignored, nor do they necessarily reflect accurate measurements of the minimum flow in a channel. Based on the hydraulic configuration of a gauge, and knowledge of the rating curve and recording policies, one can gener. ally determine the lowest discharge which can be reliably estimated and would not be recorded as a zero. The plotting-position method or the conditional probability model in Sec. 18.6.3 are reasonable procedures for fitting a probability distribution with data setscontain. ing recorded zeros. The graphical plotting position approach without a formal statistical model is often sufficient for low-flow frequency analyses. One can define visu. ally the low-flow frequency curve or estimate the parameters of a parametric distribution using probability-plot regression. .

18.9.3

Regional Estimates of Low-Flow

Statistics

Regional regression procedures are often employed at ungauged sites to estimate low-flow statistics by using basin characteristics. If no reliable regional regression equations are available, one can also consider the drainage area ratio, regional statistics, or base-flow correlation methods described below. Regional Regression Procedures. Many investigators have developed regionaJ models for the estimation oflow-flow statistics at ungauged sites using physiographic basin parameters. This methodology is discussedin Sec. 18.5.2. Unfortunately, most low-f1ow regression models have large prediction errors, as shown in Fig. 18.5.1, becausethey are unable to capture important land-surface and subsurface geologicaJ characteristics of a basin. In a few regions, efforts to regionalize low-flow statistics have been improved by including basin parameters which in some manner describe the geohydrologic response of each watershed.9.18.ls6 Conceptual watershed rnodeIs can be used to formulate regional regression models oflow-f1ow statistics.ISS Drainage Area Ratio Method. Perhaps the simplest regional approach for estima. tion of low-flow statistics is the drainage area ratio method, which would estimatea low-flow quantile Ypfor an ungauged site as Yp = ( ~ \ Xp \A x J

(18.9.1) , " ;;~

;,:1 .

whereXpis the corresponding low-flow quantile for a nearby gauging station and A andAy are the drai~a~e areasfor ~hegau~ng station and ungauged site, respectively: Seepage~ns, consistmg of a senes of discharge measurements along a river reach ~unng pe~ods ofbase ~ow, are usefu~for deteffilining the applicability of this simple lInear dramage-area discharge relatlOn.116Some studies employ a scaling factor {Ay/Ax)bto allow for lossesby using an exponent b < 1derived by regional regression. (SeeSec. 18.5.2.) Regional Statistics Methods. One can sometimes use a gauging station record to construct a monthly streamflow record at an ungauged site using y(i,}) = M(y;) + SO'i) [x(~t2;~ M(x;~l

(18.9.2)

wherey(i,)) and x(i,)) are monthly stream flows at the ungauged and nearby gauged sites,respectively, in month i and yearj; M(x;) and S(x;)are the mean and standard deviation of the observed flows at the gauged site in month i; and M(yi) and S(Y.) are the corresponding mean and standard deviation of the monthly flows at th~ ungaugedsite obtained from regional regression equations, discussed in Sec. 18.5.2. Hirsch64found that this method transfe1Tedthe characteristics of low flows from the gaugedsite to the ungauged site. BaseFlow Correlation Procedures. When baseflow measurement.\'(instantaneous or averagedaily values) can be obtained at an otherwise ungauged site, they can be correlatedwith concurrent stream flows at a nearby gauged site for which a long flow recordis available.114.116 Estimators oflow-flow moments at the ungauged site can be developedby using bivariate and multivariate regression, as well as estimates of their standard errors.l44 This is an extension to the record augmentation idea in Sec. 18.5.3.Ideally the nearby gauged site is hydrologically similar in teffils of topography, geology, and base flow recession characteristics. For a single gauged record, if regressionof concurrent daily flows at the two stations yields a model y=a+bx+E

with Var(E)=S~

(18.9.3)

estimatorsof the mean and variance of annual minimum d-day-average flows yare M(y) = a + b M(x) S2(y) = b2 S2(X) + s~

(18.9.4)

where M(x) and S2(X) are the estimators of the mean and variance of the annual minimum d-day averagesat the gauged x site. Base flow correlation procedures are subjectto considerable error when only a few discharge measurements are used to estimatethe parameters of the model in Eq. ( 18.9.3), as well as error introduced from useof a model constructed between base flows for use in relating annual minimum d-day averages.(Thus the model R2 should be at least 70 percent; seealso Refs. 56 and 144.)

18. 10 FREQUENCY ANAL YSIS OF WA TERQUALITY v ARIABLES

,

In the early 1980s,most U .S.water-quality improvement programs aimed at obvious and visible pollution sources resulting from direct point di.\'charges of sewageand

wastewatersto surface waters. This did not always lead to major improvements in the quality of receiving waters subject to non-point-source pollution loadings, conesponding to storm-water runoff and other discharges that carry sediment and pollutants from various distributed sources. The analyses of point- and nonpoint-source water-quality problems differ. Point sources are often sampled regularly, are less variable, and frequently have severe impacts during periods of low stream flow. Nonpoint sources often occur only during runoff-producing storm events, which is .when nonpoint discharges generally have their most severe effect on water quality.

18.10.1

Selection of Data and Water-Quality

Monitoring

Water -quality problems can be q uantified by a num ber of variables, including biodegradable oxygen demand (BOD) and concentrations of nitrogenous compounds, chlorophyll, metals, organic pesticides, and suspended or dissolved solids. Waterquality monitoring activities include both design and actual data acquisition, and the data's utilization. 102. 167Monitoring programs require careful attention to the definition of the population to be sampled and the sampling plan.47A common issueis the detection of trends or changesthat have occurred over time becauseof development or pollution control efforts, as discussed in Chap. 17. The statistical analysis of water-quality data is complicated by the facts that quality and quantity measurements are not always made simultaneously, the time interval between water-quality samples can be irregular, the precision with which different constituents can be measured varies, and base flow samples (from which background levels may be derived) are often unavailable in studies of non-pointsource pollution. The U .S. National Water Data Exchange provides accessto dataon ambient water quality; seeSec. 18.7.1. A list of other non-point-source water-quality data basesappears in Ref. 76. .

18.10.2

Frequency Analysis Methods and Water-Quality

Data

It is often inappropriate to use the conventional approach of selecting a single design flow for managing the quality of receiving waters. The most critical impact of pollutant loadings on receiving water quality does not necessarily occur under low flow conditions; often the shock loads associated with intennittent urban storm-water runoffare more critical.99 Nevertheless, water-quality standards are usually statedin tenns ofa maximum allowable d-day average concentration. 147The most common type of design event for the protection of aquatic life is based on the one in T-year d-day average annual low stream flOW.lO For problems with regular data collection programs yielding continuous or regularly spaced observations, traditional methods of frequency analysis can be employed to estimate exceedance probabilities for annual maxima, or the probability that monthly observations will exceed various values. Event mean concentrations (EMC) corresponding to highway storm-water runoff, combined sewer overflows, urban runoff, sewagetreatment plants, and agricultural runoff are often well approximated by lognormal distributions,39 which have been a common choice for constituent concentrations.47 Section 18.1.3 discusses advantages and disadvantagesof logarithmic transformations. Procedures in Sec. 18.3 for sel~cting an appropriate probability distribution may be employed. Investigations of the concentrations associatedwith trace substances in receiving waters are faced with a recurring problem: a substantial portion of water sample

concentrations are below the limits of detection for analytical laboratories. Measure~ c

ments

below

the

detection

limit

are

often

reported

as

"less

than

the

detection

limit"

i rather than by numerical values, or as zero.47 Such data sets are called censored data '( in the field of statistics. Probability-plot regression and maximum likelihood techIi ;,!cccniques for parameter estimation with censored data sets are discussed in Sec. 18...6 3 60,61 ), I '"" ~t~ ~~:; ,

",,'C"

i i~[

For

intermittent

roughly

to

partial

storm-water

~*Mii been

loading

problems

duration

series,

problems,

estimated

the

the

situation

discussed

average

is more

in

Sec.

recurrence

difficult

18.6.1.

interval

In

(in

and

the

years)

corresponds

context

of

of a design

urban

event

has

as

I ~Ii' i:.'

I T=

NP(C;2:

Co)

(18.10.1)

, where

N is the

.,'

concentration,

f~"

dard

Z'% '""" ~'llic

concentrations ( 18.6.3) of

~;\~!.

Chap.

~"; j:,'"" , t~' ;;;" f;~:;,

Co

18.

in

a

number

and

probability

the

runoff

can

event

of

is

rainfall

runoff

P( C ;2: Co) obtained

by

measured in observed 18.6.1 with event arrival

Sec.

21)

11

average

also

be

used

to

COMPUTER

FREQUENCY

estimate

an

fitting

runoff rate

in

a year,

C is a constituent

observation

C exceeds

a frequency

a stan-

distribution

events.37 This A. = N. Models

to

corresponds such as SWMM

the

to

Eq. (see

T directly.

PROGRAMS

ANAL

events that

FOR

--!

YSIS

!~;i ~i~c:

Many

~;j,~

standard

frequency

computations

~]t

packages.

(~'c

can be quite

involved.

1;;1

packages

perform

f,:,

flood

functions

on

However,

to

frequency

are

hand

maximum

the

simple

and

spreadsheets,

or

likelihood

Water

estimators

management

standard

analyses

relatively

calculators,

are

and

agencies

procedures

discussed

are

performed

several

in most

they

easily

geFleral-purpose other

countries

employ.

Four

with

statistical procedures

have sets

of

computer routines

for

below.

~::; ~~,

U.S.

Army

Corps

of

~:;i,

U.S.

Army

Corps

of Engineers

~;,"

support

~~ir

routines

statistical

~~;

purpose

~r~:

positions,

!~:',~ ~c" ~-;

Frequency drologic

V;;,t

and

t]~i

support,

for

and

vey

Geological

r1;itf ~j~

British

~~~i

tation

graphical

HEC

als?

display.

Center,

Flood

Studies

ogy.l0S

It also

graph,

and

I

The

package Hydrology,

NatIonal

contains

reservoir and

personal

by

several

for

training

St.,

Stop

COE

probable for can

user

vendors.

The

analyses

(Chief

Reston,

developed maximum personal

be obtai~ed

Crowmarsh

Glfford,

HECWRC and

user Sur-

Hydrologist,

U.S.

22092).

implemenInstitute

precipitation,

from

Hy-

U .S. Geological

the

computers

Hood

interface

Va.

by

plotting

Center,

95616-46897. the

general-

contacting

Support

to

includes as

curves,

by 51

Calif.

437,

well

duration

obtained

in

library

as

is a microcomputer-based

methods

calculations inf.n:nation

Davis,

I!B

Mall

Micro-FSR

Bulldmg,

Army,

private

Bullet!n

Center,

menu-driven

Maclean

the

series, be

routines

The

analyses,

improvements

analysis

routing

time

The

FORTRAN

computers.

can

of

(HECWRC).

of60

17B

statistics,

Second

Analysis

a library

Bulletin

some

Software.

flood-frequency

~~'

609

a pr~gram

Frequency

Information

with

marketed

Survey,

of the

MS-DOS standard

general

software,

provides

Flow

developed

Department

actively

~~:

~~ i~;;:

on

including

Engineering

are

Flood has

the

Investigations,

other

t!\~'

of

analysis

performing

functions

'I

~~;:

Engineers

ofHydrolunit

running

So~ware Walhngford,

hydro-

MS-DOS.

Sales,

Instit~te

Oxfordshlfe

OXIO 38800, United fax, 0491 32256

Kingdom;

phone, 0491 38800; telex, 849365 HYDROL

a;

Consolidated Frequency Analysis (CFA) Package. The package incorporates FORTRAN routines developed by Environment Canada for flood-frequency analyses in that country with MS-DOS computers. Routines allow fitting of three-parameter lognormal, Pearson, log.;.Pearson, GEV, and Wakeby distributions using moments, maximum likelihood, and sometimes probability-weighted moment analyses. Capabilities are provided for nonparametric trend and tests of independence, as well as employing maximum likelihood procedures for historical information with a single threshold. Contact Dr. Paul Pilon, lruand Waters Directorate, Water Resources Branch, Ottawa, Ontario KIA OE7, Canada. FORTRAN

Routines for Use with the Method ofL Moments.

Hosking 14describesa

set of FOR TRAN- 77 subroutines useful for analyses employing L moments, including subroutines to fit 10 different distributions. Index-flood procedures and regional diagnostic analyses are included. Contact Dr. J. R. M. Hosking, Mathematical Sciences Dept., IBM Research Division, T. J. Watson Research Center, Yorktown Heights, N. Y. 10598. The routines are available through ST A TLIB, a system for distribution of stati:)tical software by electronic mail. To obtain the software send the message "send Imoments from general" to the e-mail address: [email protected] u.ed u.

REFERENCES I. Abramowitz, M., and I. A. Stegun, Handbook of Ma[hema[ical Func[ions, Appl. Math. Ser. 55, National Bureau of Standards, U.S. Government Printing Office, Washington, D.C., 1964. 2. Adamowski, K., " A Monte Carlo Comparison ofParametric and Nonparametric Estimation ofFIood Frequencies," J. Hydrol., vol. 108, pp. 295-308, 1989. 3. Ang, A. H-S., and W. H. Tang, Probabili[y Conceptsin Engineering, Planning and Design, vol. I, Basic Principles, Wiley, New York, 1984. 4. Ameli, N., "Expected Annual Damages and Uncertainties in Flood Frequency Estimation," J. Water Resour. Plann. Manage., vol. 115, no.1, pp. 94-107, 1989. 5. Arora, K., and V. P. Singh, " A Comparative Evaluation of the Estimators of the Log Pearson Type (LP) 3 Distribution," J. Hydrol., vol. 105, no. 1/2, pp. 19-37,1989. 6. Baker, V. R., R. C. Kochel, and P. C. Patton, eds., Flood Geomorphology, Wiley, New York, 1988. 7. Bardsley, W. E., "Using Histori
Xem th�m...

Komen

Hak cipta © 2017 SLIDEMY.COM Inc.