5  Reach and catchment environmental metrics

I have calculated many environmental metrics for all reaches of the stream network. These metrics have many potential uses for stream management and research. The remaining chapters of this manual describe the derivation of the metrics and the structure of the tables holding them.

Two primary environmental metric tables hold the most commonly used metrics:

The two tables contain site and reach fields matching those fields in the subcs table. The tables are broadly separated into variables that describe the reach or its local subcatchment (subc_env) and those that describe the catchment upstream of the reach (cat_env). Metrics fall into the following classes (some of which have fields in both subc_env and cat_env.). See Figure 1.1 for the full list of fields in these two tables.

There are two almost identical versions of subcatchment (and catchment) area (in m2) in the database. area_r_m2 in the subc_env table is the sum of gridcell areas for each subcatchment in the r_site raster table. scarea in the subcs table is the area of the polygons in the table. They differ because a small number of subcatchment polygons were manually altered to correct network relationships, resulting in an imperfect match with the DEM gridcells on which the subcatchments were originally calculated. scarea is a more accurate estimate of subcatchment area, but area_r_m2 was used for calculating subcatchment statistics (e.g. impervious cover and tree cover fractions) from raster data. cat_area_m2_r in the cat_env table sums area_r_m2 for all upstream subcatchments, and carea_km2 in the subcs table similarly sums and converts to km2 scarea.

5.1 Impervious cover

Impervious area as a fraction of the upstream catchment (total imperviousness, TI) is a primary measure of urban density. Effective imperviousness (EI) is an important measure of urban stormwater drainage impacts (Walsh et al. 2022). I have calculated these metrics for 2005, 2018, and 2022 (with 2012 and 204 to follow). The fields ti_XXXX and ei_XXXX in the cat_env table record these two metrics for year XXXX. The fields tia_XXXX_m2 and eia_XXXX_m2 in the subc_env table record the total and effect area of impervious surfaces in m2 within each subcatchment. See Chapter 6 for full details.

5.2 Tree canopy cover

Tree canopy cover near streams is a strong predictor of in-stream biotic assemblages (Peterson et al. 2010; Walsh and Webb 2014), and catchment tree canopy cover is an important predictor of streamflow volumes (Zhang, Dawes, and Walker 2001). I thus calculated both tree canopy area as a fraction of stream catchment (TF), and tree canopy area weighted by flow distance to and along streams as a fraction of similarly weighted catchment area (Attenuated forest cover, AF). I have calculated these metrics for 1750 (estimated from ecological vegetation classes), 2006, 2018 and 2022 (with 2012 and 2014 to follow). The fields “tf_XXXX” and “af_XXXX” in the cat_env table, and “tfa_XXXX_m2” (total area of tree canopy cover per subcatchment) in the *subcat_env” table follow the same convention as for ti and ei. See Chapter 7 for full details.

5.3 Flow

Two estimates of mean annual runoff were calculated: post-1997, and for the period 1911-2009 for equivalence with the measure of mean discharge used in earlier studies of the region. See Chapter 8 for full details.

5.4 Flow stress: degree of regulation

A “degree of regulation” index (DoR) was calculated as a measure of flow stress. See Chapter 9 for full details.

5.5 Barriers to fish passage

Barriers to fish (and potentially other fauna) migration such as weirs, dams, and piped streams can limit the distributions of species. The number of barriers to instream migration downstream of each reach was calculated from a revised dataset of barriers included in the database. Two variants were calculated n_all_barriers_ds summing partial and full barriers, and n_full_barriers_ds summing just full barriers were added to the subc_env table. See Chapter 10 for full details.

5.6 Septic tank density

Septic tanks can be major contributors to nitrogen export from catchments (Bernhardt et al. 2008), are their density is a strong predictor of nitrogen concentrations in Melbourne’s streams (Walsh and Kunapo 2009). Sewer backlog programs around peri-urban Melbourne are aiming to replace septic systems with reticulated sewerage (Victorian Auditor-General’s Office 2006), but septic tanks are common throughout the region, even in some areas with reticulated sewerage.

I calculated catchment septic tank density (sep_per_km2) for a portion of the region and added to the cat_env table. Data were only available for the Yarra Ranges Council area (data source: Yarra Ranges Council from around the year 2000), with a small amount of further data for the Dobsons Creek catchment (data source: Yarra Valley Water 2017), and augmented estimates for parts of the Menzies Creek catchment that fall outside the Yarra Ranges Council (Figure 5.1).

Figure 5.1: The stream network colored by septic tank density for reaches where it was calculable (grey elsewhere).

5.7 Catchment geology

Many measures of catchment geology or soil type are possible, using Government spatial data such as the Geomorphological map units of Victoria (Rees et al. 2017) and Victorian Soil type mapping (Anon. 2018). Rather than generate a large number of geological or soil metrics, I calculated two example metrics from the geormophological map: proportion of catchment area with basalt geology (c_basalt), and granitic geology (c_granitic), both in the cat_env table, and provide the code I used to calculate them as a template for other similar calculations (including any polygon land-cover data for the region). The resulting distributions of catchment geology are show in Figure 5.2.

# GMU250.shp downloaded from data.vic.gov.au and stored in geol_dir
# db_m is a connection to the database as in section 1.2
geol <- terra::vect(paste0(geol_dir,"GMU250.shp"))
geol <- terra::project(geol, "epsg:7855")
geol <- terra::makeValid(geol)
# unique(geol$LITHOLOGY) #15 classes including NAs...
# sum(is.na(geol$LITHOLOGY))  #12 - all in MW region are waterbodies, so..
geol$LITHOLOGY[is.na(geol$LITHOLOGY)]<-'Water Body'
# combine...
lith_groups <- unique(geol$LITHOLOGY)
geol$lith_id <- match(geol$LITHOLOGY,lith_groups) # Basalt = 1; Granitic = 5.
r_site <- terra::rast("directory_of_choice/mwstr_v13/r_site.tif")
geol_r <- terra::rasterize(geol, r_site, field = "lith_id") # ~ 10 s
}) # ~30 s
#### Note the next zonal calculations require ~ 70 Gb RAM
basalt_site <- terra::zonal(geol_r, r_site, 
                          fun = function(x){sum(x == 1, na.rm = TRUE)}) 
}) # ~ 4.5 min
granite_site <- terra::zonal(geol_r, r_site, 
                          fun = function(x){sum(x == 5, na.rm = TRUE)}) # ~ 4.5 min
# Assume an active connection to mwstr called db_m
subcs_geol <- DBI::dbGetQuery(db_m,
          paste0("SELECT subcs.site, subcs.nextds, subcs.nus, ", 
                 "subcs.agg_order,subcs.scarea, cat_env.cat_area_m2_r ",
                 "FROM subcs JOIN cat_env ON subcs.site = cat_env.site;"))
subcs_geol <- subcs_geol[order(subcs_geol$nus, subcs_geol$agg_order),]
subcs_geol$sc_basalt <- 
  basalt_site$lith_id[match(subcs_geol$site, basalt_site$site)]
subcs_geol$sc_granite <- 
  granite_site$lith_id[match(subcs_geol$site, basalt_site$site)]
subcs_geol$c_granite <- subcs_geol$c_basalt <- NA
# All subcatchments of agg_order 1 have no subcatchments upstream. Their 
# catchment area (in km2) = their subcatchment area (in m2 * 1e-6)
subcs_geol$c_basalt[subcs_geol$agg_order == 1] <- 
  subcs_geol$sc_basalt[subcs_geol$agg_order == 1]
subcs_geol$c_granite[subcs_geol$agg_order == 1] <- 
  subcs_geol$sc_granite[subcs_geol$agg_order == 1] 
  for(i in min(which(subcs_geol$agg_order == 2)):nrow(subcs_geol)){
# Loop through all remaining subcatchments, adding its subcatchment area to the 
# (already calculated) catchment areas of subcatchments immediately upstream
    subcs_geol$c_basalt[i] <-
      sum(subcs_geol$c_basalt[subcs_geol$nextds == subcs_geol$site[i]], 
    subcs_geol$c_granite[i] <-
      sum(subcs_geol$c_granite[subcs_geol$nextds == subcs_geol$site[i]], 
    if(i %% 10000 == 0) cat(i,"\n")
}) # 1.5 min
# multiply by 25 because each gridcell is 25 m2
subcs_geol$c_granite <- subcs_geol$c_granite*25/subcs_geol$cat_area_m2_r
subcs_geol$c_basalt <- subcs_geol$c_basalt*25/subcs_geol$cat_area_m2_r

Figure 5.2: The stream network colored by A. % of catchment with Basalt geology and B. % of catchment with granitic geology

5.8 Air temperature

Mean daily air temperature at a reach is a strong predictor of mean daily stream temperature. While the diel range of air temperature is greater than stream water, temperature of both air and water tend to vary around a similar value given equivalent shading. Other factors such as shading by tree cover, slope, aspect, flow volume, and groundwater, stormwater and industrial inputs can influence stream temperature. But mean daily air temperature provides a good measure of stream temperature in the absence of anthropogenic alterations to effects such as tree cover, stormwater drainage and flow modification. It also has the advantage of being compatible with a measure by which climate-warming predictions are communicated.

Like geological variables, there are many possible temperature metrics that could be calculated for the network. Unlike the geological data, the source data from the the Australian Bureau of Meteorology Australian Gridded Climate Data (AGCD: Australian Bureau of Meteorology, n.d.) is raster-based. I thus provide the code I used to calculate mean temperature, which can be used as a template for calculating other metrics.

The following code requires specification of a directory grid_dir used to store the rasters, and a connection to the mwstr database, db_m.

# 1. Download all available AGCD data (daily maximum tmax and minimum tmin)
#    Each .nc file is a stack of 365/6 daily rasters for a year.
for(i in 1910:2022){
                     i, ".nc"),
              paste0(grid_dir, "tmax/agcd_v1-0-1_tmax_mean_r005_daily_",i, ".nc"))
                     i, ".nc"), 
              paste0(grid_dir, "tmin/agcd_v1-0-1_tmin_mean_r005_daily_",i,".nc"))

# 2. a) Crop each raster stack to the mwstr region,
#    b) calculate the mean daily maximum and mean for each year
#    c) and add to a stack of 131 rasters of annual means,
#       (one raster for each year)
region <- sf::st_read(mwstr, "region_boundary")
region <- sf::st_transform(region, crs = 4326)
mw_ext <- terra::ext(region)

for(i in 1910:2022){
    tmini <- terra::rast(paste0(grid_dir, "tmin/agcd_v1-0-1_tmin_mean_r005_daily_",i,".nc"))
    tmini <- terra::crop(tmini,mw_ext)
    tmaxi <- terra::rast(paste0(grid_dir, "tmax/agcd_v1-0-1_tmax_mean_r005_daily_",i,".nc"))
    tmaxi <- terra::crop(tmaxi,mw_ext)
    if(i == 1910){
      mw_mean_max_annual_grid <- terra::mean(tmaxi)
      mw_mean_min_annual_grid <- terra::mean(tmini)
# At this stage, seasonal means could be calculated by subsetting, e.g:
#    mw_mean_max_summer_grid <- mean(subset(xi,c(1:60,335:dim(xi)[3])))
      mw_mean_max_annual_grid <- c(mw_mean_max_annual_grid,
      mw_mean_min_annual_grid <- c(mw_mean_min_annual_grid,
})  # 24 min
names(mw_mean_max_annual_grid) <- 1910:2022
names(mw_mean_min_annual_grid) <- 1910:2022
# 3. Estimate the annual means of daily mean temperature as mean of min and max
x <- terra::sds(list(r_tmin_ann, r_tmax_ann))
r_tmean_ann <- terra::app(x, mean)

# 4. Calculate long-term mean raster using the full record of 113 years
r_meant <- terra::mean(r_tmean_ann, na.rm = TRUE)

# 5. calculate 30-year means: 
#    a) set indices
y30_1950 <- which(names(r_tmean_ann) %in% 1921:1950)
y30_1980 <- which(names(r_tmean_ann) %in% 1951:1980)
y30_2010 <- which(names(r_tmean_ann) %in% 1981:2010)
y30_2022 <- which(names(r_tmean_ann) %in% 1993:2022)

#    b) calculate 30-y mean rasters 
r_y30_1950 <- terra::sds(list(r_tmin_ann[[y30_1950]], r_tmax_ann[[y30_1950]]))
r_tmean_y30_1950 <- terra::app(r_y30_1950, mean)
r_meant_y30_1950 <- terra::mean(r_tmean_y30_1950, na.rm = TRUE)

r_y30_1980 <- terra::sds(list(r_tmin_ann[[y30_1980]], r_tmax_ann[[y30_1980]]))
r_tmean_y30_1980 <- terra::app(r_y30_1980, mean)
r_meant_y30_1980 <- terra::mean(r_tmean_y30_1980, na.rm = TRUE)

r_y30_2010 <- terra::sds(list(r_tmin_ann[[y30_2010]], r_tmax_ann[[y30_2010]]))
r_tmean_y30_2010 <- terra::app(r_y30_2010, mean)
r_meant_y30_2010 <- terra::mean(r_tmean_y30_2010, na.rm = TRUE)

r_y30_2022 <- terra::sds(list(r_tmin_ann[[y30_2022]], r_tmax_ann[[y30_2022]]))
r_tmean_y30_2022 <- terra::app(r_y30_2022, mean)
r_meant_y30_2022 <- terra::mean(r_tmean_y30_2022, na.rm = TRUE)

# # (Inspect a raster, adding the coastline for orientation)
# terra::plot(r_meant_y30_2022) 
# coast <- sf::st_read(db_m, "coast")
# coast <- terra::vect(sf::st_transform(coast, crs = 4326))
# terra::plot(coast, add = TRUE)

# 6. extract the 30-year mean values for every subc pourpoint
rmeant_30y <- c(r_meant_y30_1950,r_meant_y30_1980,
names(rmeant_30y) <- paste0("meant_y30_",c(1950,1980,2010,2022))
subc_pps <- sf::st_read(db_m, 'subc_pps')
subc_pps <- sf::st_transform(subc_pps, crs = 4326)
subc_pps <- terra::vect(subc_pps)
x <- suppressWarnings(terra::extract(rmeant_30y, subc_pps))
x <- data.frame(site = subc_pps$site, x[,-1])

The above code calculates the mean daily air temperature at the pour-point (lowest point) of each subcatchment in the network as the mean of the minimum and maximum daily temperature grids from the AGCD (Australian Bureau of Meteorology, n.d.), for four 30-year periods: to 1950, 1980, 2010 and 2022.

Across the region, mean temperatures during each of these periods were highly correlated, but warming of ~1 C\(^\circ\) over the decade is evident (Figure 5.3).

Figure 5.3: Relationship between the mean air temperature 1910-1922 to the mean temperature over four 30-year periods for every gridcell in the 5-km resolution Australian Gridded Climate Data cropped to the extent of the Port Phillip and Westernport Region.

The four means were added to the subc_env table as meant_y30_1950, meant_y30_1980 meant_y30_2010, and meant_y30_2022. Figure 5.4 shows the distribution of the last of these among ‘sampleable’ streams using the following code.

# db_m is a connection to the database as in section 1.2
streams <- sf::st_read(db_m, 
                       query = "SELECT * FROM streams WHERE sampleable = 1;")
airtemp <- DBI::dbGetQuery(db_m, "SELECT site, meant_y30_2022 FROM subc_env;")
streams$meant_y30_2022 <- 
  airtemp$meant_y30_2022[match(streams$site, airtemp$site)]

Figure 5.4: Mean air temperature from 1993 to 2022 in sampleable streams of the network and code used to create the plot

5.9 Elevation, slope and aspect

I calculated elevation (above sea-level in m, subc_env.elev_asl_m, Figure 5.5 A) and slope (percentage, subc_env.slope_perc, Figure 5.5 B) from the filled 5-m digital elevation model (table r_dem). I calculated aspect (as northness, Figure 5.5 C) from the stream lines. These three reach attributes are stored in the subc_env table.

The stream lines are accurately delineated throughout the region, having been derived from a 1-m-resolution DEM, and a more accurate estimate of slope could be achieved using that DEM. I took a hybrid approach to estimating slope of each reach by calculating the difference in elevation between the 5-m raster pourpoints for each subcatchment and its upstream subcatchment(s), and used the length of the stream line through the subcatchment as the horizontal distance.

Slope was calculated as a percentage as:

\[ 100 \cdot \frac{\Delta E}{\sqrt{{L}^2 + {\Delta E}^2}} \tag{5.1}\]

where \(\Delta E\) = the difference in elevation between upstream and downstream points and \(L\) is the total length of the stream lines (of all types) in the reach. Slopes were not calculated for headwater reaches, most of which were very short line segments. Negative slopes, of which there were 1672, almost all very small, were set to zero.

Aspect can influence how shaded a reach is with a reach flowing north-south likely to be less shaded by its banks and hillslopes than a reach flowing east-west. A more complete measure of such an effect could integrate shading calculated from a DEM and daily patterns of solar radiation, but northness is offered as an initial approximation. Northness was calculated from the angle between the most downstream and most upstream points of the line(s) of each reach. It was calculated as

\[ abs(cos(\theta)) \tag{5.2}\]

where \(\theta\) is the angle in radians from North between the two points. Thus for a reach flowing north to south (or vice-versa) northness = 1, and for a reach flowing east to west (or vice versa) northness = 0.

# db_m is a connection to the database as in section 1.2
str_slope_elev <- sf::st_read(db_m, query = 
    "SELECT streams.geom, subc_env.elev_asl_m, subc_env.slope_perc, 
     subc_env.northness, streams.reach, streams.sampleable, subcs.agg_order 
      FROM streams JOIN subc_env ON streams.site = subc_env.site 
      JOIN subcs on subc_env.site = subcs.site;")
# Order by agg_order to place largest streams on top of plot
str_slope_elev <- str_slope_elev[order(str_slope_elev$agg_order),]
# Set colors for maps
str_slope_elev$elev_class <- 
                 breaks = c(-10,10,20,40,80,160,320,640,1340)))
str_slope_elev$elev_col <- 
str_slope_elev$slope_class <- 
  as.numeric(cut(str_slope_elev$slope_perc, breaks = c(0,1,2,4,8,16,32,75)))
str_slope_elev$slope_col <- 
str_slope_elev$north_class <- 
                 breaks = c(-0.1,0.2,0.4,0.6,0.8,1.1)))
str_slope_elev$north_col <- 
# Plot maps
par(mfrow = c(2,2), mar = c(0,0,0,0))
plot(str_slope_elev$geom, col = str_slope_elev$elev_col)
title(main = "A.", adj = 0, line = -1)
legend("bottomleft", lty = 1, col = RColorBrewer::brewer.pal(8, "Spectral"), 
  legend = c("<10","10-20","20-40","40-80","80-160","160-320","320-640",">640"),
  cex = 0.7, title = "Elevation (m)")
plot(str_slope_elev$geom, col = str_slope_elev$slope_col)
title(main = "B.", adj = 0, line = -1)
legend("bottomleft", lty = 1, col = RColorBrewer::brewer.pal(7, "Spectral"), 
  legend = c("<1","1-2","2-4","4-8","8-16","16-32",">32"),
  cex = 0.7, title = "Slope (%)")
plot(str_slope_elev$geom, col = str_slope_elev$north_col)
title(main = "C.", adj = 0, line = -1)
legend("bottomleft", lty = 1, col = RColorBrewer::brewer.pal(5, "Spectral"), 
  legend = c("0-0.2","0.2-0.4","0.4-0.6","0.6-0.8","0.8-1,0"),
  cex = 0.7, title = "Northness")

Figure 5.5: The stream network colored by A. elevation (square-root-transformed), B. stream slope (4th-root-transformed), and C. northness code used to produce the maps

5.10 Strahler stream order

I calculated Strahler order (variable strahler in subc_env, a numerical measure of branching complexity) for the reaches of the network using the following code, which applies the algorithm described at https://en.wikipedia.org/wiki/Strahler_number.

subcs_s <- subcs[c("site","nextds","nus","agg_order")] 
subcs_s <- subcs_s[order(subcs_s$nus, subcs$agg_order),]
subcs_s$strahler <- NA
subcs_s$strahler[subcs_s$nus %in% c(0)] <- 1
# If there are no upstream lines strahler = 1
  for(i in min(which(is.na(subcs_s$strahler))):dim(subcs_s)[1]){
    if(sum(subcs_s$nextds == subcs_s$site[i]) == 1){
# If the site is a split in a stream line not at a confluence make strahler the 
# value of the site upstream
      subcs_s$strahler[i] <- subcs_s$strahler[subcs_s$nextds == subcs_s$site[i]]
      usi <- subcs_s[subcs_s$nextds == subcs_s$site[i],]
      subcs_s$strahler[i] <-  ifelse(sum(usi$strahler == max(usi$strahler)) == 1,
# If the node has one child with Strahler number i, and all other children have 
# Strahler numbers less than i, then the Strahler number of the node is i again.
# If the node has two or more children with Strahler number i, and no children 
# with greater number, then the Strahler number of the node is i + 1.
                                     max(usi$strahler) + 1)
  } # took ~3 minutes

To plot a map showing spatial patterns of Strahler order, run the following code (or a variant thereof as explained in (Appendix C)[app3-db_assembly.qmd])

# db_m is a connection to the database as in section 1.2
str_strahler_subcat <- sf::st_read(db_m, query = 
      paste0("SELECT streams.geom, subc_env.strahler ",
             "FROM streams JOIN subc_env ON streams.site = subc_env.site ",
             "WHERE streams.sampleable = 1;"))

5.11 Melbourne Water’s ‘sub-catchments’

An administrative measure rather than an environmental one, Melbourne Water’s ‘subcatchments’ are spatial management units used for stream management reporting, including environmental condition. I thus included them as a field in subc_env, linked to the table mw_subcatchments

The Melbourne Water Healthy Waterways Strategy (HWS) 2018 divides the region into 69 management units, which MW call “subcatchments” (although many are not hydrological subcatchments). To permit reporting for the HWS to these spatial units, each site was allocated to a “subcatchment” (the mw_subcat field) and a “catchment” (the mw_cat field). To allocate each site, I identified the most downstream sites for MW subcatchments in QGIS (saved in terminal_segs_mw_scs).

The numeric identifiers mw_subcat and mw_cat correspond to 69 ‘subcatchment’ names, and 8 ‘catchment’ names. These ‘subcatchments’ and ‘catchments’ are defined by polygons that roughly correspond to the stream network (from 1). The table mw_subcatchments contains those polygons, with fields mw_subcat, mw_cat, mw_subcat_nameand mw_cat_name, to permit matching of numeric mw_subcat and mw_cat codes to names.

To plot a map showing spatial patterns of mw_subcat and mw_subcat_name, run the following code (or a variant thereof as explained in (Appendix C)[app3-db_assembly.qmd])

# db_m is a connection to the database as in section 1.2
str_mwcat_subcat <- sf::st_read(db_m, query = 
       paste0("SELECT streams.geom, subc_env.mw_subcat, subc_env.mw_cat ",
                "FROM streams JOIN subc_env ON streams.site = subc_env.site ",
                "WHERE streams.sampleable = 1;"))

5.12 Australian postcodes

A second administrative measure that may be used to narrow searches of the database, a postcode was assigned to (the bottom of) each reach by joining the polygons of postcode areas from Australia Post with table subc_pps. Twenty reaches (mainly estuaries, pipe outlets and a few coastal tributaries on French Island) were not assigned a postcode.

To plot a map showing spatial patterns of postcodes, run the following code (or a variant thereof as explained in Appendix C)

# db_m is a connection to the database as in section 1.2
str_postcode_subcat <- sf::st_read(db_m, query = 
           paste0("SELECT streams.geom, subc_env.postcode ",
                  "FROM streams JOIN subc_env ON streams.site = subc_env.site ",
                  "WHERE streams.sampleable = 1;"))

  1. see Grant, A. (2018) Subcatchments - Healthy Waterways Strategy 2018 -Data guide and metadata. A copy of this document and the data it describes resides on the unimelb drive “wergSpatial/Projects/MW_StreamNetwork/env_compilation”↩︎