**1. Introduction**

Wildfires have always played a significant role in the evolution of various ecosystems and are the predominant natural disturbance factor in many parts of the world. They significantly influence ecological patterns and processes on a global scale. This includes vegetation distribution and structure, as well as the carbon cycle [1]. While humans and wildfires have always coexisted, changes in wildfire patterns represent an increasing threat to human lives and property. Apart from the direct implications, wildfires have also been found to contribute to the greenhouse effect through CO2 emissions, thus fostering atmospheric changes on a global level [2,3]. Research has shown that forest loss has increased substantially over the past two decades in many parts of the world, and that the underlying dynamics are strongly linked to fire activity [4]. Several studies have discovered changes in the frequency and size of wildfires and also in the length of the fire season, for example, regarding the Canadian boreal forest [5] and the Western United States [6].

In recent years, large wildfires have occurred in regions formerly unaffected by fire, such as the Arctic regions. Some regions regularly affected by fire have experienced unprecedented large-scale fire events, such as the Australian East coast, the Brazilian Amazonas region, or the state of California in the United States. Investigating the question if these recent events are part of a natural oscillation, or must instead be regarded as a result

**Citation:** Nolde, M.; Mueller, N.; Strunz, G.; Riedlinger, T. Assessment of Wildfire Activity Development Trends for Eastern Australia Using Multi-Sensor Earth Observation Data. *Remote Sens.* **2021**, *13*, 4975. https:// doi.org/10.3390/rs13244975

Academic Editor: Paolo Mazzanti

Received: 9 November 2021 Accepted: 3 December 2021 Published: 7 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of a long-term trend, is a crucial task in fire science. Studies usually analyze parameters such as the frequency of occurrence, spatial extent, and burn severity of fire incidents to derive meaningful trends [7].

To obtain insights regarding shifts in global fire activity, studies have to incorporate multi-decadal time spans and continental-scale study areas. Global, long-term burnt area datasets are readily available for analyses, the most widely used ones being the National Aeronautics and Space Administration (NASA) MCD64A1 dataset [8] and the European Space Agency (ESA) Fire\_cci BA 5.1 dataset [9,10]. A third global burnt area product from the Global Fire Emissions Database, version 4 (GFED4, [11]) is meanwhile discontinued and is only available until 2015 [12]. A semi-automatically generated product is also provided by the Joint Research Center of the European Commission (JRC) in the frame of the European Forest Fire Information System (EFFIS) [13]. However, these data are only available for Europe, Northern Africa, and the Middle East [14].

While the listed datasets all include the fire perimeter as well as the burning detection date, they do not feature information regarding the fire severity, and thus do not allow the derivation of trends in this regard. Yet, fire severity is a critical aspect of fire regimes, determining fire impacts on ecosystem attributes and associated post-fire recovery [15].

The German Aerospace Center (DLR) recently published a global, long-term burnt area dataset, which includes information regarding the burn severity, together with fire perimeter and the burning date. This dataset is closely linked to the burnt area monitoring service operated by DLR. It is maintained by the Department of Geo-Risks and Civil Security (GZS) of the German Remote Sensing Data Center (DFD). The service is based on mid-resolution Sentinel-3 Ocean and Land Color Instrument (OLCI) satellite imagery and provides burnt area information for the region of Europe twice a day in near-realtime. The service is fully automated and targeted at supporting rapid mapping activities and timely post fire damage assessment throughout Europe. A quality-optimized version called fusion product is generated after a time delay of 10 days, when additional post-event data is available.

The recently published, global dataset is build upon this fusion product. In addition to the data available for Europe, equivalent products are generated using the same methodology for North and South America, Africa, Oceania, and Asia. The methodology is briefly described consecutively, the complete description can be found in Nolde et al. [16]. As the dataset is based on Sentinel-3 data, it is only available since 2016, which is the year the first Sentinel 3 satellite was launched. For this study, the data was extended using data from the NASA MODIS MOD09/MYD09 product [17] in order to allow the derivation of trends on a longer time scale.

The complete input dataset comprises 9612 granules of the MODIS MOD09/MYD09 product in conjunction with 3503 tiles of the OLCI instrument onboard the Sentinel-3 satellite. The area of Eastern Australia has been selected as a study region, covering the states and territories of Queensland, New South Wales, the Australian Capital Territory (ACT), and Victoria, respectively.

This study region is chosen because it experienced destructive burnings in the 2019/2020 fire season, and because it also was regularly affected by wildfires in recent decades.

The year 2019 was Australia's warmest year on record so far, with significant heat waves occurring in January and December [18]. The national, average maximum temperature was as high as 43.6 ◦C, which is more than 1.8 ◦C above the long-term average [19]. In addition, 2019 has also been the driest year on record in Australia, caused by an extraordinary strong positive Indian Ocean Dipole [20]. The nationally-averaged rainfall was 40% below average, amounting to only 278 mm. As a consequence, the annual, cumulative Forest Fire Danger Index reached the highest values since the beginning of the national records in 1950 [18]. Three of the four investigated states and territories, namely, New South Wales, the ACT, and Victoria, reside in a temperate climate zone, with dryness conditions usually reaching moderate levels at most [21]. However, in 2019, these states experienced severe drought conditions, with New South Wales suffering the most severe conditions

throughout Australia. These preconditions contributed significantly to the unprecedented fire activity in December 2019 and January 2020 [19].

This study analyses the existence of stable, wildfire related trends in this region, focusing on fire severity.

#### **2. Materials and Methods**

#### *2.1. Area of Interest*

The chosen study region comprises the states of Queensland, Victoria, and New South Wales together with the Australian Capital Territory (ACT). This region is vastly heterogeneous regarding climate and vegetation cover, with a pronounced inter-annual variability. Both Queensland and New South Wales include tropical, temperate, and arid climate zones. Figure 1 shows the area of interest, together with the respective climate zones and the burnt area from the 2019/2020 wildfire season. The highlighted climate zones are the ones found to feature increasing trends regarding fire severity in this study. The full names corresponding to the climate zone abbreviations can be found in the respective tables in the result section. This study is prepared in a hierarchical manner. Trends are analyzed on a state level, as well as regarding climate zones, and finally ecological units. These are consecutively set in relation to each other.

**Figure 1.** The area of interest, comprising Queensland, Victoria, and New South Wales together with the Australian Capital Territory (ACT). The climate zones are shown additionally, as well as the burnt area for the 2019/2020 fire season. The climate zones found to feature a significant increasing trend in this study are highlighted.

The utilized climate zone map, prepared by Beck et al. [22] following the methodology of Peel et al. [23], is derived using a long-term time series of weather station data regarding monthly precipitation and air temperature. The system classifies climatic regions

into five main classes and 30 subtypes. The discrimination between classes is based on fixed thresholds addressing the seasonality of precipitation and temperature. Climate is recognized as the major driver of global vegetation. The classification is therefore regarded as an empirical mapping of biome distributions around the world. Although developed in the 19th century, it is widely used today, for example, in ecological modelling [22].

The Australian Forest Fire Danger Index [24] relies on four meteorological parameters, which have proven their reliability with regard to wildfire activity: temperature, wind speed, relative humidity, and the Drought Factor, a component representing fuel availability [25].

The latter parameter, which is strongly influenced by seasonal variations in rainfall, has been found to be pivotal for fire occurrence, together with the vegetation structure [26]. As described by Russel-Smith et al., fire frequency as well as fire extent throughout Australia are strongly influenced by rainfall seasonality. Fire occurrence is therefore most pronounced in the tropics of Northern Australia, which are intensely seasonal. Rainfall positively influences the dynamics of biomass growth, which provides the fuel for wildfire activity. Precipitation amounts above average could be attributed to large burnings in arid, central Australia [27].

Dryness, on the other hand, strongly increases the availability of the vegetation to burn [28]. It could be shown that drought conditions are associated with major fires in the forested areas of Southern Australia [29,30]. These ecosystems usually feature sufficient litter for propagation of fire at most times. The most influential factor for fire propagation is the availability of vegetation to burn, controlled by drought conditions and the weather at the time of ignition [28].

#### *2.2. Utilized Data Sources*

The study utilizes the DLR-GZS burnt area dataset, which is based on mid-resolution optical satellite data from two different sensors. First, the OLCI instrument onboard the Sentinel-3 A and B satellites of the European Copernicus Programme [31], and second, the MODIS instrument onboard the NASA Aqua/Terra satellites. Band information from the red and near-infrared (NIR) domain are utilized for the retrieval of burnt area perimeters and burn severity estimation. Imagery of the Sentinel-3 OLCI instrument can be retrieved via the Copernicus Open Access Hub [32] and the Copernicus online data access website [33]. However, the available time span for OLCI data is considerably shorter than the one regarding MODIS. Sentinel-3A was launched in 2016, with Sentinel-3B following in 2018. Therefore, the OLCI dataset is extended by imagery from the MODIS sensor for this study. This data is available for an extended time range of more than two decades, starting in late 1999 with the launch of the NASA Terra satellite. The sensor is designed to observe the ocean, atmosphere, land, and ice on the Earth's surface [34]. It features 36 discrete spectral bands with differing spatial resolutions from 250 m to 1 km. Terra's twin satellite, Aqua, was launched in 2002, carrying a second MODIS instrument. This study makes use of MODIS information provided through the MOD09A1/MYD09A1 surface reflectance product [17], which represents a cloud-free 8-day composite of bands in the visible spectrum. The data are freely available from the NASA Land Processes Distributed Active Archive Center (LP DAAC) [35]. Furthermore, thermal anomaly information derived from the MODIS MOD14A2/MYD14A2 [36] product is used as auxiliary data. Equivalent data is utilized from the Visible Infrared Imaging Radiometer Suite (VIIRS) instrument onboard the Suomi National Polar-orbiting Partnership (SUOMI-NPP) satellite [37]. Both products are available for download on the NASA Fire Information for Resource Management System (FIRMS) website [38].

In order to reduce the total data volume to a significant selection, only the Australian summer months from November to February have been analyzed in this study.

Subsequent to the derivation of burnt area perimeters and the burn severity, the results are combined with land use and land cover (LULC) data to gain insights regarding affected vegetation classes. For this purpose, the CCI-LC (Climate Change Initiative-Land Cover) product from ESA is used. It provides mid-resolution land cover information on a global scale. Moreover, land cover maps are available for each individual year, starting in 1992 [39]. To investigate the relationship between fire activity and climatic conditions, a climate zone map based on the Köppen-Geiger classification system is utilized, published by Beck et al., 2018 [22]. Furthermore, ecological units provided by the United States Geological Survey (USGS) are used [40].

Two further burnt area datasets have been incorporated within this study to validate the presented results. First, the NASA MCD64A1 dataset, featuring global burnt area information derived from MODIS imagery, with a spatial resolution of 500 m [8]. These data are made available by the University of Maryland [41] and has been widely utilized in academic research, for example, for the Brazilian Savannas [42,43]. Second, a high resolution burnt area map for the state of New South Wales is utilized, which was prepared by the Department of Planning, Industry and Environment of New South Wales/Australia. This dataset, named Google Earth Engine burnt area map (GEEBAM) [44], is based on Sentinel-2 data and makes use of manually derived thresholds from aerial photography [45]. The National Indicative Aggregated Fire Extent Datasets [46], which are published by the Australian Government, have also been taken into consideration as a reference data source. However, the GEEBAM dataset has been found to feature a higher thematic accuracy, and is therefore chosen as reference in this study.

Table 1 lists the complete set of available MOD09/MYD09 composites and Sentinel-3- OLCI scenes used for generating the DLR-GZS burnt area dataset, which this study is build upon. Besides the number of available scenes for each time range and state, the average number of cloud free observations per pixels is given, representing a measure of the interpretability of the data. As the MOD09/MYD09 data are available as an eight-day composite, the number of cloud-free observations is considerably lower when compared to Sentinel-3. In total, 9612 MODIS MOD09/MYD09 granules from both Terra and Aqua have been analyzed, together with 3503 OLCI scenes from Sentinel-3 A and B. This amounts to an entirety of 13,115 scenes for the complete study time span.


**Table 1.** Analyzed data sources, listed separately for Queensland (QLD), New South Wales (NSW), Australian Capital Territory (ACT), and Victoria (VIC).

#### *2.3. Burnt Area Derivation Methodology*

The accurate, automatic monitoring of burnt area evolution from satellite imagery represents a demanding task for the scientific community. This is mostly due to the spatialtemporal variability of the state of the Earth's surface. Inaccuracies are also introduced by the utilized sensor and related geometrical resolution [47]. Furthermore, the presence of clouds disturbs the derivation of meaningful surface features at reflective wavelengths [48]. Adding to the complexity inherent in optical sensor data, burnt areas are highly heterogeneous regarding size, shape, and spectral reflectance on the ground surface, and can thus be difficult to differentiate from shadows cast by clouds and mountains. These circumstances significantly limit the possibilities for an automated approach, especially on a large geographical scale. Many recent research activities have assessed these challenges using Machine Learning technologies. Methodologies such as Random Forests [49], Support Vector Machines [50], and Deep Learning classifiers [51,52] have been applied in this regard. However, research has mainly been focused on feature detection using high resolution satellite imagery on a small geographic scale. Due to the spatial resolution however, these satellites inherently feature a low temporal frequency, and are thus unsuited for daily monitoring purposes. Approaches utilizing mid-resolution satellite data, thus allowing a high temporal coverage for an extended geographical region, usually use a time series approach in combination with a burn-sensitive vegetation index. Thermal anomaly detections are often used as auxiliary data (see [53]). However, these approaches feature a significant amount of uncertainty. Humber et al. analyzed four global burnt area products, and concluded that the estimates of burned area vary greatly between products in terms of total area affected, the location of burning, and the timing of the burning [54]. In a similar study, Padilla et al. found that the commission error ratio was above 40% and omission error ratio was above 65% for the analyzed products [55]. Oliva et al. conducted a study investigating if thermal anomaly data could be used as a replacement for burnt area datasets, but found high omission and commission errors especially for grasslands, savannas, and agricultural areas [56].

The methodology utilized here is based on an approach developed at DLR [16]. It is primarily designed for monitoring continental-scale regions in near-real time, but can also be invoked to perform retrospective time series analysis. The derived information comprises burnt area perimeters, the date of detection, and the burn severity by means of the differential Normalized Difference Vegetation Index (*NDV Idi f f* [57], see Equations (1) and (2). As auxiliary information, the number of detections for each burnt area pixel is available, as well as the number of cloud-free satellite overpasses for each pixel.

$$NDVI = \frac{NIR - Red}{NIR + Red} \tag{1}$$

$$NDVI\_{diff} = NDVI\_{prc} - NDVI\_{post} \tag{2}$$

At its core, the method exploits the synergetic effects of data from the red/NIR and the thermal wavelengths in order to derive burnt area information. Substantial work in this regard has been performed by Fraser et al. and Li et al. for the boreal forest of Canada as early as the year of 2000 [58,59].

As a basis for the processing, mosaics of pre- and post-*NDV I* information are generated. Consecutively, the concept of Morphological Active Contours without Edges (MorphACWE [60,61]) is used to derive accurate burnt area perimeters. The method is closely related to Geodesic Active Contour Level Sets [62], which have been used for burnt area derivation [63] as well as other domains, such as crop field size estimation [64].

The term Active Contour refers to a dynamical curve, which grows starting from a set of seed pixels and converges when an optimal segmentation result is reached. For the generation of the burnt area dataset, Active Fire locations are used as seed information. This proceeding was shown to yield results of high geometric accuracy when inter-compared with the JRC/EFFIS dataset [13] as well as the NASA MCD64A1 dataset [8]. The accuracy validation, together with the detailed description of the methodology, can be found in Nolde et al. (2020) [16]. The methodology is schematically visualized in Figure 2.

**Figure 2.** Scheme of the methodology for burnt area perimeter derivation, from Nolde et al. (2020) [16].

Figure 3 shows exemplary visualizations of the data used for the preparation of this study. Sub-figure (a) visualizes the maximally detected burn severity regarding the mega fire in the greater area of Sydney/New South Wales, which was active during November and December, 2019. Sub-figure (b) shows the temporal evolution of the burn activity, by means of detection date.

**Figure 3.** (**a**) *NDV Idi f f* values regarding the mega fire in the greater area of Sydney/New South Wales, 2019/2020. (**b**) Temporal evolution of the fire activity.

#### *2.4. Validation of Burnt Area Data*

The presented results are compared against the two reference datasets, NASA MCD64A1 and GEEBAM, regarding three criteria:


The fourth criterion, which represents the true negatives (TN), refers to the percentage of area neither contained in the presented results, nor in the reference data. However, as the total size of unburnt area greatly predominates the total size of burnt area, this true negative percentage is implicitly very close to 100 percent. This is, however, mostly due to the size of the study region, so this measure does not represent a meaningful value for this kind of study. The same applies to the overall accuracy [65]. These measures have therefore been omitted in Table 2, which shows the results of the evaluation. A more suited means of measure is to calculate the average of the true positive percentage and the inverted false positive percentage (false positives subtracted from 100%, named FPinv in Table 2).

The inter-comparison of the burnt area of the presented results with the burnt area of the high resolution GEEBAM 2019/2020 data for New South Wales reveals a percentage of overlapping area of 77%, with an error of 9%. The Jaccard index, also known as Intersection over Union, yields a similarity of 70.8% [66]. To enable an evaluation of these numbers, the NASA MCD64A1 dataset is equivalently checked against the same reference, yielding 71% overlap and an error of 7%. The combination of true positives and inverted false negatives account for 84% for the presented results, and 82% for the MCD64A1 data. The results obtained in this study regarding the 2019/2020 fire season are in accordance with burnt area extent information published by Boer [67]. The data basis is therefore considered to be of satisfactory accuracy.

**Table 2.** Inter-comparison of burnt area extent with NASA MCD64A1 and GEEBAM reference datasets.


#### *2.5. Trend Derivation Methodology*

The trend derivation in this study is based on linear regression, whereby the slope of the regression line represents the actual trend. The input values for the burnt area extent analysis are the accumulated burnt area amounts per unit of investigation (state, climate zone, or ecological unit) for each year of the analyzed time span. Regarding the burn severity analysis, the input values are derived by averaging all values within the unit of investigation for each respective year. The correlation coefficient indicates to what extent the actual values are in concordance with the calculated trend line, while the RMSE (Root Mean Square Error) illustrates the error.

Finally, the 5- and 95-percentiles (named "Perc 5" and "Perc 95" in the result tables) represent the error margins, indicating how robust the result actually is. They are derived through repeatedly and randomly altering the input values within the range of the RMSE, so that the results reflect the average of a set of possible outcomes. To eliminate outliers, the 5 and 95-percentiles are used as upper and lower limits of the yielded results.

#### **3. Results**

The consecutive sub-sections show the results of the analysis on a state-wide level, as well regarding climatic zones and finally ecological units. These results are then set in relation to each other. Each subsection contains a visualization of the burnt area extent for each year within the analyzed time span, followed by the actual burn severity analysis.

#### *3.1. Fire Trends Regarding the States in the Study Area*

Figure 4 shows the total, annual extent of burnt area for the four Australian states, and territories of interest, regarding the time period of 2000 to 2020. The results derived from Sentinel-3 OLCI data are depicted as a green line, while results regarding MODIS MOD09A1/MYD09A1 data are shown in blue. As the first one of the Sentinel-3 satellites was only launched in 2016, the analysis could only be carried out for this limited timespan. The red, dotted line represents the NASA MCD64A1 [8] burnt area dataset. The latter is included as a reference, to allow an estimation of the accuracy of the presented results. Finally, the black, dotted line is the regression line, corresponding to the MODIS MOD09/MYD09 based results.

**Figure 4.** Total yearly burnt area amount in million hectares for New South Wales, Victoria, Queensland, and the Australian Capital Territory (ACT). MCD64A1 reference data is additionally visualized as dotted, red line. The black, dotted line in the subplot for Queensland represents the linear regression line.

It can be seen that the results show a high correlation of burnt area extent between the utilized DLR-GZS burnt area dataset and the MCD64A1 reference, regarding three of the four analyzed states. For Queensland, however, this is not the case. The discrepancies are due to differences in the methodologies: Unlike in the NASA dataset, burn sites where the vegetation has quickly recovered are excluded in the presented data. As stated before, wildfire is a natural phenomenon, and in many parts of Australia the regular burning of huge areas of bush and grassland vegetation is part of the natural cycle. The affected vegetation recovers quickly, some species do so even within a period of a few weeks. In the case of Queensland, these areas account for the majority of the overall burnt area. In order to confine the results to potentially harmful wildfire events, it was decided to consider burnt area in this study only when distinctive traces of burning activity could still be detected after a period of three weeks. This filtering is not performed for the generation of the MCD64A1 dataset, and as a consequence, the overall burnt area for Queensland differs considerably between the two products.

Furthermore, the figure shows that there is no linear trend in fire extent regarding New South Wales, Victoria, and the ACT. This is also reflected by Table 3, which lists the respective statistics for Figure 4. The *p*-value, symbolizing the statistical significance, illustrates that there is no linear development regarding the derived trends. Complying with the common standard, a *p*-value below 0.05 is regarded to represent statistical significance. Such a low value is found in the case of Queensland, which features a robust, increasing trend regarding the extent of wildfires over the last two decades.

**Table 3.** Trend statistics regarding the extent of burnt area from 2000 to 2020. Only Queensland shows a statistically significant upward trend (highlighted in dark gray).


A reason for the extraordinary extent of burnt area, which is especially pronounced for the state of New South Wales, is shown in Figure A1 in the Appendix A. The burnt area extent is depicted for the past 20 years, sub-divided by month of the fire season. It can be seen that the fire activity in New South Wales reached a significant level at the beginning of November, while the main activity usually only occurs towards the end of the year. The fire activity started several weeks earlier than usual in the 2019/2020 season.

While the figures above provide a general impression by showing the fire extents over the past 20 years, this study focuses on fire severity, which is analyzed consecutively. The term is defined as a measure of the degree of environmental change caused by fire [68]. It represents a critical aspect of fire regimes, indicating the impacts on ecosystems and associated post-fire recovery [15]. The respective impact ranges from the partial consumption of litter to the complete dieback of canopy trees [69].

Figure 5 shows the yearly trend of average burn severity for the four analyzed states and territories. The green line depicts the severity results regarding Sentinel-3 OLCI, while the results for the MODIS data are depicted in blue. For the latter data source, not only the burnt areas have been investigated, but also the complete complementary area, which has not been affected by fire. This represents a cross check, showing that the developments in severity are not actually caused by unrelated factors such as climate, soil or moisture related conditions. It can be seen that the states of New South Wales and Victoria feature a general upward trend, while the development is stable or even decreasing for the ACT and Queensland.

Table 4 shows that, even if the severity trend for New South Wales and Victoria is increasing, the *p*-values are too heterogeneous for these trends to be statistically significant. Only the slightly negative trend for Queensland can be considered robust. The results of the cross check are listed in Table 5. The mentioned heterogeneity is addressed in the Discussion section.

**Table 4.** Trend statistics regarding the mean severity of burnt area from 2000 to 2020. Queensland features a statistically significant trend (highlighted in dark gray), while Victoria shows a trend close to statistical significance (depicted in light gray).


**Table 5.** Cross check trend statistics regarding the mean decline in vegetation fitness of area unaffected by fire, from 2000 to 2020.


**Figure 5.** Yearly average burn severity for the period of 2000 to 2020. The line depicted in light blue represents the mean severity values for the burnt areas derived from MODIS MOD09/MYD09 data. The dashed, dark blue line shows the equivalent results for areas not affected by fire. Green color is used for Sentinel-3 OLCI burnt area results (available only since 2016). The black, dotted line is the regression line, with respect to the MODIS burnt area results.

In order to discriminate regions which account for the rising trend in New South Wales and Victoria, the study region is subdivided into regions which share similarities regarding available fuel and species composition. These are consecutively analyzed equivalent to the states. As a reasonable classification, climate zone mapping information is used. Parks et al. could show that climatic conditions, next to fuel and weather, represent a major driver of fire intensity [70]. As a next step, the climate zones featuring statistically relevant trends are further subdivided into ecological units, which are provided by the United States Geological Survey (USGS) [40]. This allows the fine-granular attribution of state-wide trends to small scale vegetation types.

The following section analyzes the climate zones contained in the study area regarding their burn severity trends. For the zones featuring developments with statistical significance, it is shown to which extent they overlap with the area of New South Wales and Victoria.

#### *3.2. Fire Trends Regarding Climate Zones*

Figure 6 shows, first and foremost, the exceptionality of the 2019/2020 wildfire events, regarding the total extent of the burnings. The calculated size reaches 6.5 million hectares, which is more than twice as high as in every other year in the analyzed time span. In addition to the total extent, a subdivision regarding the affected climate zones is shown. For this division, the Köppen–Geiger classification system published by Beck et al. (2018) [22] is used. The color scheme follows the one proposed by the authors.

**Figure 6.** Total yearly burnt area amount in million hectares, subdivided by climate zones.

Table 6 represents a listing of all climate zones in the area of interest together with their *p*-value as a measure of statistical significance. A threshold of *p* ≤ 0.05 is applied to indicate the significance, the respective table row is marked in dark gray. Climate zones featuring a *p*-value close to statistical significance (0.05 < *p* ≤ 0.1) are highlighted with light gray color.

**Table 6.** Trends regarding the size of affected area for each climate zone in the area of interest. The climate zone featuring a statistically significant trend is marked in dark gray. The one showing a trend close to statistical significance is depicted in light gray.


Note that, apart from the temperate zones featuring dry winters and warm to hot summers (Cwa/Cwb), all climate zones existing in the Eastern part of Australia are affected by wildfire.

The table shows that only one climate zone satisfies the *p*-value condition for statistical significance (Cfa: Temperate, no dry season, hot summer). The respective correlation coefficient lies in the moderate range, even if the RMSE (Root Mean Square Error) shows a considerable oscillation around the regression line. The column "Slope (%)" shows the actual inclination of this line, given in percent. This value represents the actual trend. For this climate zone, the yearly rate is 3.4%, indicating a considerable increment in fire size over recent years.

As can be seen in Figure 7, robust trends regarding fire severity can be derived for two climate zones. First, the temperate zone featuring dry and hot summers (Csa), and second, the arid desert zone featuring cold conditions throughout the year (BWk). The first one shows an inclination of 0.42% per year on average, the second one features a value of 0.11%. The zone of arid steppe with year-round cold conditions (BSk) shows the second largest positive trend inclination, but features a *p*-value just above the threshold for statistical significance.

**Figure 7.** Trends regarding fire severity for each climate zone. The lengths of the bars visualize the strength of the trend. The bars for climate zones featuring a statistically significant trend are drawn with a thick, black border. A gray border is used to identify trends close to statistical significance. The horizontal, black lines represent the error margins.

Table 7 lists the severity trend for each climate zone.


**Table 7.** Trends regarding fire severity for each climate zone, given in percent of increase/decrease. Climate zones featuring a statistically significant trend are marked in dark gray. The one showing a trend close to statistical significance is depicted in light gray.

In order to verify that the trends depicted in the above figures and tables are actually connected to fire occurrence, instead of being a general phenomenon or an effect by a cause not investigated, a cross-check has to be performed. Figure A2, which is located in the Appendix A, shows severity trends for each climate zone, where only areas are considered that have not been affected by fire. As can be expected, no general trend of increased fire severity is observable. In fact there is a generally negative development, indicating a general increase in vegetation fitness. The *p*-values and the correlation coefficients are generally low, meaning that there is no connection between vegetation fitness and the progression of years in these unaffected areas.

Table A1, which can be found in the Appendix A, shows the statistical information regarding the cross-check.

#### *3.3. Fire Trends Regarding Ecological Units*

In order to draw conclusions regarding the vegetation types causing the increasing severity trends in some of the climate zones, the analysis is also carried out on the basis of ecological units. These units feature a higher spatial and thematic resolution, and are thus better suited for analyses on a smaller scale.

For the incorporation of ecological units, the Global Ecological Land Units global dataset provided by the United States Geological Survey (USGS) [40] is used. The ecological units are a combination of bioclimate region, landform type, surficial lithology, and land cover information [71], and allow for a very high thematic resolution. This results in more than 3600 different units covering the area of interest. In order to reduce the number of units to be analyzed to an appropriate level, a subset is generated from the original data in a preceding step. This subset contains all ecological units which were affected by the 2019/2020 burnings and covered more than 1% of the burnt area. Furthermore, it comprises all ecological units covering more than 1% of the area of interest.

Figure 8 shows the analysis results of the total burnt area extent development, regarding the ecological units. It can be seen that the three most prominent units (1750, 1529, and 1664) feature either needleleaf or evergreen forest, which has not been the case in former years. Statistically significant results can, however, only be derived for the class of "Hot Wet Mountains on Non-Carbonate Sedimentary Rock with Mostly Needleleaf/Evergreen Forest" (2268). The statistical results are found in Table A2 in the Appendix A.

**Figure 8.** Total yearly burnt area amount in million hectares, subdivided by ecological units.

Figure 9 shows the severity trends regarding ecological units, equivalent to the climate zone analysis above. Similar to the burnt area extent analysis, stable trends can be derived for two classes featuring either needleleaf or evergreen forest. These classes are "Hot Wet Mountains on Non-Carbonate Sedimentary Rock with Mostly Needleleaf/Evergreen Forest (2268)" and "Warm Wet Mountains on Metamorphic Rock with Mostly Needleleaf/Evergreen Forest (1652)". Two further units feature *p*-values close to statistical significance, and are thus worth being considered: "Warm Semi-Dry Plains on Unconsolidated Sediment with Mostly Cropland (1712)" and "Hot Moist Plains on Unconsolidated Sediment with Grassland, Scrub, or Shrub (2529)". All respective statistical information is found in Table A3 in the Appendix A.

**Figure 9.** Trends regarding fire severity for the most affected ecological units. The lengths of the bars visualize the strength of the trend. The bars for ecological units featuring a statistically significant trend are drawn with a thick, black border. A gray border is used to identify trends close to statistical significance. The horizontal, black lines represent the error margins.

The increasing burn severity for some classes can be linked to a higher degree of combustion. Other reasons include the higher amount of combustible biomass, the concern of younger, healthier vegetation, or the exposure of different vegetation types.

The cross-check, conducted for areas which have not been affected by fire, shows the expected, generally negative trend. The derived fire severity trends presented above are therefore demonstrably caused by wildfire activity. Figure A3 shows the trend for each ecological unit. All respective statistical information is found in Table A4 in the Appendix A.

Figure 10 shows a section of the two ecological units with a statistically significant increase in fire severity during the analyzed time span. The figure represents a detailed

view of the North-Eastern part of New South Wales, where these two ecological units overlap with the area affected by the 2019/2020 wildfires.

**Figure 10.** Ecological units featuring a statistically significant severity trend, overlapping 2019/2020 wildfire locations in northern New South Wales.
