2019–2020 Australia Fire and Its Relationship to Hydroclimatological and Vegetation Variabilities

: Wildfire is a major concern worldwide and particularly in Australia. The 2019–2020 wildﬁres in Australia became historically signiﬁcant as they were widespread and extremely severe. Linking climate and vegetation settings to wildﬁres can provide insightful information for wildﬁre prediction, and help better understand wildﬁres behavior in the future. The goal of this research was to examine the relationship between the recent wildﬁres, various hydroclimatological variables, and satellite-retrieved vegetation indices. The analyses performed here show the uniqueness of the 2019–2020 wildﬁres. The near-surface air temperature from December 2019 to February 2020 was about 1 ◦ C higher than the 20-year mean, which increased the evaporative demand. The lack of precipitation before the wildﬁres, due to an enhanced high-pressure system over southeast Australia, prevented the soil from having enough moisture to supply the demand, and set the stage for a large amount of dry fuel that highly favored the spread of the ﬁres.


Introduction
Wildfires are one of the major environmental concerns in Australia that are widespread and occur regularly. They have significantly shaped the nature of the continent over millions of years. However, the wildfires that occur virtually every summer have devastating outcomes. The 2019-2020 wildfire season killed 33 people, destroyed over 3000 homes, razed almost 19 million hectares, and exterminated about one billion animals, including several endangered species. The 2019-2020 Australian bushfire season began with several serious uncontrolled fires in June 2019. Throughout the summer, hundreds of fires burned in southeast Australia with the major peak during December and January. In eastern and northeastern Victoria, large forest areas burned at an unstoppable rate for four weeks. There has been considerable debate regarding the underlying cause of the intensity and scale of the fires, including climate change, which brought significant international attention to this event [1][2][3].
Wildfires have been extensively studied both in Australia and around the globe. There is significant evidence that demonstrates a growth in the frequency, spatial extent, and severity of wildfires in several areas around the globe. For instance, Jolly et al. [4] used three daily global climate datasets and three fire danger indices to develop a simple annual metric of fire weather season length, and map spatiotemporal trends from 1979 to 2013. They showed that fire weather seasons have lengthened across 29.6 million km 2 (25.3%) of the Earth's vegetated surface, resulting in an 18.7% increase in the global mean fire weather season length. Such tendencies have also been linked to climate change. For example, Bradstock et al. [5] modeled fire events using FIRESCAPE to simulate landscape fires, weather and fuel dynamics, and treatments in the Sydney region, indicating that climate change affects fire regimes extensively. In addition, Liu et al. [6] investigated the trend in global wildfire potential under climate change due to the greenhouse effect using General Circulation Models (GCMs). They reported that future wildfire potential increases significantly in the United States, South America, Central Asia, southern Europe, southern Africa, and Australia.
One way to analyze the relationship between wildfires and hydroclimate is to observe signals in climate variables before and during the wildfires. The relationship between climate, meteorological conditions, and fires is well known; the most destructive fires are usually preceded by extremely high temperatures, low relative humidity, and strong winds, which combine to create ideal conditions for the rapid spread of fire [7]. In recent decades, the climate has changed by an increase in the mean temperature on a global scale [8]. Because this increase has been especially noticeable since the 1970s, it is possible to examine whether the changing climate has already affected the frequency and magnitude of wildfires [7,9]. In particular, Australia shows one of the most variable climates due to inter-annual rainfall caused by El Niño Southern Oscillation (ENSO) [10].
The main issue in using climate variables for the study of wildfires is the lack of data availability in many places, especially when long records are required. Lack of spatially-comprehensive climate data on regional variation in fire regimes, and statistical limitations associated with the temporal extent of instrumental records are among the many challenges of studying wildfires. Fortunately, remote sensing and observational products can help manage these challenges by providing both temporally and spatially comprehensive datasets. Remotely sensed data are useful for the identification of the extent and intensity of wildfires. For example, Cattau et al. [11] used Moderate Resolution Imaging Spectroradiometer (MODIS) Fire Radiative Power (FRP) data from 2003 to 2016 to show that fires, especially in the western United States, are becoming larger and more frequent over time. In a similar study, Vadrevu et al. [12] mapped the trends of fires in southeast Asia, concluding that the highest fire count occurs in Indonesia. Hence, the MODIS FRP dataset is a reliable source for identifying the extent of the fires. In addition, Vadrevu et al. [13] quantified the fire hotspots in India using the fire pixel from the Fire Information for Resource Management System (FIRMS). In addition to supporting the use of the FIRMS dataset, they found that the intensity of fires in evergreen forest regions is higher compared to other surface types.
A large-scale synoptic pattern that can be represented by 500 hPa geopotential height along with 10-m wind is often studied in relation to wildfires [14,15]. The 10-m wind brings dry air mass to wildfire-prone areas, making the environment more favorable for wildfires. Garcia-Ortega et al. [14] demonstrated that the 500 hPa geopotential height is helpful to characterize the synoptic pattern during wildfire season in Spain. They showed that wildfire events caused by lightning are closely related to the triggering of convective clouds, which is directly linked to a favorable synoptic environment for the convective organization. Zhong et al. [15], using statistical analysis (i.e., composite Empirical Orthogonal Function (EOF) and self-organizing map), found that the dominant anomalous 500 hPa synoptic pattern associated with the occurrence of large wildfire events over the northwestern United States is characterized by positive height anomalies over the western United States. They also noticed that the environment becomes more favorable for wildfires due to the dry air mass transported by the dominant anomalous 10-m winds from central Canada.
The relationship between wildfires and local hydroclimatic variables is investigated in several studies as well. Well known is the impact of the root zone Soil Moisture (SM) in the availability of water required by plants to perform photosynthesis. Long-duration droughts introduce stress on vegetations, reduce vegetation water content, and amplify the risk of wildfires. Several studies have linked SM to the risk of wildfires, including the study by Krueger et al. [16] that compared an SM-based index against the traditional ones showing a better predictive capacity. Similarly, Chaparro et al. [17] studied SM with temperature trends to predict the extent of wildfires. The impact of the precedent SM state also has a large impact on how prone to wildfires the different zones are, and SM in combination with other variables has shown potential to be used for fire prediction [18][19][20][21]. Temperature also influences the wildfires' behavior dramatically. Specifically, drought and high temperatures are strongly correlated to potential megafires in areas with dry vegetation [22].
In the same context, global climate change models predict warmer and drier conditions in the coming decades, which will increase the frequency, size, and intensity of fire events [23]. According to Balling et al. [24], the increasing temperature in the fire season, decreasing precipitation, and increasing drought conditions in the antecedent season increase the likelihood of wildfires. Recent studies have supported this relationship as well [25][26][27][28]. Precipitation is one of the most important factors controlling fuel moisture content and plant growth, and there is a high correlation between phenological cycle peak, integrated greenness, and monthly Southern Oscillation Index (SOI). Evapotranspiration is the only connecting component of water and energy balances. Due to the complex interactions between soil, water, vegetation, and atmosphere, Actual Evapotranspiration (AET) is probably one of the most impacted variables by wildfire [29][30][31]. There are limited studies that have focused on the wildfire-evapotranspiration relationship. Poon et al. [32] explored the spatial and temporal changes in evapotranspiration following the 2011 Las Conchas fire in New Mexico (United States) using the Operational Simplified Surface Energy Balance (SSEBop) model. They reported an average annual decrease of 120 mm of evapotranspiration within the regions affected by the Las Conchas fire. In addition, Vapor Pressure Deficit (VPD) can provide valuable information for drought early detection and can be used as the drought proxy for wildfire analysis [33][34][35].
Vegetation indices derived from optical remote sensing have been used to monitor aboveground biomass, vegetation density, and live fuel moisture for many regions throughout Australia [36][37][38][39]. Normalized Difference Vegetation Index (NDVI) and Leaf Area Index (LAI) can be used as proxies for vegetation properties and utilized by algorithms to retrieve values of fuel mass and water content [40,41]. NDVI and LAI have also been used directly by models to predict wildfire risk, severity, and extent [42,43]. During the growing season, high NDVI and LAI values generally indicate healthy and productive vegetation, while decreasing values generally indicate that vegetation is unhealthy and stressed due to moisture limitations or environmental disturbances that inhibit photosynthesis or reduce the growth of leaves [44,45]. In this case, areas that exhibit anomalously low NDVI and LAI can be attributed to vegetation that is drier and more flammable, thus signaling an increased risk of a wildfire [44,45]. Predominant vegetation classes must be considered when deriving live fuel moisture at large scales due to the spatial variability of plant physiology and phenology [46]. However, optical vegetation indices are limited by cloud cover and saturate quickly in dense vegetation, making it challenging to monitor fire risk for many regions throughout Australia from optical remote sensing alone [47]. Vegetation Optical Depth (VOD) from passive microwave observations experiences minimal attenuation by the atmosphere and is much slower to saturate in dense canopies [48]. While optical indices return information from the top surface of vegetation only, VOD is derived from the energy that passes through the entire vertical structure of vegetation and is more directly sensitive to the amount of water contained in both its photosynthetic and non-photosynthetic components [47,49]. Furthermore, VOD products derived for separate microwave channels can be used to distinguish and monitor specific structural features of vegetation (e.g., woody trunks and stems vs. leaves) [50,51]. Thus, VOD can be used in tandem with indices such as NDVI and LAI to enhance the measurement of biomass, live fuel moisture content, and to assess fuel conditions before wildfire events [52,53]. The main objective of this study was to understand the relationships between the 2019-2020 wildfire season and its climate and vegetation settings. This understanding can help to improve the prediction of the severity of the wildfires in the future and provide some insights into the behavior of future wildfires under climate change.

Study Area
The Region of Interest (ROI) for this study was defined as the burned areas over Australia from December 2019 to February 2020. The burned areas were obtained from MODIS-MCD64A1 collection 6 product, identified through tracking rapid changes in daily surface reflectance time-series. These changes in surface reflectance are due to charcoal and ash deposits, removal of vegetation, and alteration of vegetation structure during fire events. As the land-cover types in the burned area could be grouped into two very different main classes, including short vegetation (e.g., grasslands and open shrublands) and evergreen forest, the portion of the total burned area covered with evergreen forest (17.84%; Supplementary Materials Figure S1) was analyzed beside the total burned area. Figure 1 shows that the majority of the burned areas located in southeast Australia were related to evergreen forests (including both broadleaf and needle leaf) that were close to populated residential regions, supporting the decision to investigate the burned evergreen forests separately. The main objective of this study was to understand the relationships between the 2019-2020 wildfire season and its climate and vegetation settings. This understanding can help to improve the prediction of the severity of the wildfires in the future and provide some insights into the behavior of future wildfires under climate change.

Study Area
The Region of Interest (ROI) for this study was defined as the burned areas over Australia from December 2019 to February 2020. The burned areas were obtained from MODIS-MCD64A1 collection 6 product, identified through tracking rapid changes in daily surface reflectance time-series. These changes in surface reflectance are due to charcoal and ash deposits, removal of vegetation, and alteration of vegetation structure during fire events. As the land-cover types in the burned area could be grouped into two very different main classes, including short vegetation (e.g., grasslands and open shrublands) and evergreen forest, the portion of the total burned area covered with evergreen forest (17.84%; Supplementary Materials Figure S1) was analyzed beside the total burned area. Figure 1 shows that the majority of the burned areas located in southeast Australia were related to evergreen forests (including both broadleaf and needle leaf) that were close to populated residential regions, supporting the decision to investigate the burned evergreen forests separately.  Table 1 summarizes datasets used in this study. The FRP was retrieved from the MODIS MOD14 version 6 product, which is a daily level 2 product generated with 5-min temporal swaths and pixel size of 1 km. Fire pixel counts (FPC) data retrieved from brightness temperature was obtained from  Table 1 summarizes datasets used in this study. The FRP was retrieved from the MODIS MOD14 version 6 product, which is a daily level 2 product generated with 5-min temporal swaths and pixel size of 1 km. Fire pixel counts (FPC) data retrieved from brightness temperature was obtained Water 2020, 12, 3067 5 of 18 from the FIRMS dataset. The 500 hPa geopotential height and the 10-m winds data were obtained from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis monthly means product. The spatial resolution of the dataset is 2.5 • × 2.5 • with 17 pressure levels, available from 1948 to the present. Access to reliable long-term records of soil moisture with sufficient resolution and spatial coverage is a challenge. Point measurements, although probably more precise, are sparse and have a very limited spatial representativeness. Satellite-based remotely sensed products of soil moisture offer an adequate balance between accuracy and spatial coverage, although the temporal extent of the records is not adequate to have a climatological perspective; for instance, retrievals from the Soil Moisture Active Passive (SMAP) are only available since 2015; while Gravity Recovery and Climate Experiment (GRACE) retrievals are not limited to the root zone, their resolution is too coarse for this purpose. On the other hand, models provide good resolution and larger records, although their accuracy varies largely among different models. The Global Land Data Assimilation System (GLDAS) version 2.1 was chosen for soil moisture due to its wide use, availability, and broad validation [54]. The root zone soil moisture is available globally at a resolution of 0.25 • × 0.25 • every 3 h from NASA. GLDAS is the result of the NOAH Land Surface Model, forced with the National Oceanic and Atmospheric Administration/Global Data Assimilation System atmospheric analysis fields, the disaggregated Global Precipitation Climatology Project (GPCP) precipitation field, and the Air Force Weather Agency's AGRicultural METeorological modeling system radiation [55,56].

Datasets
Air temperature data was recollected by using European Centre for Medium-Range Weather Forecasts (ECMWF) global reanalysis (ERA5). ERA5 covers the period from 1950 to the present with Water 2020, 12, 3067 6 of 18 hourly resolution. The native horizontal grid is 30 km and resolves the atmosphere using 137 levels from the surface up to a height of 80 km. Integrated Multi-satellitE Retrievals for GPM (IMERG) V06 product was used for precipitation analysis. IMERG estimates precipitation rate by merging, interpolating, and inter-calibrating data from various microwave and infrared sensors. Using multiple data sources allows this product to offer fair spatial and temporal resolution (i.e., 0.1 The MODIS ET product was used for exploring actual and potential evapotranspiration over the ROI [57]. TerraClimate dataset was used to calculate and investigate VPD over the ROI. TerraClimate is a high-spatial-resolution dataset (1/24 • ,~4 km) for global terrestrial surfaces from 1958-present. TerraClimate uses climatically aided interpolation, combining high spatial resolution climatological normals from the WorldClim dataset, with coarser resolution time-varying data from other datasets in order to produce a monthly dataset of precipitation, maximum and minimum temperature, wind speed, vapor pressure deficit, and solar radiation.
To assess the pre-fire condition of fuels (e.g., aboveground vegetation), several vegetation products derived using optical and passive microwave remote sensing were downloaded for the years 2003-2020. NDVI data from individual MODIS Terra and Aqua products (MOD131Q1.006 and MYD131Q1.006 respectively) and LAI data from a blended MODIS Aqua and Terra product (MCD15A3H.006) were used in this study. C-Band and X-Band VOD data from the Vegetation Optical Depth Climate Archive (VODCA) (Moesinger et al., 2020) were downloaded for the years 2003-2018. Land cover classes from MODIS (MCD12Q1) were downloaded and masked to the ROI.

Method of Data Processing and Analysis
In this study, anomalies of several variables were investigated for the period 2002-2020 for September-October-November (i.e., fire season in previous years and pre-fire season in the 2019-2020 fire season; SON) and December-January-February (DJF) to assess whether climate data reveal a fire season's response to climate variabilities. The selection of these variables and indices was based on a comprehensive review of the available literature. Selected variables were categorized into four groups. The first group refers to the quantification of the wildfires through the FRP and the burned area. The second group includes geopotential height, wind speed, and ENSO. The third group corresponds to hydroclimatic variables, including root zone soil moisture, near-surface air temperature, precipitation, potential evapotranspiration (PET), actual evapotranspiration, and water vapor pressure deficit. Finally, the fourth group consists of vegetation indices, including LAI, NDVI, Enhanced Vegetation Index (EVI), and VOD.
For spatial analysis of 2019-2020 fire events, EOF also known as Principal Component Analysis (PCA) was applied to the FRP in DJF from 2002 to 2020 to obtain the dominant mode of the fire intensity and distribution during the 2019-2020 fire season. The computation involved the multiplication of the annual FRP anomaly and the vector variances of the first principal component (PC1). Significant differences in explained variance of the first two principal components indicate the existence of significant dominant modes or patterns. To make the computation manageable, the FRP datasets were regridded from 1 km to 10 km. Only FRP pixels with detection confidence greater than 40% were included in the computation following the methodology of Vadrevu et al. [12]. PCA was also performed on the geopotential height and 10-m wind to obtain the dominant mode of both variables during the 2019-2020 fire events. FPC data was used for time-series anomaly analysis in SON and DJF for each year relative to the mean FPC in the study period. This time-series allows to track the trend and change of FPC in the total burned area and the forested subset over the years.
Oceanic Nino Index (ONI) was used to identify the ENSO phase. ONI is a 3-month running mean of the Extended Reconstructed Sea Surface Temperature (ERSST) anomalies in the Niño 3.4 region (5 • N-5 • S, 120 • E-170 • W) relative to a single fixed 30-year base period [58]. Simple time-series of the ONI values were compared with that of FRP pixel counts from 2002 to 2020.
Each of the remaining variables were first masked to the extent of the fires during the 2019-2020 fire season (i.e., ROI) as shown in Figure 1. Then, masked variables were aggregated to monthly means for each pixel to finally be aggregated to a single monthly time-series containing the spatial average of all pixels in the mask. The mean annual cycle was computed by averaging each of the 12 months of the year during the period starting in January 2003 and ending in February 2020. Annual cycles of standard deviations were computed in the same way. Anomalies for each month in the record were constructed by subtracting the mean of its corresponding month from its actual value. The standardized anomalies correspond to the anomalies normalized by the standard deviation of that month. The same procedure was followed for a seasonal assessment, with the four seasons defined as DJF (December, January, February), MAM (March, April, May), JJA (June, July, August), and SON (September, October, November). The mean values for each season were computed for the period spanning from December 2002 to February 2020. The DJF season includes December of that year plus January and February of the following year. The same analysis was performed over different land-covers within the burned area of the 2019-2020 wildfire events.
VOD data were available up to 2018. Random Forest (RF) was used to complete the VOD time-series for the 2019-2020 period. RF is used widely in different applications from precipitation retrieval to dataset augmentation and is an ensemble learning method for classification, regression, and other tasks that operate by constructing a multitude of decision trees at training time and outputting the class that is the mode of the classes (classification) or mean prediction (regression) of the individual trees [59]. There are three reasons why RF was chosen for time-series augmentation. First, the predictive performance of RF can compete with the best supervised learning algorithms. Second, RF provides a reliable feature importance estimate. Third, RF offers efficient estimates of the test error without incurring the cost of repeated model training associated with cross-validation [60].

Results and Discussion
The total burned area from December 2019 to February 2020 is approximately 137,518 km 2 . Time-series of the total burned area in DJF for the period 2003-2019 shows that the recent wildfire has the highest burned area (Figure 2). In addition, the percentage of the burned areas within the ROI in previous years was investigated and the results show that the highest percentage occurred in 2009 (~6%). This indicates that most of the areas burned in 2019-2020 fire events have not been affected by wildfires since 2003.
Water 2020, 12, x FOR PEER REVIEW 7 of 19 region (5° N-5° S, 120° E-170° W) relative to a single fixed 30-year base period [58]. Simple timeseries of the ONI values were compared with that of FRP pixel counts from 2002 to 2020. Each of the remaining variables were first masked to the extent of the fires during the 2019-2020 fire season (i.e., ROI) as shown in Figure 1. Then, masked variables were aggregated to monthly means for each pixel to finally be aggregated to a single monthly time-series containing the spatial average of all pixels in the mask. The mean annual cycle was computed by averaging each of the 12 months of the year during the period starting in January 2003 and ending in February 2020. Annual cycles of standard deviations were computed in the same way. Anomalies for each month in the record were constructed by subtracting the mean of its corresponding month from its actual value. The standardized anomalies correspond to the anomalies normalized by the standard deviation of that month. The same procedure was followed for a seasonal assessment, with the four seasons defined as DJF (December, January, February), MAM (March, April, May), JJA (June, July, August), and SON (September, October, November). The mean values for each season were computed for the period spanning from December 2002 to February 2020. The DJF season includes December of that year plus January and February of the following year. The same analysis was performed over different land-covers within the burned area of the 2019-2020 wildfire events.
VOD data were available up to 2018. Random Forest (RF) was used to complete the VOD timeseries for the 2019-2020 period. RF is used widely in different applications from precipitation retrieval to dataset augmentation and is an ensemble learning method for classification, regression, and other tasks that operate by constructing a multitude of decision trees at training time and outputting the class that is the mode of the classes (classification) or mean prediction (regression) of the individual trees [59]. There are three reasons why RF was chosen for time-series augmentation. First, the predictive performance of RF can compete with the best supervised learning algorithms. Second, RF provides a reliable feature importance estimate. Third, RF offers efficient estimates of the test error without incurring the cost of repeated model training associated with cross-validation [60].

Results and Discussion
The total burned area from December 2019 to February 2020 is approximately 137,518 km 2 . Timeseries of the total burned area in DJF for the period 2003-2019 shows that the recent wildfire has the highest burned area (Figure 2). In addition, the percentage of the burned areas within the ROI in previous years was investigated and the results show that the highest percentage occurred in 2009 (~6%). This indicates that most of the areas burned in 2019-2020 fire events have not been affected by wildfires since 2003.   year, the burned area includes December of that year plus January and February of the following year. Values plotted were retrieved from MODIS MCD64A1 (collection 6) product.  During SON, which is Australia's fire season, the fire pixel count anomaly does not change significantly in the past 18 years, as shown in Figure 4    During SON, which is Australia's fire season, the fire pixel count anomaly does not change significantly in the past 18 years, as shown in Figure 4 (left). The forested areas indicate large fires in 2013 and 2019. Otherwise, the annual fire events appear to be within its climatological mean as the standardized anomaly of fire pixel count (FPC SDAnom) is close to zero. On the other hand, during DJF (Figure 4 right), there is a significant increase of FPC SDAnom in forested areas in 2019-2020, indicating a significant increase in burned area that agrees with the EOF analysis. year, the burned area includes December of that year plus January and February of the following year. Values plotted were retrieved from MODIS MCD64A1 (collection 6) product.  During SON, which is Australia's fire season, the fire pixel count anomaly does not change significantly in the past 18 years, as shown in Figure 4     higher than the median and the FPCs in SON 2019-2020 are close to the 25 percentile. This statistical data shows that the wildfires in DJF 2019-2020 were unprecedented events that had not occurred since 2002. In forested burned areas shown in Figure 5 (right), the FPC exhibits a worse situation. In DJF 2019-2020, the total FPC is beyond what has occurred in the past. In December and January, the total FPC is around 20,000, which is an outlier in the wildfire climatology of the months. In February, the total FPC is around 3,000, which is also an outlier of the month. Therefore, there is a sharp increase in FPC in a way that many forested areas that had never seen wildfires during the summer season in almost two decades were burned during DJF 2019-2020. This significant increase agrees with the EOF analysis shown in Figure 3 and FPC SDAnom in Figure 4.
Water 2020, 12, x FOR PEER REVIEW 9 of 19 Figure 5 shows the monthly boxplots of FPC over the total burned area and the forested subset from 2002 to 2020. In Figure 5 (left), the FPCs of the total burned area in the 2019-2020 wildfires, marked by black dots, are greater than the 75 percentile, although the FPC in February is only slightly higher than the median and the FPCs in SON 2019-2020 are close to the 25 percentile. This statistical data shows that the wildfires in DJF 2019-2020 were unprecedented events that had not occurred since 2002. In forested burned areas shown in Figure 5 (right), the FPC exhibits a worse situation. In DJF 2019-2020, the total FPC is beyond what has occurred in the past. In December and January, the total FPC is around 20,000, which is an outlier in the wildfire climatology of the months. In February, the total FPC is around 3,000, which is also an outlier of the month. Therefore, there is a sharp increase in FPC in a way that many forested areas that had never seen wildfires during the summer season in almost two decades were burned during DJF 2019-2020. This significant increase agrees with the EOF analysis shown in Figure 3 and FPC SDAnom in Figure 4.  Figure 6. The correlation is true from 2002 to 2011. Dry conditions during El Niño events (positive ONI) associated with suppression of rain cloud development in the anomalously cold sea surface temperature over the West Pacific cause droughts and thus trigger wildfires across Australia as well as its neighboring maritime continent. On the other hand, wet conditions during La Niña events (negative ONI) suppress the wildfires because the soil is generally wet from precipitation. However, in the years following the 2011 La Niña event, the correlation between ONI and the FRP anomalies in DJF appears to be negative. The El Niño events are not ensued by the increased wildfires and the La Niña events do not result in decreasing wildfires. Likely, ENSO is not the only driver of the Australian wildfires. A study of megadrought (2010-2019) in Chile by Garreaud et al. [61] indicates that the El Niño events, which had been believed to be the driver of the wet condition in Chile, had little influence in maintaining the precipitation deficit during the megadrought period, suggesting that other factors contribute to the hydrological circulation around the Pacific.  Figure 6. The correlation is true from 2002 to 2011. Dry conditions during El Niño events (positive ONI) associated with suppression of rain cloud development in the anomalously cold sea surface temperature over the West Pacific cause droughts and thus trigger wildfires across Australia as well as its neighboring maritime continent. On the other hand, wet conditions during La Niña events (negative ONI) suppress the wildfires because the soil is generally wet from precipitation. However, in the years following the 2011 La Niña event, the correlation between ONI and the FRP anomalies in DJF appears to be negative. The El Niño events are not ensued by the increased wildfires and the La Niña events do not result in decreasing wildfires. Likely, ENSO is not the only driver of the Australian wildfires. A study of megadrought (2010-2019) in Chile by Garreaud et al. [61] indicates that the El Niño events, which had been believed to be the driver of the wet condition in Chile, had little influence in maintaining the precipitation deficit during the megadrought period, suggesting that other factors contribute to the hydrological circulation around the Pacific.
Climatological change of the geopotential height can create a favorable condition for fire events. Persistent high-pressure systems trigger high temperatures and suppress evapotranspiration as well as a convective cloud organization. High wind near the surface may contribute to the severity, and enlarge the extent of the fires, especially in dry environments. The regression map of the dominant mode (Figure 7) is the product of the 500 hPa geopotential height anomaly and the vector variance of PC1 and shows the statistically significant pattern of the geopotential height. The second dominant mode (PC2) is also plotted. The EOF analysis shows that the high-pressure system was the dominant mode over southeast Australia in DJF both in 2002-2019 and 2002-2020. This high-pressure system may have contributed to the low precipitation and drier surface condition that is favorable for fire occurrence. The more intense high-pressure system in the latter period suggests that the environment in DJF 2019-2020 was more favorable for wildfire events. Though the intense fire event location was more to the south of the high-pressure system center, both the FRP EOF and geopotential EOF indicate that southeast Australia is the most impacted area. Climatological change of the geopotential height can create a favorable condition for fire events. Persistent high-pressure systems trigger high temperatures and suppress evapotranspiration as well as a convective cloud organization. High wind near the surface may contribute to the severity, and enlarge the extent of the fires, especially in dry environments. The regression map of the dominant mode (Figure 7) is the product of the 500 hPa geopotential height anomaly and the vector variance of PC1 and shows the statistically significant pattern of the geopotential height. The second dominant mode (PC2) is also plotted. The EOF analysis shows that the high-pressure system was the dominant mode over southeast Australia in DJF both in 2002-2019 and 2002-2020. This high-pressure system may have contributed to the low precipitation and drier surface condition that is favorable for fire occurrence. The more intense high-pressure system in the latter period suggests that the environment in DJF 2019-2020 was more favorable for wildfire events. Though the intense fire event location was more to the south of the high-pressure system center, both the FRP EOF and geopotential EOF indicate that southeast Australia is the most impacted area.
It was also found that the mean 10-m wind vector in DJF from 2002 to 2020 was southerly over southern Australia ( Figure S2). It means that the wind has transported colder and drier air mass to the continent, which makes the environment more favorable for wildfires during DJF. However, the calculations did not show any notable anomaly in wind speed over southern Australia. The anomaly analysis exhibits only about 0.2 to 0.6 m s −1 increase, relative to an 18-year wind speed average as shown in Figure S2. Soil moisture depends on complex hydrometeorological interactions mainly driven by precipitation, temperature, and evapotranspiration and modulated by the characteristics of the soil. Conversely, SM is the main driver of the vegetation water content; hence, impacting the risk of wildfire events [62]. The response of SM to its forcing variables is smoother than any other variable as SM acts as a reservoir, which also adds memory to the system [63,64]. In Figure 8, the standardized anomalies of SM show a wet period from 2010 to 2012, but after that, a recurrent drier state has been the norm. Since June 2017, only negative anomalies have been recorded with 2019 and the first It was also found that the mean 10-m wind vector in DJF from 2002 to 2020 was southerly over southern Australia ( Figure S2). It means that the wind has transported colder and drier air mass to the continent, which makes the environment more favorable for wildfires during DJF. However, the calculations did not show any notable anomaly in wind speed over southern Australia. The anomaly analysis exhibits only about 0.2 to 0.6 m s −1 increase, relative to an 18-year wind speed average as shown in Figure S2.
Soil moisture depends on complex hydrometeorological interactions mainly driven by precipitation, temperature, and evapotranspiration and modulated by the characteristics of the soil. Conversely, SM is the main driver of the vegetation water content; hence, impacting the risk of wildfire events [62]. The response of SM to its forcing variables is smoother than any other variable as SM acts as a reservoir, which also adds memory to the system [63,64]. In Figure 8, the standardized anomalies of SM show a wet period from 2010 to 2012, but after that, a recurrent drier state has been the norm. Since June 2017, only negative anomalies have been recorded with 2019 and the first months of 2020 being more than two times the standard deviation below the normal. During December 2019, a record anomaly of about −70 mm was recorded. In addition, steady drying was observed from the monthly anomalies time-series since 2011 (Figure 9). The very low amounts of precipitation during December 2019, during a particularly warm growing season of the vegetation, were the main drivers of the record low soil moisture. When the season before the peak of forest fires is analyzed (SON; Figure 8), a rapid decrease in soil moisture during the last three years coinciding with the three low record values for the period of analysis can be noticed. Even a more extreme decrease was seen for the DJF season ( Figure 9). This extremely low soil moisture before and during the current wildfire events is arguably one of the responsible factors for a quick spread of any active fire, as also reported by Farahmand et al. [20,21] over the United States. Only slight differences were observed when the area was limited to that covered by evergreen forests.
Water 2020, 12, x FOR PEER REVIEW 12 of 19 The standardized monthly mean temperature anomalies show an upward trend from 2010 to 2019 ( Figure 9); positive values indicate that the temperature was higher than the means from 2002 to 2019 (i.e., almost 1 °C above the mean for the total area (ROI) and forested subset during the 2019-2020 fire season). In fact, the highest air temperature in Australia was registered in 2019 that may suggest a relationship between rising temperature and wildfire events.
The annual precipitation in 2019 over the study area has been exceptionally low, reaching about 530 mm (mean value for the 2002-2019 period is 750 mm), which may indicate the role of precipitation in the recent fire events. Precipitation decrease in 2019 was intensified in November (25 mm, 62 mm normal) and December (21 mm, 89 mm normal), compared to the previous years (Figures 8 and 9). Figure 8 shows that every month (except March) in 2019 has negative anomalies setting the dry conditions preceding the fire events that occurred in December 2019 and continued until February 2020. The same pattern for precipitation anomalies was observed over burned evergreen forest areas in southeast Australia. Higher precipitation in March 2020 (109 mm, 97 mm normal) may highlight the important role of precipitation in extinguishing fires as well.
Based on Figures 8 and 9, actual evapotranspiration (AET), potential evapotranspiration (PET), and water vapor deficit (VPD) have changed significantly in recent years. There was a decreasing trend in AET, while the PET and VPD increased in recent years. Increasing air temperature and decreasing precipitation, which are the most dominant factors of wildfires, can justify these trends. PET and VPD have increased, mainly due to high air temperature in the region. In contrast, AET has decreased due to less biomass and more water stress on the ground because of wildfires. Changes in the DJF are more notable than SON. It partly shows that the main impacts of wildfires on AET, PET, and VPD have occurred in DJF. pressure deficit (VPD). The plotted values are spatial averages for the total burned area (i.e., ROI; left) and forested subset (right). Monthly anomalies for NDVI, EVI, LAI, and VOD generally show positive anomalies early in 2019 (as well as previous years) changing to negative anomalies in the final months of 2019 ( Figure  10). This may indicate an increase in vegetation early in 2019, followed by vegetation water stress during the months approaching the start of the fire, which suggests that other factors (e.g., drought) are creating prime conditions for increased fuel load in the form of moisture-limited vegetation. For optical indices (i.e., NDVI, EVI, LAI), negative anomaly values in the final months of 2019 are seen in both the forested and total burned areas, reaching more than 5.4 and 6.8 times the standard deviation below the normal for the total burned areas and the forested burned areas, respectively. For these same months, VOD anomalies in the total burned area show slightly increasing values in contrast to the anomalies shown by optical vegetation variables. Because increasing VOD typically is associated with increasing biomass or vegetation water content, these anomaly values indicate a decrease in fire likelihood for the total burned area. However, monthly VOD anomalies approaching the end of 2019 for the forested area are consistent with the decreasing trend seen with optical vegetation indices. Differences in monthly VOD anomalies between the total and forested burned areas suggest high sensitivity of VOD to different land-covers or vegetation types [47]. Seasonal anomalies ( Figure 11  The standardized monthly mean temperature anomalies show an upward trend from 2010 to 2019 ( Figure 9); positive values indicate that the temperature was higher than the means from 2002 to 2019 (i.e., almost 1 • C above the mean for the total area (ROI) and forested subset during the 2019-2020 fire season). In fact, the highest air temperature in Australia was registered in 2019 that may suggest a relationship between rising temperature and wildfire events.
The annual precipitation in 2019 over the study area has been exceptionally low, reaching about 530 mm (mean value for the 2002-2019 period is 750 mm), which may indicate the role of precipitation in the recent fire events. Precipitation decrease in 2019 was intensified in November (25 mm, 62 mm normal) and December (21 mm, 89 mm normal), compared to the previous years (Figures 8 and 9). Figure 8 shows that every month (except March) in 2019 has negative anomalies setting the dry conditions preceding the fire events that occurred in December 2019 and continued until February 2020. The same pattern for precipitation anomalies was observed over burned evergreen forest areas in southeast Australia. Higher precipitation in March 2020 (109 mm, 97 mm normal) may highlight the important role of precipitation in extinguishing fires as well.
Based on Figures 8 and 9, actual evapotranspiration (AET), potential evapotranspiration (PET), and water vapor deficit (VPD) have changed significantly in recent years. There was a decreasing trend in AET, while the PET and VPD increased in recent years. Increasing air temperature and decreasing precipitation, which are the most dominant factors of wildfires, can justify these trends. PET and VPD have increased, mainly due to high air temperature in the region. In contrast, AET has decreased due to less biomass and more water stress on the ground because of wildfires. Changes in the DJF are more notable than SON. It partly shows that the main impacts of wildfires on AET, PET, and VPD have occurred in DJF.
Monthly anomalies for NDVI, EVI, LAI, and VOD generally show positive anomalies early in 2019 (as well as previous years) changing to negative anomalies in the final months of 2019 ( Figure 10). This may indicate an increase in vegetation early in 2019, followed by vegetation water stress during the months approaching the start of the fire, which suggests that other factors (e.g., drought) are creating prime conditions for increased fuel load in the form of moisture-limited vegetation. For optical indices (i.e., NDVI, EVI, LAI), negative anomaly values in the final months of 2019 are seen in both the forested and total burned areas, reaching more than 5.4 and 6.8 times the standard deviation below the normal for the total burned areas and the forested burned areas, respectively. For these same months, VOD anomalies in the total burned area show slightly increasing values in contrast to the anomalies shown by optical vegetation variables. Because increasing VOD typically is associated with increasing biomass or vegetation water content, these anomaly values indicate a decrease in fire likelihood for the total burned area. However, monthly VOD anomalies approaching the end of 2019 for the forested area are consistent with the decreasing trend seen with optical vegetation indices. Differences in monthly VOD anomalies between the total and forested burned areas suggest high sensitivity of VOD to different land-covers or vegetation types [47]. Seasonal anomalies ( Figure 11

Conclusions
Linking the climate setting to the wildfires can provide insightful information to predict wildfires and understand their future behavior. Burned areas over Australia have been steadily increasing in the past two decades. In this context, the 2019-2020 fire season has set a new high record. The severity of these wildfires was high enough to substantially change the dominant mode in the EOF analysis of the fire radiative power. The large-scale analysis was consistent with the expectations, showing that the ENSO modulation in the precipitation patterns over Australia has a profound impact on the dry/wet conditions preceding the fire seasons. In the same way, the geopotential height anomalies showed a dominant mode with a high-pressure system over southeastern Australia, which was enhanced during the 2019-2020 wildfire events.
Consistent behavior of the examined variables highlights the role of the steady warming and drying as one of the main drivers of the severity and large extent of the fires in the 2019-2020 season. The 2019-2020 DJF season also showed a record high temperature anomaly with almost 1 • C above the normal since 2003. The elevated temperature led to an increased evapotranspiration demand as seen through the potential evapotranspiration and vapor pressure deficit positive anomalies. In contrast, this excess of demand was not satisfied due to the lack of precipitation in the previous years, yielding negative anomalies of actual evapotranspiration. This resulted in a very severe deficit of soil moisture, which led to an increase in the dryness of the vegetation and setting the conditions for faster expansion of the wildfires. Vegetations showed large greenness well before the 2019-2020 event, rapidly decreased after October 2019. This suggests the existence of large amounts of vegetation that rapidly became a great source of dry fuel, feeding the wildfires, and enhancing their severity. The pixel count of burned areas in the forested portion was 20 times larger than the standard deviation of burned pixels from previous years, indicating that the 2019-2020 event was very abnormal.
In this study, several climate variables and satellite-driven vegetation indices and their anomalies were investigated to assess whether climate data reveal a fire season response to climatic variability. The combination of remote sensing and in situ observations enabled us to investigate various important variables. However, there remain other useful variables (e.g., crown bulk density and canopy base height) that were not considered in this study mainly due to the limited access to reliable data.