Category Archives: R

Graphing a water balance

The water balance for a urban catchment equates the change in storage during a certain period, with the difference between water inputs (precipitation and mains water) and water outputs (evaporation, stormwater runoff and wastewater discharge).

\Delta s = (P+I) - (E_a + R_s + R_w)
where:

\Delta s change in catchment storage
P precipitation
I  imported water
E_a actual evaportranspiration
R_s stormwater runoff
R_w wastewater discharge

Mitchell et al., (2003) provides data on the water balance for Curtin, ACT for 1979 to 1996.  The water balance for the average, wettest and driest years are shown in the table below.

Water-balance_table.png

When presenting financial statements, a common approach is to use a waterfall chart which shows how the components of a financial balance contribute to an overall result.  Here I’ve used a waterfall chart to show the water balance for Curtin for the driest and wettest year as reported by Mitchell et al., (2003).

Water_balance_driest

Water_balance_wettest

Figure 1: Water balance for Curtin, ACT in (A) and driest and (B) the wettest years as estimated by Mitchell et al., (2003).

Does this approach to visualising a water balance help understanding?  A few things stand out:

  • In the driest year, more water was input from the mains than from rainfall
  • In the driest year, actual evapotranspiration was larger than rainfall and mains inputs.
  • Evapotranspiration and stormwater change with climate, with large variation between the wet and dry years.  Wastewater doesn’t change all that much.
  • Precipitation is highly variable, ranging from 247 mm to 914 mm.

There is a guide to making a waterfall chart in excel here.  The R code to produce the graphs shown in this blog is available as a gist, which draws on this blog.

References

Mitchell, V. G., T. A. McMahon and R. G. Mein (2003) Components of the Total Water Balance of an Urban Catchment. Environmental Management 32(6): 735-746. (link)

Munging rating tables

The Victorian water monitoring site includes rating tables for stream gauges but they are in a format that is not easy to work with.   An example is shown in Figure 1 below.

Rating table

Figure 1: Extract of rating table

The following steps can be used to extract and convert the data into a useable format.

1. Download and save rating table.  Click the button shown to get the rating table as a text file.

Rating_table_download

Figure 2: Save the rating table

2. Re-format the data to create columns of levels and flows.  You’ll need to use your favourite tool for this munging step.  An example using R is available as a gist.

3. Plot and compare with the online version

Rating_table_404216_warehouse

Figure 2: Rating plot (source: data.water.vic.gov.au/monitoring.htm?ppbm=404216&rs&1&rscf_org)

Rating_table_404216_R

Figure 3: Rating plot using data from rating table

4. Save as a csv file for further use.

R code is available here.

Related posts:

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.

FFA_plot

 

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
  • Plot a secondary axis showing the AEP as 1 in X years
  • Use the Probit transformation for the AEP values

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.

RORB_RoC

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.

RoC_impervious

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.

Hydro-Comp

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.

TP_Melbourne_6h_1pc

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).

RE_hydro

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.

RORB.png

FInure 3: RORB output

RORB2.png

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}}

Rain_excess.jpeg

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

TemporalPatter_gt30y6h

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

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.

Hyeto_ILCL

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 6.4.3.1.

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.

Rainfall_excess_hydro

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.

Rainfall_excess_hydro_RORB

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.

equal-area-slope-ave-slope

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.

equal-area-slope

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.

References

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)

https://nicgreeneng.wordpress.com/2016/07/03/equal-area-slope-tool/