Category Archives: R

Flood frequency plots using ggplot

This post provides a recipe for making plots like the one below using ggplot2 in R.  Although it looks simple, there are a few tricky aspects:

  • Superscripts in y-axis labels
  • Probability scale on x-axis
  • Labelling points on the x-axis that are different to the plotted values i.e. we are plotting the normal quantile values but labelling them as percentages
  • Adding a title to the legend
  • Adding labels to the legend
  • Positioning the legend on the plot
  • Choosing colours for the lines
  • Using commas as a thousand separator.



Code is available as a gist, which also shows how to:

  • Enter data using the tribble function, which is convenient for small data sets
  • Change the format of data to one observation per row using the tidyr::gather function.
  • Use a log scale on the y-axis.

Links for more information:

Modelling impervious surfaces in RORB – II

This blog builds on the previous post; looking at the runoff coefficient approach to modelling losses and the implications for representing impervious surfaces in the RORB model.

In addition to the IL/CL model discussed in the previous post, RORB can be run using an initial loss / runoff coefficient model, where the runoff coefficient specifies the proportion of rainfall lost in each time step after the initial loss is satisfied.  This reason these different loss models are of interest is that the new version of Australian Rainfall and Runoff is recommending that the IL/CL model is used in place of the runoff coefficient model (Book 5, Section 3.3.1).  In some areas, modelling approaches will need to change and this will have implications for flood estimates.

The runoff coefficient loss model is selected as shown in Figure 1.


Figure 1: A runoff coefficient loss model can be selected in RORB

The user inputs the runoff coefficient, C, for a pervious surface.  For an impervious surface, there is no opportunity to specify the runoff coefficient which is hard-wired in RORB as 0.9.  For mixed sub-areas, the runoff coefficient is scaled, the equations from the RORB manual are:

C_i = F_iC_{imp} +(1-F_i)C_{perv}, \qquad C_{perv} \le C_{imp} \qquad \mathrm{Equation  \;3.5}
C_i = C_{imp}, \qquad C_{perv} > C_{imp}\qquad\qquad \mathrm{Equation \; 3.6}

Where Ci is the runoff coefficient for the ith sub-area.

Example: For a fraction impervious, F_i = 0.6 and C_{perv} = 0.5
C_i = 0.6 \times 0.9 +(1-0.6) \times 0.5 = 0.74

The initial loss is calculated as as a weighted average of the pervious and impervious initial losses as shown in the previous post.  The impervious initial loss is always set to zero in RORB.

Let’s do the calculations for a 100% impervious surface.  RORB will set I\!L = 0 and C = 0.9.  Using the 6 hour, 1% rainfall as before, the rainfall excess hyetograph is shown in Figure 2.


Figure 2: Rainfall excess hyetograph for an impervious surface using the runoff coefficient model.  RORB sets IL to zero and the the runoff coefficient to 0.9 so 10% of rain is lost at each time step

Example calculation:

As explained in the previous post, the rainfall between 1.5 hour and 2 hour is 19.4 mm.  With a runoff coefficient of 0.9, the rainfall excess will be: 0.9 x 19.4 = 17.5 mm.

The rainfall excess hydrograph from a 10 km2 impervious sub-area can be calculated from the rainfall excess hyetograph using the method described in the previous post. The peak flow corresponding to the 17.5 mm rainfall peak is 97.2 m3s-1 (see the previous post for sample calculations).

The key point is that we have changed the peak flow from an impervious surface, just by changing the loss model.  With the IL/CL model, both initial and continuing loss for a 100% impervious surface are hard-wired to zero. The peak runoff was 107.8 m3s-1. For the runoff coefficient model, initial loss is hard-wired to zero, but the runoff coefficient is hard-wired to 0.9, i.e. we have some loss from the impervious surface. This changes the hydrograph as shown in Figure 3.


Figure 3: Comparison of rainfall excess hydrographs from a 100% impervious surface; same rainfall, different loss model

The value of the runoff coefficient for an impervious surface is noted in the RORB manual:

The impervious area runoff coefficient Cimp is set by the program to 0.9, reflecting the fact that losses occur even on nominally impervious surfaces in urban areas.

This is reasonable, but inconsistent with the treatment of the continuing loss when the IL/CL loss model is used.  In this case, CL is hard-wired to zero so there are no losses from impervious surfaces; a feature of RORB for modellers to be aware of.

Also note Equation 3.6 above.  This suggests that if the user inputs a runoff coefficient larger than the impervious coefficient (i.e. larger than 0.9) then a value of 0.9 will be used.  This isn’t actually implemented.  If a runoff coefficient of 1 is input, there is a direct conversion of rainfall to runoff i.e. there is no loss.  It is even possible to input runoff coefficients greater than 1.

Equation 3.6 may just be the result of a typo.  Some experimenting suggests the behaviour in the model is represented the combination of equation 3.5, above and the following in place of equation 3.6:

C_i = C_{perv}, \qquad C_{perv} > C_{imp}

That is, the runoff coefficient for an impervious surface is 0.9 unless the runoff coefficient input by the user is larger than 0.9.

Calculations are available via a gist.

Modelling impervious surfaces in RORB

The previous post looked at rainfall excess hydrographs; here I explore how these hydrographs change when modelling impervious surfaces in RORB.  This post focusses on the initial loss/continuing loss modelling approach.

Usually, losses are reduced for impervious compared to pervious surfaces and RORB sets both initial and continuing loss to zero if a surface is 100% impervious.

As an example, consider the 6 hour 1% rainfall for Melbourne, which is 83.4 mm.  If we use the ARR1987 temporal pattern (see the previous post), the hyetograph is as shown in Figure 1.


Figure 1: The 6 hour 1% rainfall (83.4 mm) multiplied by the ARR1987 temporal pattern.  For an impervious surface, RORB sets both initial and continuing loss to zero

Example calculation:

In the ARR1987 temporal pattern, the time period between 1.5 and 2 hours has 23.3 percent of the rain.  The total rainfall is 83.4 mm so the rain in this period is 83.4 x 23.3% = 19.43 mm which is consistent with Figure 1.

The corresponding rainfall excess hydrograph, for an area of 10 km2, which is 100% impervious, is shown in Figure 2 (Note that Areal Reduction Factors have not been used).


Figure 2: rainfall excess hydrograph for an area of 10 km2

Example calculation:

The instantaneous flow at a 2 hours will be

\frac{1}{3.6} \times \frac{1}{0.5} \times 19.4 \times 10 = 107.8 \mathrm{m^3 s^{-1}}

To explain factors at the start of the equation, 1/3.6 is for unit conversion, 1/0.5 is because the temporal pattern has a 0.5 hour time step.   The RORB output matches the calculations (Figure 3).

By default, RORB will show the rainfall excess hyetograph above the calculated hydrograph but this is  based on the initial and continuing loss as provided by the user.  In this case, I’ve specified IL = 10 mm and CL = 2 mm/h for the pervious areas.   These losses, and the hyetograph, are misleading where a sub-catchment has some impervious component.  In this case, for a 100% impervious sub-area, both IL and CL are set to zero by the program.  It would be best not to display the misleading hyetograph, which can be turned off as shown in Figure 4.


FInure 3: RORB output


Figure 4: The hyetograph can be toggled off using the button outlined in pink

If a sub-area is a combination of both impervious and pervious surfaces, this must be specified to RORB as a Fraction Impervious (Fi).  The initial and continuing losses are scaled based on this fraction.

IL_i = (1 - F_i) I\!L_{perv}
CL_i = (1 - F_i) C\!L_{perv}

Where I\!L_{perv} and C\!L_{perv} are the initial and continuing losses for pervious areas as input by the user.

For example, if the pervious value of IL is set to 10 mm and CL to 2 mm/h, then for a sub-area with a Fraction Impervious value of 60%, the initial and continuing losses will be:

I\!L = (1 - 0.6) \times 10 = 4 \; \mathrm{mm}
C\!L = (1 - 0.6) \times 2 = 0.8 \; \mathrm{mm/h}

The continuing loss is 0.8 mm/h which is 0.4 mm per 30 min time step.

Running the model with these parameters results in a rainfall excess hydrograph as shown in Figure 5.  Note that the start of the rise of the hydrograph is delayed because of the initial loss.  The peak is reduced by a small amount (from  108 cumec to 106 cumec because of the continuing loss).

Example calculation, flow peak:

\frac{1}{3.6} \times \frac{1}{0.5} \times (19.4 - 0.4) \times 10 = 105.6 \mathrm{m^3 s^{-1}}


Figure 5: Rainfall excess hydrograph

For a real catchment with a 60% fraction impervious, we would expect some early runoff from the impervious surfaces that would provide flow directly into the urban drainage system.  RORB doesn’t model this process, which may not matter, depending on the application, but as modellers we need to be aware of this limitation.

Calculations are available as a gist.

Rainfall excess hydrograph

Ever wondered what a ‘rainfall excess’ hydrograph is and how they are calculated?  Then read on.

‘Rainfall excess’ is the rainfall left over after the initial and continuing loss are removed.  Rainfall excess hydrographs are used in the runoff-routing program RORB.  The RORB manual (Section 3.3.4) describes them as follows:

In catchment studies, the program calculates hyetographs for all sub-areas.  After deducting losses, it converts the hyetograph ordinates to ‘hydrographs’ of rainfall-excess on the sub-areas, in m3/s, and interprets the average ‘discharge’ during a time increments as an instantaneous discharge at the end of the time increment.

Lets look at an example.  I’m using the methods from the 1987 version of Australian Rainfall and Runoff so I can compare results with calculations in RORB.

1. Choose a design rainfall depth

I’m working on a catchment in Gippsland where the 1% AEP 6 hour rainfall is of interest.  Rainfall IFD data is available from the Bureau of Meteorology via this  link.

For the site of interest, the 1% (100-year), 6 hour rainfall depth is 90.9 mm.

2. Select a temporal pattern

Temporal patterns are available in Australian Rainfall and Runoff Volume 2, Table 3.2.  Gippsland is in zone 1 and ARI is > 30 years so we need the bottom row from the table below.   This shows the percentage of the rainfall depth in each 30 min time period


Applying the temporal pattern to the design rainfall depth results in the following hyetograph.


Figure 1: Design rainfall hyetograph

3. Remove the losses

Calculate the rainfall excess hyetograph by removing the initial loss and continuing loss.  For this example,

  • IL = 10 mm and
  • CL = 2 mm/h.

Note that the continuing loss is 2 mm/h and the time step of the hyetograph is 0.5 h so 1 mm is lost per time step.

The rainfall excess hyetograph is shown in Figure 2.


Figure 2: Rainfall excess hyetograph

4. Convert to a hydrograph

The procedure to convert a rainfall excess hyetograph to a rainfall excess hydrograph is explained in the quote at the start of the blog.  We need to:

  • Multiply the rainfall excess by the catchment area (converts rainfall to a volume)
  • Divide by the time step (to calculate volume per unit time)
  • Ensure flow is allocated to the correct time step – the rainfall during a time step produces the instantaneous flow at the end of the time step
  • Ensure the units are correct – calculated flow is is m3/s, rainfall is in mm and catchment area is in km2.

There is also a discussion of this in ARR2016 Book 5, Chapter

Example calculation: in this case, the sub-catchment area is 78.7 km2.  The rainfall in the 3rd time step,  between 1 hour and 1.5 hour, is 8.9 mm so the flow at the end of this time step will be:

Q = \frac{1}{3600} \times  10^{-3} \times 10^6 \times \frac{8.999}{0.5} \times 78.7 = 393.46 \; \mathrm{m^3s^{-1}}

The rainfall excess hydrograph is shown in Figure 3.


Figure 3: Rainfall excess hydrograph

4. Comparison with RORB

Figure 4 shows the rainfall excess hydrograph as calculated by RORB.  The answers look close and I’ve confirmed this by looking at the calculated values.


Figure 4: Rainfall excess hydrograph as calculated by RORB

Calculations are available as a gist.

On the calculation of equal area slope

As noted in the previous post, the equal area slope was adopted for use with the Bransby Williams time of concentration formula in the 1987 version of Australian Rainfall and Runoff to: “…give a better indication of flow response times, especially where there are large variations of slope within a catchment” (ARR 1987 Boo IV, Section 1.3.2(d)).  The equal area slope is also used as part of flood estimation in New Zealand (NZ, 1980; Auckland Regional Council, 1999), in the Papua New Guinea Flood Estimation Manual (SMEC, 1990) – where it is used in the estimation of overland flow times and runoff coefficients for use in the rational method – and is discussed in the Handbook of Hydrology (Pilgrim and Cordery, 1993).

The equal area slope is the slope of a straight line drawn on a profile of a stream such that the line passes through the outlet and has the same area under and above the stream profile.

An alternative to the equal area slope, as used by Bransby Williams, is the average slope.  The differences between the equal area and average slopes are highlighted in the Figure 1 below.


Figure 1: Average slope and equal area slope (McDermott and Pilgrim, 1982, page 28)

I haven’t been able to find the history of the equal area slope.  I’m guessing that the average slope was found to be too steep for use in hydrologic calculations.  The equal area slope may have been considered more representative and was easy to calculate in pre-computer times.  The procedure, perhaps to be undertaken by a draftsman, is  specified in NZ (1980) (see Figure 2).

The method involves the calculation of the slope of the hypothetical line AC, which is so positioned that the enclosed areas above and below it, i.e. areas X and Y, are equal. The procedure is to planimeter the total area under the longitudinal profile.  This area Ad, equals the area of the triangle ABC.


Figure 2: Definition diagram for the calculation of the equal area slope (NZ, 1980)

A_d = \frac{1}{2}AB \times BC

A_d = \frac{1}{2}L \times h

\therefore h= \frac{2A_d}{L}

(Point C is known at the ‘equal area slope ordinate’).

Hence the equal area slope S_a is given by

S_{ea} = \frac{h}{L} = \frac{2A_d}{L^2}

When the units for the elevation and length in the diagram above are used:

S_{ea} = \frac{2A_d}{1000L^2} \;\; \mathrm{m/m}

Notice that nothing iterative is required, we are just calculating the area under the profile and then working out the triangular area to match.

This simple calculation procedure seems to have been forgotten, as some modern approaches suggest an iterative procedure is necessary, for example, in this excel equal area slope tool.

An R function to calculate the equal area slope is available as a gist.


Auckland Regional Council (1999) Guidelines for stormwater runoff modelling in the Auckland Region. (link)

McDermott, G. E. and Pilgrim, D. H. (1982) Design flood estimation for small catchments in New South Wales.  Department of National Development and Energy.  Australian Water Resources Council.  Research Project No. 78/104. Australian Government Publishing Service, Canberra (link)

Pilgrim, D. H. and Cordery, I. (1993) Flood runoff. In: Maidment, D. R. (ed) Handbook of Hydrology.  McGraw Hill.

SMEC (1990) Papua New Guinea Flood Estimation Manual. Department of Environment and Conservation, Bureau of Water Resources. (link)

NZ (1980)A method for estimating design peak discharge: Technical Memorandum No. 61. New Zealand. Ministry of Works and Development. Water and Soil Division. Planning and Technical Services; National Water and Soil Conservation Organisation (N.Z). (1980)  (link to catalog entry) (link to document)

Load calculations: error approximations

David Fox has written a very useful article about calculating the uncertainty in estimates of sediment and nutrient loads.

A standard expression for the variance of a variable is:

Mean of the square minus the square of the mean

In this case:

\mathrm{Var[L]} = E[L^2] - (E[L])^2

Where L is the load and E is the expectation.

Fox (2005) derives expressions for the required expectations based on the assumption that flow and concentration is described by the bivariate lognormal distribution.

E[L] = \exp{ \{ (\mu_C +\mu_Q) +\dfrac{1}{2(1-\rho^2) } \left( (\rho \sigma_Q+ \sigma_C)^2 + (\rho \sigma_C + \sigma_Q)^2 - 2 \rho(\rho \sigma_Q + \sigma_C)(\rho \sigma_C + \sigma_Q) \right) \}}

E[L^2] = \exp{ \{ 2(\mu_C +\mu_Q) +\dfrac{2}{(1-\rho^2)} \left( (\rho \sigma_Q + \sigma_C)^2 + (\rho \sigma_C + \sigma_Q)^2 - 2\rho(\rho\sigma_Q + \sigma_C)(\rho \sigma_C + \sigma_Q) \right ) \} }

An R function to calculate these values is as follows:

# mQ = mean log(flow)
# mC = mean log(concentration)
# sQ = standard deviation log(flow)
# sC = standard deviation log(concentration)
# r = correlation between log(flow) and log(concentration)

Var_L = function(mQ, mC, sQ, sC, r){
  EL = exp(  (mC + mQ) + (1/(2*(1-r^2))) * ((r*sQ + sC)^2 + (r*sC + sQ)^2 - 2*r*(r*sQ + sC)*(r*sC + sQ))) 
  EL2 = exp( 2*(mC + mQ) + (2/(1-r^2)) * ((r*sQ + sC)^2 + (r*sC + sQ)^2 - 2*r*(r*sQ + sC)*(r*sC + sQ)))
  EL2 - (EL)^2

Sample calculation

Fox (2005) provides data for a sample calculation, we can use this data to test the functions to calculate these expectations.

Parameter log(flow) log(concentration)
mean 2.5561 -0.02834
standard deviation 0.6706 0.8008

Correlation between log(flow) and log(concentration) = 0.482.

Using these, values Var[L] = 3132.297 which is consistent with the Fox’s calculations.

Code is available as a gist.


Fox, D. R. (2005) Protocols for the optimal measurement and estimation of nutrient loads: error approximations. Australian Centre for Environmetrics.  University of Melbourne. link

Hydstra rating plots explained

When viewing discharge data at the Victorian Water Monitoring website, one option is to see a Gaugings vs. Ratings Plot (Figure 1).


Figure 1: A gauging v rating plot can be created for each gauge

An example for the Avon River at Stratford is shown in Figure 2.  This plot may make perfect sense for a hydrographer but was challenging for me to understand, and I imagine many others.

First, check out the y axis label:

  • 100 – Stream Water Level (m) plotted against CTF of 1.225 in metres

Also look at the  labels on the y-axis tick marks:

  • 1.226, 1.235, 2.225, 11.225.

Figure 2: Avon River at Stratford, gauging vs. ratings plot

What this is referring to is the rating equation that has been used to establish the rating curve.  This is likely to be in the form:

Q = C(h - h_o)^m \quad \mathrm{Equation \; 1}

Where: Q is discharge,  C is a constant, h is gauge level,  ho is the cease to flow (CTF) level, m is a constant exponent with typical value 1.5 to 2.5.  These three parameters are fitted based on gauging data.

For the Avon River at Stratford, ho is 1.225 m.  This explains the comment about CTF in the y-axis label.  The ‘100’ in the y-axis label is a value used in Hydstra to refer to water levels. There is a corresponding ‘141’ in the x-axis label, which Hydstra uses to refer to discharge which is created from levels by use of the rating curve.  A list of Hydstra codes is here.

Now to the y-axis labels.  These are based on log-cycles with the addition of a constant representing the CTF level.

  • 101   + 1.225 = 11.225
  • 100  + 1.225 = 2.225
  • 10-1 + 1.225 = 1.325
  • 10-2 + 1.225 = 1.235
  • 10-3 + 1.225 = 1.226

So, Figure 2 is a plot of the following:

  • \log_{10}(Q) on the x axis
  • \log_{10}(L - 1.225) on the y axis.

Where, Q is the discharge and L is the stream water level (the gauge reading).  These transformations relate to equation 1.

Stream gaugings are available from the Australian Bureau of Meteorology website  Rating curves are also available from the website and from the Victorian warehouse

I’ve plotted the gaugings and the rating curves in the Hydstra format in Figure 3 below. This graph is the same as Figure 2 with the added value of colouring points based on when data was collected.  Clearly there is a lot of scatter in the low end gaugings and its not just based on when the data were collected.  There are ‘red’ points (recent gaugings) above and below the current rating.  At higher discharges, it looks like the water levels are increasing for a given flow.


Figure 3: Plot of ratings and gaugings for the Avon River at Stratford.  Colours indicate when gaugings were made.

The upshot is, that it is possible to construct these Hydstra style plots given gauging data and an estimate of the cease to flow level.

Code to create this plot is available as a gist.