Assessment of post-wildfire erosion risk and effects on water quality in south-western Australia

Additional keywords: Peel–Harvey estuary, RUSLE, water catchments.

Water quality impacts from wildfire have been documented in catchments globally across various fire-prone regions including south-eastern Australia, the Mediterranean, South Africa, the western United States of America and Canada (see for example White et al. 2006; Rhoades et al. 2011; Smith et al. 2011a; Oliver et al. 2012; Bladon et al. 2014; Dahm et al. 2015). Recent wildfires in eastern Australia (proximal to the population centres of Canberra and Sydney) have demonstrated the potential to elicit adverse water quality responses. Water quality impacts include increased suspended sediment, nutrient concentrations and other constituents associated with fine sediment (Smith et al. 2011a). These impacts can have adverse effects on aquatic ecosystems and challenge water supply systems. In south-eastern Australia, there has been a substantial investment in research programs aimed at building capacity to predict post-fire erosion and associated water quality impacts (AFAC 2017). In fire-prone forests of south-western Australia, this research is lacking such that the impacts of wildfire on the quality of water resources are unknown and land managers have few means to prioritise mitigation efforts.

Sources of pollutants in forested watershed catchments include hillslopes, channels and floodplains as well as different land-use activities and fire, either prescribed or wildfire. Differences in landform, soil, vegetation and rainfall regimes are all likely to result in distinct post-fire response domains in terms of frequency and magnitude of erosion. A key factor is landform and the degree of connectivity between hillslopes and streams. In the eastern uplands of Victoria, for instance, the dissected uplands are characterised by high relief and long, steep slopes that converge in first-order ephemeral drainage lines that are effective at delivering sediment from hillslopes to streams during runoff events (Lane et al. 2006; Smith et al. 2011b). In contrast, the hillslopes in the Hawkesbury sandstone escarpments near Sydney are largely disconnected from streams owing to high infiltration rates in colluvial deposits between cliffs and the stream network (Blake et al. 2009). The Darling Ranges in south-western Australia represent another distinct landform, with dissected low-relief escarpments that drain onto a coastal plain.

The escarpments along the south-western coast of Australia are heavily forested, receive reasonably high rainfall and form important catchments for water supply reservoirs, irrigation and wetlands on the coastal plains. Within the escarpments, the forested catchments deliver high water quality to streams and thus help support aquatic ecosystems that are unique to this region (Bunn and Davies 1990). These high-value catchments have experienced declining water yield since 1975 due to ongoing rainfall deficit that may be linked to climate change in the region (Hope et al. 2006). The decline in surface water has seen the region becoming increasingly reliant on groundwater (McFarlane et al. 2012). However, surface water and catchment management remain critically important; the numerous large reservoirs that are used to supply Perth, the major population centre of Western Australia, and surrounding regional communities with water receive catchment water but also treated groundwater and desalinated water that is pumped into them as holding capacity. Furthermore, the water that reaches the coastal plains is used directly for irrigation and ultimately ends up in wetlands and estuaries that are of high ecological, social and economic value (Hale and Butcher 2007). Erosion in the escarpments may affect all these values through increased total suspended sediment (TSS) and turbidity but also increased concentrations of dissolved constituents (Armstrong et al. 2005) associated with post-wildfire erosion. On the Swan Coastal Plain, fire-impacted organic-rich wetland sediments have been shown to generate a plume of acidic metal-rich shallow groundwater that can also have downstream consequences for users and receiving environments (Blake 2013).

In January 2016, a large and intense wildfire, the Waroona wildfire, burnt 69 165 ha of forested and agricultural landscapes in the Harvey Catchment south of Perth, Western Australia, and destroyed 181 dwellings, 166 of which were in the historic town of Yarloop (Ferguson 2016) (Fig. 1). The wildfire is said to have started from lightning strikes and burnt over the course of several days, mostly during extreme fire weather on the first day when temperatures were 37°C with wind gusts exceeding 50 km h⁻¹ (McCaw et al. 2016). Fanned by north-westerly, then north-easterly winds, the fire initially burnt through jarrah forest (Eucalyptus marginata) in the escarpment before entering the coastal plains and burning through irrigated agricultural land, ultimately reaching wetlands and dune systems on the coast itself. The potential for post-fire impacts on both soil erosion and water quality represents a management issue for water supply reservoirs, as well as coastal dune systems and the Peel–Harvey Estuary that together make up a significant Ramsar-listed suite of wetlands (Hale and Butcher 2007).

Fig. 1.  Map of the Swan Coastal Plain and Darling Ranges region showing the dune and lakes system to the west, the plains themselves with their intensive network of major drains (lower catchment region), the escarpment (centre and east), and the reservoirs of the upper catchments (upper catchment region). For the purposes of this study, the escarpment forms the boundary between the lower and upper catchment regions of the Harvey River basin. Insert (a) shows a regional perspective of the study site; the extent of the Waroona 2016 fires and the burn severity distribution and place names are also shown. Insert (b) shows the location of the historical water quality (WQ) sampling sites (solid black stars; 1 = site no. 6131335, 2 = site no. 613014, 3 = site no. 613027, 4 = site no. 613052, 5 = site no. 613053, 6 = site no. 613031). Insert (c) shows the distribution of the erosion field sampling sites (solid black squares).

Models for predicting soil erosion are an essential decision-making tool for the rapid allocation of erosion mitigation measures following fire. Although empirical soil erosion models such as the Universal Soil Loss Equation (USLE) and the revised version, RUSLE, were initially developed for agricultural applications, over time, their use has been adapted to forested catchments (e.g. Rulli et al. 2013; Karamesouti et al. 2016). The advantage of RUSLE over other empirical, semi-empirical and physical models is its simplicity, as it requires less field data than other models. When implemented within geographic information systems (GIS) and remote sensing (RS) technologies, RUSLE provides rapid, albeit coarse, post-fire estimations of annual soil erosion. RUSLE has been found to be suitable as a rapid measure of determining high-erosion-risk areas and for determining the effectiveness of mitigation (Vieira et al. 2018).

Set against this background and the likelihood of similar events occurring again in the future, this study explores the risk of enhanced erosion and reduced water quality caused by wildfire in the Darling Ranges of Western Australia in the Harvey River catchment. The specific aims in the present study were to first parameterise the RUSLE (Renard et al. 1991) and determine erosion risks across contrasting landforms (escarpment, dunes and coastal plains) and the variable fire severity in the Harvey catchment. The second aim was to determine whether fire impacts could be detected at pre-existing river monitoring stations where water quality parameters were being measured at regular intervals before and after the wildfire, and to see whether this could be related to the RUSLE outputs.

Situated in south-western Australia (see Fig. 1a), the Peel–Harvey Catchment encompasses all the lands that drain into the Peel–Harvey Estuary. This includes the catchments of the Serpentine, Murray and Harvey rivers and covers an area of ~11 300 km², incorporating two distinct geological provinces across which the fire burnt, namely the Archaean Yilgarn Block of the upper catchment to the east, and the lower catchment situated on the Perth Basin to the west. The upper catchment comprises two geomorphic units, the Darling Plateau and the Darling Scarp. The lower catchment has three geomorphic units, the Swan Coastal Plain, Pinjarra Plain and Ridge Hill Shelf.

The Darling Scarp runs approximately parallel to the coast. Inland, the upper ‘slopes’ of the Darling Plateau are gradual inclines of no more than 10%, but on the escarpment itself, there are steep (up to 30%), but relatively short hillslopes. The elevation of the Scarp is ~100 m above mean sea level. The upper catchment largely comprises Archaean granite–gneiss rocks with intrusions of dolerite, typically overlain with laterite (Playford et al. 1976). The Ridge Hill Shelf, situated at the base of the Scarp, consists of slopes formed by colluvium and has a surface layer of ferruginous, cemented laterite (Davidson 1995). The Pinjarra Plain is primarily of alluvial origin comprising a series of clay overlain with leached sands (Davidson 1995). From this plain to the coast, the Swan Coastal Plain is composed of a series of aeolian deposits of dune systems (Quindalup, Spearwood and Bassendean Dunes), arranged approximately parallel to the coast and decreasing in geological age (late Pleistocene to Holocene) east to west where they continue to form. The dune systems consist of coloured sands, with leaching representative of their geological age, ranging from pale-grey (highly leached) to yellow–brown (slightly leached) to white siliceous (Quindalup) closer to the coast (McArthur and Bettenay 1974).

The Darling Scarp is predominantly remnant jarrah (Eucaluptus marginata) and marri (Corymbia calophylla) open forests. Wetter areas contain bullich (Eucalyptus megacarpa) and blackbutt (Eucalyptus patens) stands. The plateau is mostly state forest that has been subject to periodic timber harvesting; bauxite mining is also undertaken in the region, leading to fragmentation of the forest ecosystem. In the area of the wildfire, the Scarp is a mixture of state forest (jarrah forest, wandoo (Eucalyptus wandoo) woodland and shrubland), water catchments for reservoirs, and freehold land that has largely been cleared for cattle grazing and town sites. Substantial quantities of annual grasses, both native and introduced species, are grown on the Scarp and lower catchment regions.

The Swan Coastal Plain is characterised by low-relief topography. As largely freehold land, it has been extensively cleared for a range of agriculture practices including beef production, dryland and irrigated cropping, and horticulture. Original vegetation included woodlands of Eucalyptus, Banksia and Allocasuarina with Melaleuca predominant in wetlands and swamps. Much of the remaining vegetation is highly fragmented, degraded and subject to invasive annual grasses.

The fire history of the upper catchment is reasonably well known (Peace et al. 2017). Remnant forest areas have not been subject to fire for several decades, with fuel age dating back to the last time broad-scale aerial fuel reduction burning was undertaken in the late 1980s and early 1990s. The exception was a small fire (111 km²) in 2006 that affected several rehabilitated mine sites (McCaw et al. 2016). The burn history of the lower catchment is less well known. Bushland to the east and south of the Waroona town site was burnt by an unplanned fire in January 2015 (McCaw et al. 2016). A small bush remnant in the Yarloop town site was subjected to a prescribed burn in May 2015. The bushland areas surrounding State-owned pine plantations, adjoining the Forrest Highway, have been subjected to prescribed burning. Yalgorup National Park had not experienced a fire for more than 20 years (McCaw et al. 2016).

The region has a mediterranean-type climate with hot dry summers (December to March) and cool wet winters (June–August) (see Fig. 2). The region has a long-term average annual rainfall of ~700 mm near the coast, increasing to 1000 mm on and about the Scarp and then decreasing to ~400 mm at the eastern boundary of the catchment. The winter rainfall provides flow for the three major river systems Serpentine, Murray and Harvey and their associated tributaries. The Serpentine and Murray flow to the Peel inlet and the Harvey River flows into the Harvey Estuary at its southern end (Kelsey et al. 2011). These rivers and their tributaries originate on the Darling Plateau and are well-defined watercourses that traverse the plateau and dissect the Scarp. As these watercourses reach the Swan Coastal Plain, they intersect with numerous wetlands and floodplains and are less well defined (Hale and Butcher 2007).

Fig. 2.  2016 monthly, long-term (1935–2016) mean monthly and mid-term (1990–2016) mean monthly rainfall for the Bureau of Meteorology climatological station in Waroona. Mid-term rainfall data are typical of those used in the RUSLE model. Monthly streamflow in the 1-year post-fire period in relation to long-term (1975–2016) and pre-fire (2011–15) annual mean monthly stream flows.

The application of RUSLE in this study was based on the premise of providing a rapid post-fire estimation of erosion risk areas within the burnt catchment. The objectives were to identify spatial variability in erosion risk and to determine how landscape attributes (e.g. slope, rainfall erosivity) and fire severity contribute to that variability. There are insufficient data from SW Australia to parameterise a physical-based model. To this end, the availability of data dictated model selection and parameterisation.

Burn severity was estimated using the delta Normalized Burn Ratio (dNBR). This is a multitemporal change detection method differencing multispectral single-band indices. Landsat 8 Operational Land Imager (OLI)-derived NBR data calibrated for surface reflectance (30-m resolution) were obtained courtesy of the United States Geological Survey (USGS) Earth Resources and Science Centre (EROS). The NBR indicates which areas of a landscape were burnt, as well as the severity of the burn, using the equation NBR = (NIR – SWIR)/(NIR + SWIR) where NIR is the near-infrared band and SWIR is the short-wave infrared band (Key and Benson 2006). The difference between the pre-fire (15 December 2015) and post-fire (17 February 2016) NBR scenes was calculated to determine the severity of the burn relative to the pre-existing condition of the landscape (Key and Benson 2006). Each pixel was assigned to a burn severity category based on those proposed by Key and Benson (2006). Field verification of burn severity classifications occurred in March 2016. Ten sites per severity classification were identified and field measurements of fire severity assessed using indicators derived from Parson et al. (2010) and Keeley (2009).

This section describes an approach to modelling post-fire erosion risk using the RUSLE, a soil-loss model used for determining hillslope erosion rates. Rulli et al. (2013) produced a version of this model based on the original one developed by Wischmeier and Smith (1965) and revised by Renard et al. (1997) for determining post-fire soil loss from a defined catchment in a mediterranean climate. RUSLE has also been used in other studies from burnt areas (Fernández et al. 2010; Karamesouti et al. 2016; Sankey et al. 2017).

Soil loss (A), in tonnes per hectare per year (in t ha⁻¹ year⁻¹), in the RUSLE is calculated from rainfall erosivity (R), soil erodibility (K), slope and drainage area (LS) and land cover (C).

The equation was implemented at a spatial resolution of 30 m. An overview of the data input for RUSLE can be seen in Table 1.

Table 1.  List of RUSLE factors and subfactors with columns describing the methods, data sources, details of data (equation, source, date, resolution, etc.) and references related to each factor
DAFWA, formerly Department of Agriculture and Food Western Australia, now Department of Primary Industries and Regional Development; n/a, not applicable

Rainfall erosivity is the erosive force of rainfall. It is defined as the intensity and amount of rainfall experienced in a particular region. Erosivity is defined as the mean annual sum of storm erosion index values (EI₃₀), which is a measure of the storm kinetic energy and the maximum 30-min rainfall intensity. Accordingly, such rainfall measurements must be recorded at temporal intervals of less than 30 min.

Daily rainfall data were obtained for 23 Bureau of Meteorology (BoM) (2018) and 10 Department of Primary Industries and Regional Development (2018) monitoring sites within or surrounding the burnt area, with rainfall data spanning the years 1990–2016. EI₃₀ values for each month (j), were estimated following the method proposed by Lu and Yu (2002) using the equation:

where R_d denotes the daily rainfall amount, R₀ the threshold rainfall amount to generate runoff, and n is the number of days with rainfall amount in excess of R₀ in the month, and α, β, η and ω are model parameters. The sinusoidal function with a fundamental frequency F = 1/12 is used to describe the seasonal variation of the coefficient, that is, seasonal variation of rainfall erosivity for a given amount of daily rainfall. The parameter ω is assigned a value of π/6, indicating that for a given amount of daily rainfall, the corresponding rainfall intensity is highest in January, when the temperature is the highest for most regions of the (Australian) continent (Lu and Yu 2002). It is noted that, in the region, most rainfall occurs during winter but these rains have a reduced capacity to initiate significant erosion compared with intense summer rains. During winter, rainfall, which is more frequent but of lower intensity, has time to infiltrate the soil profile, thus reducing the capacity to generate soil runoff. As recommended by Lu and Yu (2002), R₀ was set to 0 mm, therefore, parameter β = 1.49 and parameter η = 0.29, and parameter α was calculated as follows:

where M_R corresponds to mean annual rainfall and Ψ corresponds to mean summer rainfall (November to April).

To account for natural climatic variation, the calculation of both EI₃₀ and R-factor values requires continuous (uninterrupted) daily rainfall data over many years (20 years according to Wischmeier and Smith 1978), which was not the case with the data obtained for some sites. Where rainfall measures were reported following accumulation over 2–5 days, rainfall was split evenly between the days. For BoM sites with missing daily rainfall data and where these sites were situated within close proximity to Department of Agriculture and Food, Western Australia (DAFWA) sites, DAFWA daily rainfall measures were used where there were missing BOM data for those years. For years in which EI₃₀ values were calculated for 9–11 months, values for the missing months were interpolated via tension spline (weight = 5, number of points = 5; negative values set to 0) (Mitášová and Mitáš 1993). Lastly, the historic R-factor surface was generated by taking the mean of the yearly R-factor values calculated at each respective site with >15 years of data (n = 16), and interpolating the values across the study area via tension spline (weight = 5, number of points = 5). Owing to a flooding event in December 2012 in which >100 mm of rainfall was recorded at Harvey resulting in extreme EI₃₀ and R-factor values, an additional R-factor layer was generated excluding 2012 from each site’s historic mean calculation. Rainfall erosivity (R-factor) ranged between 500 and 1188, with the histogram showing a strong binomial distribution that matches the orographic features of the study area; lower values were experienced on the Swan Coastal Plain, and higher values on the escarpment. Generally, rainfall increased from west to east (Fig. 3).

Fig. 3.  RUSLE factors including (a) R-factor: rainfall erosivity, and (b) histogram of R-factor value distribution; (c) K-factor: soil erodibility, and (d) histogram of K-factor value distribution; (e) LS-factor: hillslope length and variation in gradient, and (f) histogram of LS-factor value distribution.

Soil erodibility (K-factor) was determined for unburnt soils. Although K-factor parameters may be affected by fire, this was not taken into account in applying the RUSLE in this study, as we were interested in comparative rather than absolute erosion values. The K-factor was calculated using the equation proposed by Auerswald et al. (2014), which expands Wischmeier and Smith’s K-factor equation (Wischmeier and Smith 1978) adjusting for extreme values in silt content, soil erodibility, organic matter and rock fragments. The parameters required for the equation were mass fraction as a percentage of

1. very fine sand and silt,

2. clay, and

3. organic matter in the fine earth fraction;

4. soil structure index,

5. permeability index, and

6. fraction of the soil surface covered with rock fragments.

Parameters (b), (c) and (f) were derived from the relevant National Soil and Landscape Grid National Soil Attribute Maps (NSAMs) (~90-m spatial resolution) (Viscarra Rossel et al. 2014). The NSAMs use existing digital soil maps to calculate a variance-weighted mean for each pixel. As the K-factor only relates to the top 15 cm of the soil profile (Pringle et al. 2013), the 0–5- and 5–15-cm NSAMs attributes were used in calculations. As the NSAMs contain organic carbon values and not organic matter (parameter (c)), each pixel with an organic carbon value within the K-factor layer was multiplied by 1.72 as suggested by Wischmeier et al. (1971).

Parameter (e) was derived from soil subsystem data supplied by DAFWA using the attribute that describes the capacity of the soil to transmit water based on the least permeable layer in the top 50 cm of the soil profile. Permeability index values were assigned to each DAFWA permeability class based on those listed by Landcom (2004). Each polygon was then assigned a permeability index value weighted by the proportion of each soil type present in that polygon.

As the NSAMs contain information for sand and not fine sand (parameter (a)) as required in the RUSLE, parameter (a) was derived from both the NSAMs and DAFWA soil data. For each soil type within each DAFWA layer, the soil code was matched with those listed by Schoknecht and Pathan (2013). Where soil types were reported by Schoknecht and Pathan (2013) as consisting of fine sands (with or without other textures), it was assumed that all of the sand within these soil types could be classed as fine sand. A new percentage fine sand layer was generated by multiplying the percentage of fine sand soil types within each (rasterised) polygon of the DAFWA soil data by the NSAM percentage of sand for each pixel.

As soil structure data (parameter (d)), as used by Lu et al. (2003), were not readily available, a uniform value of 3 was applied across the study area. Lu et al. (2003) show that the assigned class is representative of major soils having medium to coarse aggregates with diameters between 2 and 10 mm.

Soil erodibility (K-factor) values ranged between 0.01 and 0.79, and the histogram of data was also binomial; for this factor, the lower values, constituting most of the data points, were from the escarpment and the irrigation district beyond its footslopes, whereas the higher values were from the sands of the Swan Coastal Plain, particularly the dune areas closer to the coast (Fig. 3).

The length–slope factor (LS-factor) is a unitless measure of the contribution of slope length and steepness of angle for each cell. The LS-factor was calculated based on the method by Pelton et al. (2014) using the equation derived by Moore and Burch (1986) as follows.

A_s is the area of the cell with a width of 25 m; 22.13 is the standardised length [L] of a RUSLE cell; m and n are adjustable cell coefficients relating to soil variants and erosion factors respectively.

Length and slope (LS-factor) values ranged from 0 to 23.5. The overwhelming majority of values were low, exemplifying the generally flat terrain of the Swan Coastal Plain. The highest values (also seen as a small peak in the histogram) were seen on the escarpment itself, demarcating the western edge. The dune system at the coast appears as three faint bands parallel to the coast (Fig. 3).

The land cover factor (C-factor) layer was derived from the DAFWA land-use dataset, with estimated values assigned to each land use based on values developed by Panagos et al. (2015). Australian Land Use and Management (ALUM) classifications were matched with the most similar land uses listed by Panagos et al. (2015) based on description and desktop-based visual inspection of high-resolution RS imagery, with the middle value of the C-factor range from Panagos et al. (2015) applied to each respective ALUM class. The C-factor is representative of the land cover, including native vegetation, plantations and crops and their management and how this land cover might vary in erosion rates when compared with bare, tilled soil (Kinnell 2010). Values are assigned to land cover types and range between 0 and 1, with 1 representing high erosion potential associated with bare, tilled soil (Renard et al. 1997).

As no land uses listed by Panagos et al. (2015) were close matches with the ALUM subcategories ‘Production forestry’ and ‘Plantation forestry’, the C-factor values for these land uses were estimated based on those listed by Bonilla et al. (2010) and Hartcher and Post (2005).

The cover management factor was modified to include the effects of the burn. Change in burn severity (dNBR) ranged between –2 and 2. More severe fire is shown on the western area of the escarpment, and another patch among the dune systems closer to the coast. Between and around these areas were less severe burns. In fact, most (97%) dNBR values were in the range 0–0.5, indicating low-severity burns; values less than 0.1 typically indicate unburnt conditions. The histogram shows a skewed distribution with a low burn severity median value for the majority of the burnt area (Fig. 4).

Fig. 4.  RUSLE inputs for (a) pre-fire C-factor (land cover) distribution, and (b) histogram showing frequency of pre-fire C-factor values; (c) image of dNBR (delta Normalized Burn Ratio) results, and (d) histogram showing frequency of dNBR values; (e) image of post-fire C-factor results, and (f) histogram showing the frequency of post-fire C-factor values.

The change in burn severity, or change in the NBR, was used as an input to the C-factor. Burn severity was factored into the C-factor layer based on values used by Terranova et al. (2009) for different burn severity classes. The burn severity layer was reclassified, in which pixels with a high burn severity were assigned a value of 0.2, moderate–high burn severity a value of 0.1, moderate–low burn severity a value of 0.05 and low burn severity a value of 0.01. This layer was merged with the C-factor layer to produce the final post-fire C-factor layer. The maximum pixel value from either of the merged layers was retained in the final post-fire C-factor layer. Cover (C-factor) values, after taking into account fire severity, ranged between 0.01 and 0.3, mainly reflecting where native vegetation has been cleared or converted to land for pasture and grazing on the Swan Coastal Plain and in the irrigation areas. Peaks of frequencies occurred around the low- and medium-cover management values; smaller peaks occurred at low–medium and high values, with the latter reflecting the vegetated areas on the western escarpment and in the dune system at the coast (Fig. 4).

The outputs from RUSLE were examined to determine if the predicted erosion rates were consistent with observations of post-fire erosion at the hillslopes scale. There are no data available for the Harvey catchment. However, studies of background erosion rates and fire impacts elsewhere in the region can provide some insights into the reliability of the predictions. A field survey was carried out to determine the degree with which RUSLE represents differences in erosion potential quantified from a rapid visual field assessment technique developed by Morris et al. (2014). Fifteen sample sites were selected to cover a range of erosion potentials as predicted by RUSLE and based on accessibility (see Fig. 1c). At each sample site, 10 locations were sampled, and classified into one of five erosion categories ranging from 0 (least erosion) to 5 (most significant erosion). The coordinates for the initial location at each site were chosen randomly, and then nine further sites were selected at least 30 m apart (as the spatial resolution of the RUSLE is 30 m), aligned parallel to the hillslope. Locations situated below roads were not chosen because of their interference on surface runoff processes. The median of the 10 locations was then correlated with the mean RUSLE output (t ha⁻¹ year⁻¹) for the sampled area.

The relationships among RUSLE, slope and burn severity were examined to determine how these variables contribute to variation in erosion risk. Fire severity was classified into four categories (unburnt: typically –99 to 100), 2 (low severity: 101–350), 3 (moderate severity: 351–500) and 4 (high severity: >500). Slope was subdivided into five categories (all values in degrees): 0–5, 6–10, 11–15, 16–25, >25. The mean of RUSLE outputs was calculated for all combinations of these categories. The effect of slope and dNBR on RUSLE predictions was quantified with an ANOVA in Matlab. The residuals of the RUSLE were normally distributed when using the untransformed data in the ANOVA, so no transformation was necessary.

Impacts of wildfire on water quality were assessed using data collected from six monitoring stations within the Department of Water and Environmental Regulation’s extensive monitoring network (Department of Water and Environmental Regulation 2018). Samples were collected following standard protocols outlined in Department of Water (2009a, 2009b). The sites were selected because they met the criteria of having at least 5 years of continuous pre-fire water quality data and continuous post-fire data for 2016. Table 2 presents the characteristics of the six sites and their catchments, showing that four of the sites had all or part of their catchments burnt. Three of these were nested within the Harvey River Catchment, which was monitored near the outlet. Two control catchments were entirely outside the burn area and not hydrologically connected to surface waters directly impacted by the wildfire.

Table 2.  Characteristics of water quality monitoring sites used for pre-fire v. post-fire comparisons
The table shows the design criterion for each site, its latitude and longitude, along with characteristics of the subcatchment that is drained to the site. The total area is given in kilometres squared; and area and land-use categories drawn from Department of Water Peel Harvey Catchment Nutrient Report 2015 Update 2017 data sheets: Conservation and Natural (CN), Plantation (P), Grazing (including beef and dairy cattle, mixed) (G), Built (including housing, industry, roads) (B). Percentage of the subcatchment burnt is calculated from the area of the burn extent that intersects with each individual catchment area

Concentrations of total nitrogen (TN), total organic nitrogen (TON), ammonia/ammonium-N (abbreviated to NH₃), nitrate and nitrite (NO_x), total phosphorus (TP), and soluble reactive phosphorus (SRP), together with physicochemical parameters of temperature, dissolved oxygen, pH and TSSs were analysed at all of these stations (which were always sampled at the same time of day).

These pre-and post-fire data were analysed by comparing means for each parameter at each site using Welch’s t-test. A principal component analysis for each site was performed to determine whether together all the measured parameters were indicating a pre-fire to post-fire shift in multivariate space. Following this, each parameter was analysed separately using multiple linear models, in which each measurement was standardised and its deviation from the site averaged, and plotted over time.

Pre-fire RUSLE output (i.e. dNBR set to 0) can be seen in Figure 5a. Post-fire, the RUSLE produced values that ranged from 0 to 94.5 t ha⁻¹ year⁻¹ (Fig. 5b). Parts of the burn area with higher RUSLE values (>1 t ha⁻¹ year⁻¹) are the ridges of two dunes in the coastal areas, the western portion of the escarpment, and patches of the hinterland forested areas east of the escarpment in the Darling Ranges. When compared with unburnt conditions (Fig. 5a), these impacts of the wildfire on RUSLE in these areas are visible. In general, the most widespread change is a shift in erosion risk from values of the order of 0.01–0.1 t ha⁻¹ year⁻¹ to 0.1–1 t ha⁻¹ year⁻¹. There are localised areas where the erosion risk has shifted two orders of magnitude, from 0.01–0.1 to 1–10 t ha⁻¹ year⁻¹.

Fig. 5.  RUSLE output (t ha⁻¹ year⁻¹) for the Peel–Harvey catchment (a) without, and (b) with fire. White areas indicate ‘no data’.

Data from the erosion field survey were summarised for each site by the median erosion category among the sampled locations. These most typical erosion conditions of the site could then be related to the corresponding RUSLE prediction to evaluate model performance in terms of providing an indication of erosion susceptibility (Fig. 6). The results show an exponential relationship between erosion categories (from Morris et al. 2014) and the RUSLE predictions. Erosion categories ranging from 0 to 2.5 correspond with erosion rates ranging from ~0.05 to 1 t ha⁻¹ year⁻¹. The exponential relationship suggests that erosion magnitude scales non-linearly with the erosion categories from Morris et al. (2014). In fact, for the equation relating the RUSLE prediction and erosion category (E_cat) () (Fig. 6), each erosion category corresponds to an order of magnitude increase in RUSLE estimates. This result is reasonable because the categories discriminate between erosion types that represent threshold shifts, which are associated with large changes in erosion rates (e.g. rainsplash to interill, and interill to rill erosion).

Fig. 6.  Predicted RUSLE output (t ha⁻¹ year⁻¹) v. the median of erosion category (0–5) measured from 10 samples in 30 × 30-m plots centred on the location of cells in the grid of RUSLE outputs. The transition in erosion processes with erosion categories from Morris et al. (2014) are listed along the top x-axis.

When the variation in post-fire erosion risk, shown as a function of slope and fire severity, is split according to the strongly bimodal distributions of soil erodibility (K) and rain erosivity (R) factors (Fig. 7), three inferences can be made. First, in areas with steep slopes, there is a strong non-linear effect on erosion risk with fire severity, where predicted erosion is higher. Further analysis, using analysis of variance, confirms slope as the first-order control on erosion risk (Table 3). Second, the soil erodibility K in the sand dunes is substantial, yielding an order of magnitude larger erosion risks (Fig. 7a v. Fig. 7b). Third, the effect of differences in rainfall erosivity (R) between the escarpment and flood plain is almost negligible compared with the influence of slope.

Fig. 7.  Variation in post-fire erosion risk shown as a function of slope and fire severity and split according to soil erodibility (K) and rain erosivity (R) factors. (a) Variation in RUSLE output with slope and dNBR (delta Normalized Burn Ratio) for K ≤ 0.4. (b) Variation in RUSLE output with slope and dNBR for K > 0.4. (c) Variation in RUSLE output with slope and dNBR for R ≤ 825. (d) Variation in RUSLE output with slope and dNBR for R > 825.

Table 3.  Result from analysis of variance with fire severity (dNBR, delta Normalized Burn Ratio) and slope as independent variables and RUSLE output as the dependent variable
sq., squares

Using the subcatchments where water quality data were collected pre-and post fire, data revealed differences in mean slope, mean fire severity and mean RUSLE outputs calculated with and without fire (Table 4). Reference subcatchments were confirmed (showing no differences with and without fire); each of the four test subcatchments had elevated outputs post fire, where the magnitude of differences for RUSLE outputs with and without the fire (dA) was mostly proportional to mean fire severity (Table 4).

Table 4.  Characteristics of sites used for pre-fire v. post-fire comparisons, and post-fire measurements
The table shows the design criterion for each site, mean slope derived from RUSLE input data, mean dNBR (delta Normalized Burn Ratio) (fire severity). Mean RUSLE output (t ha⁻¹ year⁻¹) without and with the fire is produced by manipulating the cover parameter in the equation. dA is the difference pre- and post-fire, derived as the difference between the two

A comparison of nutrient and suspended solids concentrations between pre-fire and post-fire monitoring at each of the six sites demonstrates that the differences between pre-fire and post-fire concentrations are minimal for nearly all parameters at all sites (Fig. 8); only a few exceptions to this general rule are apparent. For both the reference sites, both sites entirely within the burnt area and one of the sites downstream of the burnt area, some nitrogen variables were different at the 0.05 level of significance (Fig. 8). The similarities in nutrient and suspended solid concentrations, even for these nitrogen variables, given that they were recorded for all sites irrespective of where they were located in relation to the fire, suggest that the monthly water quality monitoring has not detected an effect of the fire.

Fig. 8.  Box plots showing means (triangles), medians (horizontal lines), interquartile range (excluding outliers; whiskers), and outliers (defined as all points that are three scaled median absolute deviations (MADs) away from the median; circles) for pre- and post-fire water quality measurements total nitrogen (Tot N), total organic nitrogen (Org N), ammonia/ammonium (abbreviated to NH₃), nitrate plus nitrite (NO_x), total phosphorus (Tot P), soluble reactive phosphorus (SRP) and total suspended solids (TSSs) for six Department of Water and Environmental Regulation Water Information Network (WIN) sampling sites (the locations of each sampling site are shown in Fig. 1b and Table 1). Data are for the months of January to August from 2010 to 2015 inclusive. Post-fire averages were calculated for the corresponding months in 2016. WIN sites 613027 and 613031 are Reference – outside the burn area with none of the catchment burnt; sites 613052 and 6131335 are outside the burn area but downstream in a catchment that is partly burnt; 613014 and 613053 are within the burnt area with much of their catchment burnt. Asterisks immediately beneath a parameter in each graph denote significant differences between pre-fire and post-fire values (using a t-test).

These interpretations are supported by the principal component analysis (PCA): in multivariate terms, and including the physicochemical parameters (for example Fig. 9), the post-fire values were not distinct from the pre-fire values. Looking at each parameter separately shows that the variability encountered post fire cannot be readily distinguished from the variability experienced pre-fire.

Fig. 9.  Principal component analysis (PCA) plots for pre- and post-fire water quality results for nine water quality and physicochemical parameters (dissolved oxygen [DO], temperature, pH, total suspended solids, soluble reactive phosphorus, total phosphorus, nitrite plus nitrate, ammonia/ammonium and total nitrogen records) from January 2011 to December 2016 for six Department of Water and Environmental Regulation Water Information Network (WIN) sampling sites (the locations of each sampling site are shown in Fig. 1b and Table 1). WIN sites (a) 613027 and (b) 613031 are Reference – outside the burn area with none of the catchment burnt; sites (c) 613052 and (d) 6131335 are outside the burn area but downstream in a catchment that is partly burnt; (e) 613014 and (f) 613053 are within the burn area with much of their catchment burnt.

Both measured and modelled erosion rates in undisturbed forests in Australia are typically in the range 0.01–0.1 t ha⁻¹ year⁻¹ (Lu et al. 2003; Lane et al. 2006; Tomkins et al. 2007; Hancock et al. 2017) and RUSLE has been shown in the past to reproduce these values with reasonable accuracy (Lu et al. 2003). In our implementation of the RUSLE, the erosion rate in the forested headwaters and vegetated coastal dunes is typically <0.1 t ha⁻¹ year⁻¹, indicating that the model is performing well. The fire effect on erosion is represented through adjustments to C, which is a factor representing vegetation cover. Although there is a distinct lack of region-specific C-factor data for this study, the use of values for vegetation in similar climates provided comparable C results with other similar studies. Typically, in the forested headwaters, the burn effect resulted in C changing from ~0.028 (in the absence of fire) to 0.1 (moderate severity) or 0.2 (high severity). This change in C has meant that modelled erosion rates increase by a factor of 3.5 or 7 for moderate- and high-severity fire respectively. In the forested headwaters and coastal dunes, this change in erosion potential due to wildfire has resulted in modelled erosion rates of 0.1–5 t ha⁻¹ year⁻¹ for hillslopes typical of the area. On the coastal plains, the LS factor was very low (owing to flat terrain) and hillslope erosion was therefore negligible irrespective of fire, with the exception of the dune system near the coast, because the dunes have moderate slopes with sandy, poorly structured soils that are prone to erosion if denuded of vegetation by fire.

The modelled post-fire erosion rate on hillslopes is relatively low compared with the modelled post-fire erosion rate using RUSLE in other studies where values for the post-fire erosion rate were modelled for moderate slopes. Karamesouti et al. (2016), for instance, reported the post-fire erosion on hillslopes with a gradient of 8–15% was 80.4 t ha⁻¹ year⁻¹, which is an order of magnitude higher than modelled post-fire erosion for similar slopes in the present study. Terranova et al. (2009) found that 30–40% of an area had modelled post-fire erosion rates >7.5 t ha⁻¹ year⁻¹ after fire. In our study, the area with modelled erosion rates >7.5 t ha⁻¹ year⁻¹ was very small, partly owing to the large proportion of coastal plain but also owing to relatively low modelled erosion rates on hillslopes above the plain. In contrast, the range of modelled post-fire erosion rates and the modelled impact of the fire are very similar to the values reported by Rulli et al. (2013).

Erosion was not monitored so a quantitative test of the accuracy of predictions made with RUSLE is not possible this way. However, the modelled erosion rates can be compared with the measured erosion rate values from hillslope experiments elsewhere. In small catchments (0.2–0.3 ha) with similar attributes to those in the headwaters of the Harvey River, the erosion rates after the 2009 catastrophic Black Saturday Fire in south-eastern Australia were 7–10 t ha⁻¹ in the first year after the fire (Noske et al. 2016). These values are in the upper range of the typical erosion rate modelled on steep hillslopes in the Waroona burn area. On small plots (8 m²), Blong et al. (1982) measured erosion rates between 2.5 and 8.2 t ha⁻¹ year⁻¹ in the first year after fire in sandstone geology. This rate is in the upper range of our modelled results. The study by Blong et al. (1982), however, concerns erosion on small plots, and scale effects (see for example Amore et al. 2004 and Parsons et al. 2006) mean that their values are therefore likely to be higher than hillslope-scale erosion predicted by RUSLE. The range of modelled values in our study is very similar to those modelled with RUSLE (0.5–5 t ha⁻¹ year⁻¹) and measured on small plots (0.5–3 t ha⁻¹ year⁻¹) in the Warrumbungles by Yang et al. (2018). In summary, it appears that the RUSLE produces reasonable predictions of hillslope erosion under unburnt conditions, but that it may underestimate erosion in areas where slopes are steep and fire severity high. The tendency for RUSLE to underestimate has been documented in other work (Larsen and MacDonald 2007) and may be due to fire-related changes in soil erodibility, K, which are not captured in the model. In a relative sense, the assessment of erosion severity using the method of Morris et al. (2014) indicates that RUSLE is effective at discriminating between areas of different erosion risk. Hillslopes with low erosion predictions showed few signs of erosion, while locations with comparably high RUSLE predictions displayed significant erosion, including sheetwash and rills.

In this study, we implemented RUSLE in a manner that will predict annual erosion and not individual erosion events. By implementing it this way, and including monthly rainfall variation, we represent variation; important elements of the regional rainfall regime that cause high-magnitude erosion events were highlighted. Using this approach, there appears to be a quite modest increase in predicted fire-induced erosion, compared with what has been predicted and/or measured in other Australian locations. However, this requires validation through field trials of post-fire sediment flows. The water quality measurements on the coastal plains collected during a 1-year period following the burn also suggest that the impacts of fire were relatively low. However, given the episodic nature of post-fire erosion, our approach to modelling and measurement may underestimate post-fire impacts. The magnitude of a post-fire erosion/yield event is dependent on the interaction of several factors such as fire severity, vegetation type, soil type and topography. Variability in erosion due to stochasticity in rainfall within the burnt landscape is particularly important (Shakesby and Doerr 2006; Robichaud et al. 2007; Cannon et al. 2010; Moody 2012). Through the rainfall erosivity parameter, RUSLE incorporates the long-term average effect of rainfall events on erosion. However, the model does not predict individual erosion events, where erosion rates are often dictated by thresholds associated with short bursts of intense rainfall. Thus, by using long-term monthly rainfall erosivity, our approach to modelling erosion with RUSLE is likely to underestimate the erosion produced in the case of an intense rainfall event. Recent developments in modelling rainfall erosivity offer opportunities for improved representation of event-based erosion rates (Yang et al. 2018). Additionally, RUSLE, by its very design, does not consider connectivity between eroding hillslopes and downstream waterways. The model indicates the exposure of catchment to erosion, but is limited in its ability to predict impacts beyond the hillslope scale, where hillslope and channel coupling is important (Smith and Dragovich 2008; Smith et al. 2011c). Similarly, the results from the water quality sampling are subject to issues related to the episodic nature of post-fire water impacts. The sampling approach, whereby samples were collected at regular intervals, may have missed important events that occurred between sampling campaigns. The wildfire signals in water quality parameters are often tightly coupled to the runoff-generating mechanism (Murphy et al. 2018). A more targeted event-based water quality sampling regime and hillslope erosion experiment are needed to obtain more robust estimates of fire impact water quality parameters. Despite these limitations, the results from the RUSLE model (Fig. 5) and associated analyses (Figs 6 and 7) provide the first insights into wildfire impacts on erosion and water quality in SW Australia. It is an important contribution that will help in understanding the spatial distribution of erosion risk and help to prioritise mitigation and future research activities.

An important oversight in our understanding of post-fire response was the coastal system, where RUSLE indicates erosion risk is very large. High erosion rates were predicted for the dune ridges and adjacent freshwater wetland systems at the coast. Our original focus was on what we perceived were the important water resources of the escarpment and the Darling Ranges, like the reservoirs, the irrigation infrastructure on the Swan Coastal Plain and the likely impact they were to have on flows and inputs for the Harvey Estuary. The aquatic ecosystems, such as the freshwater wetland systems, close to the coast deserve closer attention, and the RUSLE outputs provide us with the evidence needed to add to our foci for such a study.

The data from water quality monitoring stations point to the fact that monthly snapshot sampling is temporally too coarse to pick up subtle and event-related effects of fires on water quality. Added to this, the locations of the pre-existing sites were distal to the escarpment where fire intensities, slopes and rainfall erosivities were higher, thereby distal to the zone of contaminant production, resulting in a lower likelihood of detecting an effect. Post-fire high-magnitude water quality events tend to be episodic with quite short time scales of minutes to hours (Cannon et al. 2010; Moody 2012), requiring continuous monitoring at strategic locations to aid in their detection. More importantly, RUSLE (empirical models), while proven to accurately predict erosive response in a catchment, may not fully capture hydrological responses. Vieira et al. (2018) have shown that semi-empirical and physical models may be more accurate in predicting water contamination and warrant further investigation in low-relief catchments.

The hydrology of the Swan Coastal Plain adds complexity to sampling for post-fire water quality impacts. The low relief, porous sands, land use and the direct hydrological connection with groundwater mean that trying to discern a post-fire water quality signal is complex and requires a multifaceted approach to sampling methodology. This is particularly the case where a signal may be more subdued and interspersed with water quality signals from agriculture, such as nutrient exports (McKergow et al. 2003) or other disturbance events. Typically, most post-fire erosion modelling is done at the hillslope scale in moderate- to high-topographic-relief catchments. Little, if any work has been done post fire in low-relief landscapes and the interface between hillslope and coastal plain.

Understanding contaminant movement throughout catchments in the south-west of Western Australia must not ignore the numerous contaminant pathways, including overland flow, streamflow infiltration and groundwater interaction. Importantly, and related to the Swan Coastal Plain, we did not investigate the quality of groundwater. Future investigations could compare pre- v. post-fire records for existing regularly monitored groundwater bores in both the lower and upper catchments of the region.

The wildfire in 2016 in the catchment of the Harvey River in south-western Australia raised immediate questions about its impact on water resources. The ability to quantify overland sediment movement throughout a catchment, coupled with sediment and nutrient transport within waterways, will provide a model for identifying critical resources that require protection and for informing fire management practices. The existence of remote sensing data for burn severity, land capability data for modelling soil erosion, and routinely collected water quality monitoring data provided an opportunity to undertake a post-fire investigation of this impact.

Overall, RUSLE modelling predicted patchy erosion, consistent with a moderate hillslope catchment experiencing a mediterranean climate, and yielding movement of moderate sediment loads compared with more mountainous catchments elsewhere in Australia. Nevertheless, hot-spots of erosion were predicted for intersections of steep terrain and high fire severity that were confined to two particular landforms: forested headwaters and coastal dunes.

We anticipated water quality impacts in the upper catchment would reflect the rainfall interaction with the burnt landscape, and show elevated nutrients levels and higher total suspended solids. While water quality data showed conspicuous seasonal patterns, sampling frequency was temporally too coarse to detect the nature, magnitude and frequency of post-fire erosion events. This was particularly so given that the monitoring sites were determined well before the fire, and were distal to the zone of contamination predicted by hot-spots in the RUSLE outputs. We did not anticipate the high erosion in coastal dunes and we were unable to consider associated water quality impacts due to an absence of regularly monitored surface water and groundwater in the vicinity of the hot-spots.

Application of the RUSLE in this study identified deficiencies in available data for the factors used to predict post-fire catchment erosion. In particular there is a lack of post-fire vegetation, soil data, and patterns of fire frequency and the drivers of burn severity that are particular to the south-west region of Western Australia. Added to this is the unique geomorphological setting of the region: the interplay between a topographically subdued and ancient escarpment, a broad coastal plain and active dune systems on the coast itself. Future modelling, addressing these features, particularly given an awareness of climate-related vulnerabilities, is warranted.

The authors declare that they have no conflicts of interest.