In ilapros/ilaprosUtils: A Miscellaneous Collection of Functions I Tend to Use

knitr::opts_chunk$set(echo = TRUE, fig.retina=2, fig.width=7, fig.height=5)

This document presents how to use R to combine historical flow data and systematic gauged records to estimate flood frequency curves. The functions discussed in the document are mostly contained in the ilaprosUtils package developed by Ilaria Prosdocimi. Functions to include historical data under a Bayesian inference approach are included in the external and independently developed package nsRFA, which is not discussed further in here. The reader interested in the Bayesian formulation of the estimation procedure is encouraged to read the help page of the package ?nsRFA::BayesianMCMC function and the vignette Local and regional analyses with BayesianMCMC. When using the nsRFA::BayesianMCMC function there might be sometimes issues in the convergence of the MCMC chains, something that can not easily be fixed within the nsRFA sets of functions.

To illustrate the different methods it is assumed that some 40 years of gauged flow are available at at a certain gauging station.

par(mar = c(2, 2, 2, 0.5))
syst <- data.frame(Year = seq(1976,2015), 
                   Flow = c(300.346, 527.381, 410.964, 224.615, 370.688, 275.477, 
                            345.358, 439.179, 117.119, 371.371, 386.094, 458.024, 
                            286.608, 264.741, 262.419, 453.842, 715.718, 480.023, 
                            316.092, 479.976, 457.083, 317.254, 378.389, 249.409,
                            384.964, 228.965, 291.549, 370.726, 235.923, 229.726, 
                            376.540, 326.176, 553.550, 318.610, 105.712, 309.667, 
                            375.903, 523.235, 257.639, 229.343))
with(syst, plot(Year,Flow,type="h"))

Assuming that the data come from a GLO distribution (the default distribution for British flow series), estimates for the distribution parameter based on the at-site series can be obtained wither via the L-Moment or the maximum likelihood approach as shown below.

lmom::samlmu(syst$Flow) ## at-site L-moments
## GLO parameter estimates via L-moment approach
systOnlyLmom <- lmom::pelglo(lmom::samlmu(syst$Flow)) 
## GLO parameter estimates via Maximum Likelihood approach
systOnlyMaxL <- ilaprosUtils::glod.fit(syst$Flow, show = FALSE)
systOnlyMaxL$conv ## if 0 model has converged
systOnlyLmom; systOnlyMaxL$mle ## not very different

The resulting flood frequency curves are shown below (notice the use of the retPlot function from the ilaprosUtils package).

## ?ilaprosUtils::retPlot
rlSOML <- ilaprosUtils::retPlot(systOnlyMaxL, pch = 16, sign.alpha = 0)
lines(log(rlSOML$p/(1-rlSOML$p)), 
      lmom::quaglo(f = rlSOML$p, para = systOnlyLmom), col = 2)

Imagine that some historical information was available for the area of interest and we knew that three floods had happened in the 250 years before the beginning of the systematic record, namely in year 1770, 1850 and 1920, with peak flow values assessed to be of size 840, 820 and 870 $m^3/s$ respectively. From the careful analysis of historical records it is deemed that these three flow records are the biggest and only events above 800 $m^3/s$ which have happened in the historical period.

The historical data are combined with the systematic record in the Figure below.

hist <- data.frame(Year = c(1770, 1850, 1920),  ## only needed for plotting
                   Flow = c(840, 820, 870))
with(syst, plot(Year, Flow, type="h", xlim =  c(1976-250, 2015), ylim = c(0,900)))
lines(hist$Year, hist$Flow, type="h", col = 2)
abline(h = 800, col = 4, lty = 4) ## perception threshold

To include the historical data in the estimation procedures the typical L-moments and maximum likelihood approaches need to be slightly modified as outlined in Wang (1990) and Stedinger and Cohn (1986). L-moments based on the Partial Probability Weighted Moment approach can be obtained using the ilaprosUtils::lmom.hist.fit function, while the extension for the maximum likelihood can be applied using ilaprosUtils::glo.hist.fit (at time of writing only the GLO and GEV distribution were implemented for the maximum likelihood approach - see ilaprosUtils::gev.hist.fit).

ll <- ilaprosUtils::lmom.hist.fit(
            xdat = c(hist$Flow,syst$Flow), ## all data combined
            h = 250, ## historical period length
            X0 = 800, ## perception threshold
            k = 3, ## number of historical peaks
            nmom = 4) ## number of L-mom parameters 

## GLO parameter estimates via L-moment approach
withHistLmom <- lmom::pelglo(ll); rm(ll)
## GLO parameter estimates via Maximum Likelihood approach
withHistMaxL <- ilaprosUtils::glo.hist.fit(
            xdat = c(hist$Flow,syst$Flow), ## all data combined
            h = 250, ## historical period length
            X0 = 800, ## perception threshold
            k = 3, ## number of historical peaks
            show = FALSE)
withHistMaxL$conv ## if 0 model has converged - check!
withHistLmom; withHistMaxL$mle ## somewhat different

The resulting flood frequency curves are shown in the Figure below, with the dashed lines indicating the estimates obtained when including the historical information in the estimation procedure.

## ?ilaprosUtils::retPlot
## retPlot recognises that some data are from an historical record
## displayed as filled squares
rlWHML <- ilaprosUtils::retPlot(withHistMaxL, 
                p = c(seq(0.01,0.92,l=50),seq(0.92,.998,l=20)),
                pch = 16, sign.alpha = 0, lty = 2)
## Sysy. only  max lik
lines(log(rlWHML$p/(1-rlWHML$p)), 
      lmom::quaglo(f = rlWHML$p, para = systOnlyMaxL$mle), col = 1)
## Sysyonly - Lmoments
lines(log(rlWHML$p/(1-rlWHML$p)), 
      lmom::quaglo(f = rlWHML$p, para = systOnlyLmom), col = 2)
## With Historical  info - Lmoments
lines(log(rlWHML$p/(1-rlWHML$p)), 
      lmom::quaglo(f = rlWHML$p, para = withHistLmom), col = 2, lty = 2)

legend("topleft", col = c(1,1,2,2), lty = c(1,2,1,2), bty = "n",
       legend = c("Syst. only - ML", "With Hist. - ML",
                  "Syst. only - Lmom", "With Hist. - Lmom"))
legend("bottomright", pch = c(16,15), bty = "n",
       legend = c("Systematic data", "Historical data"))

In some cases it might happen that the actual values of the historical peak flows can not be assessed with sufficient confidence and only the information on the threshold exceedance is known. That is to say it is only known that a flood larger than X0 has occurred but the size of the peak cannot be assessed: this case is referred in Stedinger and Cohn (1986) as binomial censoring. The Partial Probability Weighted Moments approach was not extended to cope with this binomial censored data, while the maximum likelihood approach can easily include such information, as shown below.

## GLO parameter estimates via Maximum Likelihood approach
## Under Binomial censoring
withBinCMaxL <- ilaprosUtils::glo.hist.fit(
            ### The first points only need to be larger than X0           
            xdat = c(801, 801, 801 ,syst$Flow), ## all data combined
            h = 250, ## historical period length
            X0 = 800, ## perception threshold
            k = 3, ## number of historical peaks
            binomialcens = TRUE,
            show = FALSE)
withBinCMaxL$conv ## if 0 model has converged - check!
withBinCMaxL$mle; withHistMaxL$mle ## little difference

\

References

Stedinger, J. R. and Cohn, T. A. (1986). Flood Frequency Analysis With Historical and Paleoflood Information. Water Resources Research, 22, 785--793. \

Wang, Q.J.(1990). Unbiased estimation of probability weighted moments and partial probability weighted moments from systematic and historical flood information and their application. Journal of hydrology, 120, 115--124.

ilapros/ilaprosUtils documentation built on April 6, 2023, 4:44 a.m.