Envelope curve for record Australian Rainfall

I couldn’t find an envelope curve for Australian record rainfall so made one as shown below.

record rain-au

For anyone who is interested, the details are as follows.

Record rainfalls for Australia for each state and territory for a range of durations are summarised on the Bureau of Meteorology web site (link).

We can use the readHTMLTable function from the XML package to scape the data off the web.  The following code does this and then combines the data for all jurisdictions into a single dataframe.

# Read in the tables for each jurisdiction
au.rain <- readHTMLTable(my.url,
header = FALSE,
colClasses = c('character', 'character', 'numeric', 'numeric', 'character'),
encoding = "UTF-8", )

# set header = FALSE so column names are ignored, there are small inconsistencies
# in column names which was causing trouble
# set encoding to "UTF-8" which seems to be what the Bureau is using

# name the juridication for each data frame
names(au.rain) <- c('ACT', 'NSW', 'NT', 'QLD', 'SA', 'Tas', 'Vic', 'WA')
# bind into a single data frame
au.rain <- do.call('rbind', au.rain)

# add row names as a column (so we can find where values come from)

au.rain <- au.rain %>% mutate(source.row = row.names(au.rain))

# name columns

names(au.rain)[1:5] = c('location', 'date', 'duration.min', 'rainfall.mm', 'source')

Next steps are to:

  • Get the maximum rainfall for each duration (and store this in the dataframe au.rain.max)
  • Calculate duration in hours
  • Remove any records where there is already a larger depth for a shorter duration, this will leave behind the extreme values.
# get maximum rainfall for each duration

au.rain.max <- au.rain %>%
group_by(duration.min) %>%
slice(which.max(rainfall.mm))

# Calculate duration in hours and add to data frame

au.rain.max <- au.rain.max %>%
mutate(duration.hour = as.vector(duration.min/60))

# create an increasing series
# start with the first value and only accept the next value if its larger

au.rain.max <- au.rain.max %>%
ungroup() %>%
mutate(running.max = cummax(rainfall.mm)) %>%
filter(rainfall.mm >= running.max)

The remaining record rainfalls are in the table below and are available as a csv file here.

Location Date Duration (min) Rainfall (mm)
Adelaide, SA 3 May 1942 2 11
Warawarralong, NSW 25 Jan 1904 3 27
Toorooka, NSW 8 Mar 1964 6 51
Wandoona, NSW 17 Feb 1983 10 86
Nobby, QLD 11 Mar 1939 20 203
Florence, QLD 4 Dec 1920 60 230
Gunyerwarilli, QLD 20 Nov 1937 80 270
Upper Ross, QLD 3 Mar 1946 120 305
Kumbia, QLD 14 Mar 1941 180 356
Mt Kumba, QLD 13 Mar 1941 210 406
Binbee, QLD 6 Jan 1980 270 607
Wongawilli, NSW 18 Feb 1984 540 653
Wongawilli, NSW 18 Feb 1984 720 717
Duck Creek, QLD 20 Jan 1970 840 864
Bellenden Ker Top, QLD 4 Jan 1979 1440 960
Bellenden Ker Top, QLD 4 Jan 1979 1800 1330
Bellenden Ker Top, QLD 5 Jan 1979 2520 1577
Bellenden Ker Top, QLD 5 Jan 1979 2880 1947
Bellenden Ker Top, QLD 5 Jan 1979 4320 2517
Bellenden Ker Top, QLD 8 Jan 1979 11520 3847

 

The next step is to calculate the envelope curve – a straight line on a log-log plot of rainfall against duration that provides an upper bound of the record rainfall depths.  The way this has been done by others (e.g. Castellarin et al., 2009)  is to regress rainfall against duration, to obtain the slope of the curve, and adjust the intercept so that the envelope is greater than or equal to all recorded rainfalls.  Here I use robust regression to obtain the slope estimate which limits the influence of outlying data values.

m.1 <- MASS::rlm(log(rainfall.mm) ~ log(duration.hour), data = au.rain.max)
coef(m.1)
# (Intercept) log(duration.hour)
# 5.2959970 0.5696566

# So the slope we need is 0.5696566
# The line needs to pass through the 5th data point which determines the intercept
au.rain.max[5, ]

# location date duration.min rainfall.mm source source.row duration.hr duration.hour running.max
# 1 Nobby 11 Mar 1939 20 203 QLD.10 0.3333333 0.3333333 203

env.intercept <- as.vector(203/((1/3)^coef(m.1)[2]))
env.intercept
# [1] 379.5695

# Envelope curve
#380*d^0.57

So our first cut at an envelope curve for Australian rainfall is:

R = 380 D^{0.57}

Where, R is rainfall depth in mm and D is duration in hours.

Record rainfalls and envelope curves are also available the world (WMO, 2009) and for New Zealand (Griffiths et al., 2014).  These are shown on the figure below.  The envelope curve I proposed for Australia (the green line) looks much too steep as it crosses the world curve.  The New Zealand curve also looks a little steep.

Rain-comb

Record rainfalls and envelope curves for the world, Australia and New Zealand

The plot was made as follows:

#remove(list = objects())

# World Rainfall Records and Envelope curve
# download CSV file extracted from Annex II of WMO (2009)

# WMO (2009) Manual on estimation of probable maximum precipitation (PMP). WMO-No. 1045. 
# World Meteorological Organization. http://www.wmo.int/pages/prog/hwrp/publications/PMP/WMO%201045%20en.pdf


library(repmis)
library(dplyr)

my.url <- 'https://dl.dropboxusercontent.com/u/10963448/RecordRainfall-World.csv'

world.rain <- tbl_df(repmis::source_data(my.url, skip = 3, sep = ','))

world.rain <- world.rain %>%
 mutate(duration.hour = `Duration (min)`/60) # work in hours


# World Envelope Curve
# WMO (2009) 

# R = 491 D^0.452

env.world <- function(d) {
 491 * d^0.452
}

# Range of durations

# from 0.167 hours to 120 hours
d.seq.world <- seq(from = min(world.rain$duration.hour), to = max(world.rain$duration.hour), length.out = 100) 



# New Zealand rainfall records
# from Griffiths, G. A., A. I. McKerchar and C. P. Pearson (2014). 
#"Towards prediction of extreme rainfalls in New Zealand." 
# Journal of Hydrology (NZ) 53(1): 41-52.

nz.rain <- data_frame(
 rainfall = c(34, 134, 428, 566, 758, 1049, 2927, 18405),
 duration.hour = c(0.167, 1, 8.5, 12, 24, 48, 730, 8760)
)


# Envelope curve for NZ
# Griffiths envelope curve
# from Griffiths, G. A., A. I. McKerchar and C. P. Pearson (2014)

env.griffiths <- function(d) {
 145 * d^0.55
}

# from 0.167 hours to 8760 hours (1 years)
d.seq.griffiths <- seq(from = 0.167, to = 8760, length.out = 100) 

# Australian data and envelope curve 

# loading these in from dropbox

my.url <- 'https://dl.dropboxusercontent.com/u/10963448/RecordRainfall-au.csv'
au.rain.max <- tbl_df(repmis::source_data(my.url, skip = 3, sep = ','))


au.rain.max <- au.rain.max %>%
 mutate(duration.hour = `Duration (min)`/60) # work in hours


# Australian Envelope Curve

# R = 380 D^0.57

env.au <- function(d) {
 380 * d^0.57
}

# from 0.033 hours to 192 hours 
d.seq.au <- seq(from = 0.033, to = 192, length.out = 100) 


# Plot
par(oma = c(0,3,0,0))

plot(`Depth (mm)` ~ duration.hour, data = world.rain, 
 log = 'xy',
 xlab = 'Duration (hr)',
 ylab = '',
 xaxt = 'n',
 yaxt = 'n',
 pch = 21,
 bg = 'red')

my.labels = format(axTicks(1), big.mark = ',', scientific = FALSE, digits = 0, trim = TRUE)

axis(side = 1, at = axTicks(1), labels = my.labels)
axis(side = 2, at = axTicks(2), labels = format(axTicks(2), big.mark = ','), las = 2)

mtext(side = 2, text = 'Rainfall (mm)', line = 4 )


points(rainfall ~ duration.hour, data = nz.rain, pch = 21, bg = 'grey', col = 'black') # add NZ records

points(`Depth (mm)` ~ duration.hour, data = au.rain, pch = 21, bg = 'green', col = 'black') # add Australian records


lines(d.seq.griffiths, env.griffiths(d.seq.griffiths), lty = 3, col = 'grey', lwd = 2) # Griffiths envelope curve (NZ)
lines(d.seq.world, env.world(d.seq.world), lty = 2, col = 'red', lwd = 2) # World envelope curve
lines(d.seq.au, env.au(d.seq.au), lty = 2, col = 'green', lwd = 2) # Australian envelope curve

legend('topleft', 
 lty = 2, 
 pch = 21,
 col = c('red', 'green', 'grey'), 
 legend = c('World', 'Australia', 'New Zealand'), 
 pt.bg = c('red', 'green','grey'),
 bty = 'n'
)

text(0.5, 1000, col = 'red', lwd = 2, expression(paste('R = 491 ',D^0.452))) # World
text(50, 7000, col = 'dark green', lwd = 2, expression(paste('R = 380 ',D^0.57))) # Australia
text(200, 500, col = grey(0.2), lwd = 2, expression(paste('R = 145 ',D^0.55))) # New Zealand



Clearly, the envelope curve for Australia needs to be flatter.  It seems reasonable to adopt the slope of the world curve for the Australian envelope and adjust the intercept to pass through the highest rainfalls.

The resulting curve is:

R = 364.2 D^{0.452}

When plotted, this looks Ok.

Calculations are as follows:

# Review Australian envelope curve
# Use the slope of the world envelope curve
# The critical point that the envelope curve passes through is from row 19 of the Australian data
#
au.rain.max[19, ]
# location date duration.min rainfall.mm source source.row duration.hour running.max
# 1 Bellenden Ker Top 5 Jan 1979 4320 2517 BoM QLD.41 72 2517

env.intercept <- as.vector(2517/(72^0.452))
env.intercept
# [1] 364.2243

Now, plot the data, envelope curve and annotate.

par(oma = c(6,3,0,0))
plot(`Depth (mm)` ~ duration.hour, data = au.rain.max, 
 xlim = c(0.01, 200),
 ylim = c(1, 5000),
 log = 'xy',
 xlab = '',
 ylab = '',
 xaxt = 'n',
 yaxt = 'n',
 pch = 21, 
 bg = 'grey')



min.values <- c(1, 5, 10, 30, 60, 120, 6*60, 12*60, 24*60, 48*60, 96*60, 4*24*60, 7*24*60 )
min.ticks <- min.values/60
min.labels <- format(min.values, big.mark = ',', trim = TRUE)
axis(side = 1, at = min.ticks, labels = min.labels) # label with minutes
mtext(side = 1, text = 'min', line = 1, adj = c(0,1))

hr.labels <- format(min.values/60, big.mark = ',', trim = TRUE, digits = 0)
hr.labels[as.numeric(hr.labels) < 1] <- NA
mtext(side = 1, text = 'hour', line = 3, adj = c(0, 1))

axis(side = 1, at = min.ticks, labels = hr.labels, line = 3) # label with hours

day.labels <- format(min.values/(60*24), big.mark = ',', trim = TRUE, digits = 0)
day.labels[as.numeric(day.labels) < 1] <- NA

axis(side = 1, at = min.ticks, labels = day.labels, line = 6) # label with hours
mtext(side = 1, text = 'day', line = 6, adj = c(0,1))

mtext(side = 1, text = 'Duration', line = 8 )

axis(side = 2, at = axTicks(2), labels = format(axTicks(2), big.mark = ','), las = 2)

mtext(side = 2, text = 'Rainfall (mm)', line = 4 )



# Add envelope curve

env.au <- function(d) {
 364.22 * d^0.452
 }


d.seq.au <- seq(from = 0.033, to = 200, length.out = 100)

lines(d.seq.au, env.au(d.seq.au), lty = 2, col = 'blue', lwd = 2)
text(0.06, 300, expression(paste('R = 364.2 ',D^0.452, ' (D in hours)')))


Code is also available as a gist.

Bibliography

Arekhi, S. (2013)  Preparing the envelop curve of precipitation (case study: Semnan province). International Journal of Agriculture and Crop Sciences.  5(22): 2768-2772. (link to pdf)

Arekhi, S. (2012) A study to determine the envelope curve of precipitation (case study: Semnan Province). Advances in Environmental Biology 6(1): 210-214. (link to pdf)

Griffiths, G. A., McKerchar, A. I. and Pearson, C. P. (2014) Towards prediction of extreme rainfalls in New Zealand.  Journal of Hydrology (New Zealand) 53(1):41-52. (link to abstract at Informit)

Hershfield, D. M. (1961) Estimating the probable maximum precipitation.  J. Hyd. Div. Am. Soc. Civil Engineering 87: 99-116

Wang, G., L. I. Bao-gui and J.-l. Wang (2006) World’s greatest known point rainfalls and their enveloping curve formula. Advances in Water Science 17(6): 824-829. (link to journal)

WMO (2009) Manual on estimation of probable maximum precipitation (PMP). WMO-No. 1045.  World Meteorological Organization. (link to pdf)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s