Study area
The Sydney Basin bioregion covers approximately 45,000 km²along the eastern seaboard of southeastern Australia (Figure 1). Soils are mainly of low fertility, derived from sandstone and shale parent material (DPI 2017). Elevation ranges from sea level to over 1200 m. Mean annual temperature and rainfall ranges between 11°C–18°C and 600 mm to 1500 mm, respectively, as a function of both altitude and distance from the coast (e.g. stations 068024, 063292, 066214, 067027, 061394, http://www.bom.gov.au/climate/data/, 08/07/2021). DSF dominates ridgetops and WSF dominates gullies, with most forested land occurring within the National Park estate (see Online Resource 1 for details on vegetation types). Typical fire intervals are between 5–20 years in DSF and 20–100 years in WSF (Murphy et al. 2013)and fires are generally of mixed fire severity(Collins et al. 2021). Fires occurring at very short intervals (<5 years) are likely to result in changes to forest structure and composition, when compared to areas burnt at longer intervals (Arno & Allison-Bunnell 2002; Lewis et al. 2012; Cawson et al. 2017). Most sclerophyll forests in the study region have burned 1–3 times since the 1970s, when reliable fire history records began (Bradstock et al. 2009; Price & Bradstock 2010), with a smaller proportion (~10%) burning in excess of three times (Hammill, Tasker & Barker 2013). The occurrence of extreme wildfires has increased across the broader temperate biome in response to warmer drier conditions (Sharples et al. 2016; Abram et al. 2021; Collins et al. 2021), with large wildfires (>100, 000 ha) occurring across the study area in 1993/94, 2001–03, 2006/07, 2013/14 and again in 2019/20 following our study.
Fire history
The study focused on areas of the Sydney basin that were most recently burnt by large fires in October 2013 (Figure 1), to control for the potentially confounding effects of time since fire. The fires occurred in four different sub-regions across the Sydney Basin, with two sub-regions experiencing mild/moderate drought (MD) and two experiencing severe drought (SD; Fig. 1). Fire frequency was calculated for each sub-region as the number of fires that occurred between 1993 and 2013 and was categorized as low (1–2 fires; LF) or high (3 or more fires; HF). Fire history data were obtained from the New South Wales National Parks and Wildlife Service (NPWS 2016).
Drought severity
Drought severity was quantified by using the Standardized Precipitation-Evapotranspiration Index (SPEI, Vicente-Serrano et al. 2010). SPEI is an index of the climatic water balance (i.e. precipitation minus potential evapotranspiration)expressed as the number of standard deviation units from average values calculated over a 30 year period (e.g. 1980–2010). Negative SPEI values indicate increased water deficit relative to long-term conditions, while positive values indicate surplus water availability relative to long-term conditions. We calculated SPEI at a 6-monthly time scale, which is sufficient for detecting drought stress in temperate eucalypt forests (Pook 1986; Pook, Gill & Moore 1997). Slette et al. (2019) suggest SPEI values between -1 and +1 can be treated as falling within the range of normal climatic variability, while values below -1 represent progressively more severe drought conditions; however, values at or below -0.5 have been considered to be representative of drought conditions in temperate Australia (Ma et al. 2015). We partitioned sub-regions in our study into either mild/moderate drought (MD; SPEI = 0 to -1.4) or severe drought (SD; SPEI = < -1.4). This threshold was chosen as it divided the study area into two approximately equal-sized and climatically coherent regions, each contain substantial DSF and WSF populations. For example values below -1.4 were restricted to typically drier, warmer areas at lower elevation in the north-west of the study region while values above -1.4 occurred along the coastal fringe and in cooler areas with higher elevation (Fig. 1).
In the six months preceding the 2013/14 fire season, drought severity varied considerably within the study region, with fires in the northern sub-regions burning under severe drought and fires in the southern sub-regions burning under mild/moderate drought (Fig. 1).Most areas returned to low drought/normal conditions in the six months following October 2013. For each sub-region, SPEI was calculated using spatially gridded climate data at 0.05° x 0.05° resolution for the 6-month period prior to and after the 2013 fire. Field sites were only placed in areas where post-fire SPEI had returned to normal/near normal, to avoid the confounding effects of post-fire drought on juvenile mortality and recruitment. Thus, there was substantial variability in pre-fire drought severity and minimal variation in post-fire drought conditions across the study regions. Climatic data used to calculate SPEI was obtained from the SILO database (SILO 2019). SPEI was calculated using the ‘SPEI’ package in R (Vicente-Serrano et al. 2010).
Study design
The study design incorporated drought severity (MD; SD), fire frequency (LF; HF) and vegetation type (DSF; WSF) in a fully factorial manner, with 14 replicate sites per treatment combination (n = 112). The 112 sites were evenly distributed across the four sub-regions (28 sites per sub-region) to obtain sufficient spatial variability in drought severity (Fig. 1). Sub-regions occupied narrow bounds of mean annual temperature and rainfall (±2°C and 200 mm across sites within each sub-region) to control for climatic variability. All DSF sites were last burnt at moderate-high severity, with a high amount of scorching and consumption of canopy foliage (canopy 70–100% burnt; severity classes 3–5, Hammill & Bradstock 2006), whereas all WSF sites were burnt at low-moderate severity, with a mix of unburnt and scorched canopy foliage (canopy <70% burnt; severity classes 1–2, Hammill & Bradstock 2006). Because the topography of the study region limits the prevalence of high severity fires in gullies (Bradstock et al. 2010), fire severity could not be matched between vegetation types. Consequently, we have contrasted the ‘common’ fire severity patterns between vegetation types: high severity in DSF and low/moderate severity in WSF (Bradstock et al. 2010).
A 50 m x 5 m plot was established at each site following established forest measurement protocols (McElhinny et al. 2006; McElhinny et al. 2005). DSF plots were confined to the top of ridges along contours, whereas WSF plots were confined to gully bottoms or lower slopes, adjacent to creeks along contours (see Online Resource 1 for examples of typical sites). Plot aspect varied between sites to minimize aspect bias. Plots were selected randomly within a few kilometers of access roads and within the treatment levels identified in a GIS. Plots were placed at least 50 m from roads and trails to avoid edge effects and at least 300–500 m apart to reduce the effect of spatial autocorrelation. Plots were surveyed between February 2018 and July 2018.
Field methods
All juvenile trees between 2.5 and 10 cm diameter at breast height over bark (DBH) were identified and individually measured within the 50 m x 5 m plot. Species were identified using the keys provided by Klaphake (2012) and Brooker and Kleinig (1999). Juvenile stems that arose from dead stems >10 cm DBH were included, while juvenile stems that arose from live stems >10 cm DBH were not included, i.e. when trees were multi-stemmed, the largest living stem was used to determine that maximum DBH. When stems were closely-spaced, a 1 m long steel rod (4 mm diameter) was used to probe between stems to determine whether they were connected by a sub-surface lignotuber. To determine whether a juvenile was a new post-fire seedling or a surviving resprout (Fig. 2), the base of the stem was excavated of soil and manually checked for lignotuber presence. The DBH of each juvenile stem was measured over bark at 1.3 m above the ground on the uphill side of the tree.
Mortality was defined as a dead standing stem or downed stem/associated stump representing an individual that had died due to the most recent fire (Fig. 2). Stumps and downed stems were only included if they were determined to be a product of the most recent fire based on criteria described by Gordon et al. (2018) and Roxburgh et al. (2006); and had most likely been felled via fire scar formation and collapse, evidenced by a fire scar/break point. Mortality of juveniles smaller than 2.5 cm DBH was unable to be determined, as charred stems in this size-class looked similar to other non-eucalypt plant genera.
Data analysis
We fitted Bayesian regression models to analyse the influence of fire frequency and drought severity on each of the following response variables: the probability of juvenile mortality; the number of post-fire seedlings; the post-fire replacement balance (number of seedlings minus the number of dead juveniles); and the total post-fire juvenile abundance (number seedlings plus the number of surviving resprouts). For all models, the single predictor was a four-level categorical variable giving the combination of fire frequency (low versus high) and drought severity (mild/moderate versus severe).
Juvenile mortality was modelled as a Bernoulli process via a logit-link function. A weighting term was included to account for different plot sizes between standing stems and downed stems. The number of post-fire seedlings was modelled as following a negative binomial distribution parameterized in terms of mean and dispersion. A hierarchical model was fitted in which, for each combination of drought and fire classes, the priors for distribution parameters were informed by overall priors. We chose this model structure to ensure more reliable inferences given the relatively small number of sites within each combination of fire and drought classes, and the occurrence of several large outlier values in the data. Since the fitted negative binomial distributions could be strongly right-tailed, we monitored posterior median values rather than posterior means. The model for post-fire replacement balance followed a similar hierarchical structure as that for the number of post-fire seedlings. However, since the data included negative values, it was treated as continuous and modelled using a location-scale t-distribution, with the mean and standard deviation parameters specific to each drought and fire combination, and a global shape parameter learned by the model to reduce the influence presence of several large outliers in the data. While this method monitored posterior means rather than posterior medians, it shared the same intent as in the other models, i.e. to obtain an estimation of the central tendency/most likely values.
Models were fitted via Markov Chain Monte Carlo using R version 3.5.0 (R Core Team 2019). The juvenile mortality model was fitted using the brms package (Bürkner 2018). The models of post-fire seedlings and post-fire replacement balance were fitted using the JAGS program (Plummer 2003) via the runjags package (Denwood 2016).
For each model, we sampled four Markov chains, each consisting of at least 5000 model iterations. We assessed model convergence using the diagnostic of Gelman and Rubin (1992) and checked for acceptable levels of serial autocorrelation. Separate Markov Chains for each model were then combined into a matrix of samples from the joint posterior distribution of model parameters, which we subsequently used to derive predictions of probabilities/tree count per site among the treatments (Kruschke 2015; Suzuki 2019). We then calculated the difference between selected comparisons by arithmetically generating a distribution of differences that could be used to inform interpretation of the magnitude of differences among treatment combinations. Hence, these calculations are referred to in the results as ‘median posterior difference’, i.e. the median value of summarised difference calculations. Credible intervals were calculated as highest posterior density intervals (HPDI), in order to display the central 50% of model predictions and lower/upper 95% bounds of model predictions.
Data for all models, with the exception of the juvenile mortality model, were aggregated by site (DSF, n = 56; WSF, n = 56). We modelled DSF and WSF independently due to confounding by fire severity. The data and R scripts used to generate the results are provided online in a data repository (https://github.c-om/erb418/EB.Ch3.scripts) and secondary results summaries can be found in Online Resource 2.