The Spatiotemporal Characteristics of Wildfires across Australia and Their Connections to Extreme Climate Based on a Combined Hydrological Drought Index

Abstract

With the frequent occurrence of extreme climates around the world, the frequency of regional wildfires is also on the rise, which poses a serious threat to the safety of human life, property, and regional ecosystems. To investigate the role of extreme climates in the occurrence and spread of wildfires, we combined precipitation, evapotranspiration, soil moisture (SM), maximum temperature (MT), relative humidity, plant canopy water, vapor pressure deficit, and a combined hydrological drought index based on six Gravity Recovery and Climate Experiment (GRACE) and its follow-on (GRACE-FO) products to study the relationship between climate change and wildfires across Australia between 2003 and 2020. The results show that Australia’s wildfires are mainly concentrated in the northern region, with a small number being distributed along the southeastern coast. The high burned months are September (2.5941 × 10⁶ ha), October (4.9939 × 10⁶ ha), and November (3.8781 × 10⁶ ha), while the years with a larger burned area are 2011 (79.95 × 10⁶ ha) and 2012 (78.33 × 10⁶ ha) during the study period. On a seasonal scale, the terrestrial water storage change and the hydrometeorological factors have the strong correlations with burned area, while for only the drought index, SM and MT are strongly related to burned area on an interannual scale. By comparing the data between the high burned and normal years, the impact of droughts on wildfires is achieved through two aspects: (1) the creation of a dry atmospheric environment, and (2) the accumulation of natural combustibles. Extreme climates affect wildfires through the occurrence of droughts. Among them, the El Niño–Southern Oscillation has the greatest impact on drought in Australia, followed by the Pacific Decadal Oscillation and the Indian Ocean Dipole (correlation coefficients are −0.33, −0.31, and −0.23, respectively), but there is little difference among the three. The proposed hydrological drought index in our study has the potential to provide an early warning of regional wildfires. Our results have a certain reference significance for comprehensively understanding the impact mechanism of extreme climates on regional wildfires and for establishing an early warning system for regional wildfires.

1. Introduction

Wildfire is a devastating natural disaster. Due to global warming, the likelihood of extreme weather has increased, and continuous high local temperature and drought have greatly increased the risk of wildfires [1]. Australia is one of the most fire-prone regions in the world, and its wildfire situation is complex [2,3]. The 2019/2020 Australian wildfire is believed to be the largest wildfire in the Australian history, which resulted in serious socioeconomic, environmental, and ecological impacts [4]. The direct consequence of this wildfire was the loss of 33 lives, the destruction of over 3000 houses, and more than 800 hectares of burned land [5]. The heavy smoke released by the wildfire lasted for several weeks, covering many major cities in Australia and bringing adverse effects on people’s health [6]. The long-term impacts on the local environment and ecology are difficult to measure. This wildfire generated a heated debate on the impact of climate change on wildfire risk. Scientific studies indicate that climate change does raise the risk of wildfires in Australia and the rest of the world [7,8,9]. Therefore, discussing and analyzing the role of climate change in wildfires is helpful to understand the formation mechanism of wildfires and to realize the early warning of regional wildfires.

The effects of climate change on wildfire are reflected on different time scales. On the long-term scale, climate change will alter regional climate, resulting in a trend of local aridification, which will lead to the changes in local vegetation types; on the short-term scale, drastic changes in climate factors will create environmental conditions that are prone to wildfires and provide sufficient combustibles for wildfires [10]. The scientific research demonstrating the probability of wildfires shows an increasing trend in several regions around the globe. This trend may be due to climate change [11]. Williams et al. [12] quantified the possible impact of climate change on wildfires based on the weather changes during Australian wildfires and the McArthur Forest Fire Danger Index, and the investigation showed that global warming does have a significant impact on the rising risk of wildfires in Australia. Westerling et al. [13] revealed the response of wildfires to climate change based on the gridded monthly fire-and-burned-area dataset from the U.S. Forest Service, the Bureau of Land Management, the National Park Service, and the Bureau of Indian Fire. Liu et al. [14] indicated that global wildfire risk is on the rise due to greenhouse effect according to the General Circulation Models. Researchers have usually analyzed the relationship between wildfires and hydroclimate changes mainly by observing and studying the abnormal signals in hydrometeorological variances before and during wildfires. Most wildfires are preceded by high temperatures, low relative humidity, and strong winds, and these factors are conducive to the occurrence and spread of wildfires [15,16,17]; this view has become the consensus. In addition to creating an external environment conducive to wildfires, climate change can also provide sufficient fuel for wildfires by affecting plant production [18,19,20]. Soil moisture (SM) has an important impact on the water content in vegetation. Long-duration droughts cause severe SM deficit, decrease vegetation moisture, and increase the number of combustibles [21]. Chaparro et al. [22] used SM as a key indicator to measure the impact of climate change on wildfires and combined it with temperature data to predict wildfires in the Iberian Peninsula. Krueger et al. [23] pointed out that SM is closely related to the probability of wildfires in Oklahoma, USA, and this impact has a significant seasonality. Jensen et al. [24] demonstrated that there is a strong correlation between SM and wildfires in the United States, and this result can be used to assess the likelihood of wildfire occurrence.

From the above analysis, we can see that hydrological variances have a close relationship with the occurrence and spread of wildfires. Although ground meteorological stations can collect high-precision point data, they suffer from insufficient coverage and uneven distribution, especially in sparsely populated regions. Therefore, they cannot meet the needs of large-scale regional research [25,26]. Remote sensing technology can obtain high-resolution gridded data in a large-scale region, but it is limited to surface water and shallow SM [27]. Therefore, we urgently need a new technique that can monitor all types of terrestrial water storage change (TWSC) for a comprehensive assessment of the impact of climate change on terrestrial hydrology.

The Gravity Recovery and Climate Experiment (GRACE) mission, which has been implemented since 2002, can collect high-precision Earth’s time-variance gravity field information. This information can be translated into the total TWSC [28]. Therefore, the GRACE data have unique advantages in assessing local hydrological conditions. Many scholars have applied the GRACE solutions to detect and quantify regional hydrological drought, for example, in the Amazon River basin [29,30], southwest China [31], Australia [32], and United States [33]. Due to the great success of the GRACE mission, its follow-on (GRACE-FO) satellites were launched in May 2018. So far, these two missions have provided nearly 20 years of observational data, which facilitates the study of the long-term changes of regional climate and hydrology. Some scholars have begun to use the GRACE data to discuss the connection between hydroclimate variances and regional wildfires. Chen et al. [34] applied the GRACE TWSC data to study the climate and hydrological anomalies before and during wildfires in the Amazon River basin from 2002 to 2011. The results suggested that the GRACE TWSC data hold the potential for wildfire forecasting in the basin. Jensen et al. [24] estimated the monthly SM data based on the GRACE observations in the United States during the period between 2003 and 2013, compared them with the wildfire data of the same period, and found that the wet and dry states of soil have a significant correlation with wildfire risk. Farahmand et al. [35] used the SM data derived from the GRACE satellites and vapor pressure deficit (VPD) from the Atmospheric Infrared Sounder mission to predict fire-prone months and regions and explained that the hydrological information plays a very important role in improving early warning predictions of wildfire. Cui et al. [36] evaluated the abnormal changes in various hydroclimatic factors based on the GRACE TWSC data before forest fires in Yunnan, Southwest China, during the period between 2003 and 2016. The results indicated that the GRACE TWSC data are more sensitive to the hydrological environment before a forest fire than the PPT. However, the above studies only used the GRACE data and did not combine the GRACE and the GRACE-FO data to form a longer-time series of observations for monitoring and discussing the hydroclimate conditions before and during wildfires. Compared to the TWSC, the drought index can better reflect the wet and dry conditions of a region [37]. However, to best of our knowledge, there is currently no literature focusing on the connection between the GRACE/GRACE-FO-based drought index (GRACE-DSI) and regional wildfires.

Therefore, our research goal is to combine six GRACE and GRACE-FO solutions to construct a 19-year drought index time series, and this drought index time series is used to study the connection between extreme climates and wildfires.

2. Study Area

Australia is located in the Oceania, and its area is 7.69 million km² (Figure 1). It is the smallest one among the six continents in the world. There are a variety of climate conditions spread across Mainland Australia, ranging from humid tropical condition to dry temperate continental condition [38]. Australia has a hot and dry climate and a flat terrain. The daily average minimum temperature is 280.15 K, the maximum temperature is 288.15 K, and the average annual PPT is 355 mm. There are many hills and few mountains, and the three species mainly belong to eucalyptus, which contains a lot of oil. It is affected by many factors. Australia is one of the regions that are most prone to forest fires in the world. Due to less precipitation (PPT), severe drought events are more prone to occur in Australia. The climate over Australia is more susceptible to the changes in the temperature and pressure of the surrounding oceans. Therefore, the influences of the El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD) on Australia are very significant [39].

3. Data and Method

3.1. Data

3.1.1. GRACE/GRACE-FO Data

The GRACE/GRACE-FO RL06 spherical harmonic (SH) solutions, which are used to estimate the monthly gridded TWSC data, are from the Center for Space Research at the University of Texas at Austin (CSR), the Helmholtz Centre Potsdam–German Research Centre for Geosciences (GFZ), the Jet Propulsion Laboratory (JPL), and the Institute of Geodesy at the Graz University of Technology (ITSG), respectively. We need to preprocess the SH solutions to improve the accuracy of the gridded TWSC data. Firstly, the degree-1 and C20 coefficients of the SH solutions were replaced with the results from Swenson et al. [40], and the C20 coefficients were estimated by using satellite laser ranging [41]. Then, the combined filtering (300 km Gaussian filter and P3M6 polynomial filter) was used to weaken the influence of high-frequency and correlated errors in the SH solutions [42]. Finally, we used the scale factor approach to restore the lost signals due to order truncation and filter smoothing [43]. Aside from the SH solutions, we could also extract the monthly gridded TWSC data directly from the Mascon solutions without any additional data processing. The Mascon solutions were provided by the CSR and the JPL for our study. Therefore, we used the four SH solutions and two Mascon solutions to estimate the monthly TWSC gridded data in Mainland Australia from January 2003 to December 2020. Because the GRACE and GRACE-FO data are the same type of data, we unified the GRACE and GRACE-FO data as the GRACE data. For convenience, the four GRACE SH solutions and the two Mascon solutions were termed as CSR-SH, GFZ-SH, JPL-SH, ITSG-SH, CSR-M, and JPL-M.

Due to an 11-month data gap between the GRACE and GRACE-FO missions, we used the reconstructed TWSC data to fill this gap. The data were estimated based on PPT, land surface temperature, sea surface temperature, and climate index data by using principal component analysis, least squares regression, and multiple linear regression methods. The specific technical details of the PCA-LS-MLR method can be found in Li et al. [44].

3.1.2. Burned Area Data

The monthly burned-area gridded data from the Climate Change Initiative project of the European Space Agency are based on surface reflectance in the near infrared band obtained from the Moderate-Resolution Imaging Spectroradiometer (MODIS) instrument on board of the Terra satellite, as well as active fire information from the same sensor of the Terra and Aqua satellites [45,46]. The gridded data with a spatial resolution of 0.25° include the sum of the burned area, the standard error, the fraction of the burned area, the fraction of the observed area, the number of patches, and the sum of the burned area for each land cover class.

3.1.3. In Situ Climate Data

In our study, we obtained the PPT, maximum temperature (MT), relative humidity (RH), and VPD gridded data from the Scientific Information for Land Owners (SILO). This dataset mainly includes Australian climate data from 1889 to the present, which provides daily datasets for a range of climate variances. The SILO products include the station time series and spatial gridded data, which are based on the observation data from the Bureau of Meteorology database. For point data, interpolated or derived values are used where observations are missing. The gridded data are spatially interpolated from the observations. Since the temporal resolution of the data in our study is monthly, we averaged the SILO products on a monthly basis to maintain uniformity.

3.1.4. Other Hydrometeorological Data

The ET gridded data, with a spatial resolution of 0.25° × 0.25°, were obtained from the Global Land Evaporation Amsterdam Model (GLEAM) v3.6a. The GLEAM v3.6a provides various components of land evaporation, including open-water evaporation, bared-soil evaporation, transpiration, interception loss, sublimation, potential evaporation, surface and root-zone soil moisture, and evaporative stress conditions [47,48]. The SM and plant canopy water (PCW) gridded data were obtained from the Global Land Data Assimilation System (GLDAS) 2.1 model. This model is a hydrological model generated by the Goddard Space Flight Center at the National Aeronautics and Space Administration and the National Center for Environmental Prediction at the National Oceanic and Atmospheric [49]. The monthly SM data with a spatial resolution of 1° × 1° provide the sum of these 4-layer SM data in the model. In our study, the monthly ET, SM, and PCW data are from 2003 to 2020.

3.1.5. Standardized Precipitation Evapotranspiration Index Data

The Standardized Precipitation Evapotranspiration Index (SPEI) is estimated based on the cumulative sum of difference between PPT and potential ET on different time scales, which calculation approach is the same as the Standardized Precipitation Index (SPI) [50]. The SPEI has three different time scales: 3 months, 6 months, and 12 months. The 3-month scale is associated with variations in soil moisture and represents agricultural drought; the 6-month scale is associated with variations in streamflow and represents hydrological drought; and the 12-month scale is associated with variations in groundwater storage and represents hydrogeological drought [51]. The monthly SPEI gridded data span the period from January 2003 to December 2018, which spatial resolution is 1° × 1° in our study.

3.1.6. Extreme Climate Index Data

The ENSO, the Pacific Decadal Oscillation (POD), and the IOD are extreme climate phenomena caused by an anomaly in sea surface temperature difference in the Pacific and Indian Oceans, respectively [52,53,54]. Since Australia is sandwiched between the Pacific and Indian Oceans, it is significantly affected by ENSO, the POD, and the IOD. Therefore, we must consider the impact of the ENSO and the IOD on wildfires in Mainland Australia. The monthly ENSO, POD, and Indian Ocean Dipole Mode Index (DMI) was provided by the National Oceanic and Atmospheric Administration. The DMI is defined as the difference in the average sea surface temperature anomaly between the Tropical Western Indian Ocean and the Equatorial Southeast Indian Ocean [55]. The above index data in our study are from 2003 to 2020.

3.2. Method

3.2.1. Data Fusion

Due to the discrepancies in the GRACE solutions from different institutions, it may lead to inconsistent TWSC results, which can adversely affect our results [56]. Therefore, we applied the generalized three-cornered hat (GTCH) and least squares method to fuse the TWSC results from the six GRACE solutions to improve the reliability of our results. Among them, the GTCH method was used to determine the uncertainties of different GRACE solutions, which advantage is that it does not need any prior information [57,58]. The technical details about the GTCH method can be found in Long et al. [59] and Cui et al. [36].

We fused the TWSC results from the six GRACE solutions by using the least squares method. The expression is as follows:

$$X={\displaystyle \sum _{i=1}^N{p}_i•X_i\left(i=1,2,\cdots ,N\right)}$$

(1)

where $X$ is the fused result, $X_i$ is the $i$th dataset, and $p_i$ is the corresponding weight. The weights are estimated based on the relative variances of each dataset:

$$p_i=\frac{1/{\sigma }_{ii}}{{\displaystyle \sum _{n=1}^{N}1/{\sigma }_{nn}}}$$

(2)

where ${\sigma }_{ii}\left(i=1,2,\cdots ,N\right)$ is the relative variance of the $i$th dataset derived from the GTCH method. In our study, the different datasets on each grid point were fused.

3.2.2. GRACE-DSI

To obtain the GRACE-DSI data, the time series of the GRACE TWSC results were standardized [60].

$$GRACE-DS{I}_{i,j}=\frac{TWS{C}_{i,j}-TWS{C}_j^{mean}}{{\sigma }_j}$$

(3)

where $TWS{C}_{i,j}$ is TWSC in the $i$ year $j$ month. $TWS{C}_j^{mean}$ and ${\sigma }_i$ are the average and standard deviation of TWSC in $j$ month, respectively. The missing data in the time series of the GRACE TWSC results were filled based on the cubic spline interpolation. The drought events can be classified according to the GRACE-DSI value (

Table 1. GRACE-DSI drought grade classification.

Type	GRACE-DSI	Type	GRACE-DSI
Exceptional Drought	$\le -2.0$	Moderate Drought	−1.3~−0.8
Extreme Drought	−2.0~−1.6	Light Drought	−0.8~−0.5
Severe Drought	−1.6~−1.3	No Drought	$\ge -0.5$

) [31,61].

3.2.3. Composite Analysis

In climate change research, composite analysis is used to determine the significant characteristics of special years. In our study, the time series of burned area was standardized. Then, the special years were fixed according to the rule of greater than 0.5 times or less than 0.5 times of the standard deviation based on the standardized time series of burned area, which were marked as high- or low-fire years, respectively. We can compare the average of meteorological and hydrological variances during these special years with the ones of all years to discuss the performance of these variances during high- and low-fire years.

3.2.4. The Correlation Analysis and Delay Months

In hydrology research, the response between two climatic variances is not immediate and usually delayed for a certain period of time. Therefore, to discuss the relationship between two climatic variances, we do not calculate their correlation coefficient (CC), but the maximum CC and delay months. There are two independent time series, $x_1$ and $x_2$, and the expressions of their CC $\rho \left(\tau \right)$ and delay months $\tau$ are as follows [62,63]:

$$\rho \left(\tau \right)=\frac{{\sigma }_{12}\left(\tau \right)}{\sqrt{{\sigma }_{11}{\sigma }_{12}}}$$

(4)

where ${\sigma }_{11}$ and ${\sigma }_{22}$ are the variances of $x_1$ and $x_2$, respectively, and ${\sigma }_{12}$ is the covariance of $x_1$ and $x_2$. When $\left|\rho \left(\tau \right)\right|$ is maximum ($\left|\rho \left(\tau \right)\right|\le 1$), $\tau$ is the corresponding delay months ($\left|\tau \right|\le 12$).

3.2.5. Nash–Sutcliffe Efficiency (NSE)

The Nash–Sutcliffe coefficient [64] has been widely applied to evaluate the performance of simulated time series against observation series in hydrology research. The specific expressions can be found in Ref. [65]. $X^{mean}$ is the mean value of the observation series, and $n$ is the number of observations. The value range $NSE$ is from $-\infty$ to 1. $NSE=1$ means that the two-time series fits perfectly; $0<NSE\le 1$ suggests that the performance of the simulated series is acceptable; and $NSE\le 0$ explains that the simulated series is worse [66]. In this study, we used the NSE as an indicator to evaluate the relationship between the GRACE-DSI (observation series) and the traditional drought indices (simulated series).

3.2.6. Standard Precipitation Index (SPI)

The SPI is a meteorological drought index, which is based on long-term PPT observations and is used to characterize the probability of PPT over a period of time [67]. The SPI is widely used in regional drought monitoring due to its simple calculation and reliable results. The probability density function of gamma distribution for long-term PPT data is first calculated and then converted to a normal distribution to have a SPI mean value at 0. The specific expressions can be found in Ref. [68]. According to different research purposes, the SPI of different time scales can be calculated. The SPI-3, SPI-6, and SPI-12 represent meteorological drought, agricultural drought, and hydrological drought, respectively.

4. Results

4.1. GRACE-DSI Construction

Figure 2 shows the spatial distribution of uncertainties in the TWSC results from the six GRACE solutions in Australia. All the results related to uncertainties show a similar spatial distribution, that is, larger uncertainties are mainly concentrated in northern Australia. In this region, the uncertainties of the SH solutions (less than 4 cm) are smaller than those of the Mascon solutions (between 5 cm and 8 cm). The JPL-SH and the ITSG-SH show larger uncertainties (4~4.5 cm) than the CSR-SH and the GFZ-SH, while the red region (greater than 6.5 cm) in Figure 2f (JPL-M) is larger than the one in Figure 2e (CSR-M). We sorted the grid value of uncertainties in the six TWSC results in Australia in ascending order and applied the median one to represent the overall uncertainty level of these TWSC results. According to the ascending order of uncertainties in the TWSC results, the arrangement of the six GRACE solutions is the CSR-SH (17.5668 mm), the ITSG-SH (18.2371 mm), the GFZ-SH (19.7522 mm), the JPL-SH (19.9046 mm), the CSR-M (21.1316 mm), and the JPL-M (26.2787 mm).

To improve the accuracy of the TWSC results, we fused the TWSC results from the six GRACE solutions. The spatial distribution map of uncertainties in the fused results is plotted in Figure 3. All the uncertainties are less than 3 cm in the study region, and the maximum value appears in north Australia (greater than 2 cm). When comparing Figure 2 and Figure 3, the uncertainties in the fused results are significantly smaller than those of any single solution. In addition, the level of uncertainties in the fused results (9.4597 mm) is much smaller than those of the six GRACE solutions.

We compared the time series of the TWSC results to the six GRACE solutions and the fused results in Australia during the period between 2003 and 2020 (Figure 4). These time series have significant seasonal variation and the same change trend. The fused results have the greatest correlation coefficient with the CSR-SH (0.9824), followed by the ITSG-SH (0.9786), the GFZ-SH (0.9778), the JPL-SH (0.9775), the CSR-M (0.9722), and the JPL-M (0.9716). This indicates that the fused results have a good consistency with the TWSC results from the six single solutions. When comparing Figure 3 and Figure 4, the fused TWSC results can better ensure the reliability of our study. Thus, we used them to construct the GRACE-DSI. Figure 5 shows the time series of the GRACE-DSI and the traditional drought indices (SPEI-3, SPEI-6, SPEI-12, SPI-3, SPI-6, and SPI-12). The seven drought indices have a similar change trend. Except for the SPEI-3, the GRACE-DSI has a strong correlation with the other five drought indices (the correlation coefficients are greater than 0.5). The GRACE-DSI has the strongest correlation coefficient with the SPI-12 (0.74), followed by the SPEI-12 (0.68), the SPI-6 (0.65), the SPEI-6 (0.54), the SPI-3 (0.52), and the SPEI-3 (0.40). In addition, the NSEs between the GRACE-DSI and the SPEI-3, the SPEI-6, the SPEI-12, the SPI-3, the SPI-6, and the SPI-12 are 0.01, 0.17, 0.32, 0.07, 0.30, and 0.49, respectively, which indicates that the GRACE-DSI is within the acceptable range. The above results explain that the GRACE-DSI can be applied to monitor droughts in Australia.

4.2. Spatiotemporal Distribution of Burned Area

Figure 6 shows the spatial distribution of monthly average burned area from 2003 to 2020 in Australia. In terms of spatial distribution, wildfires in Australia are mainly concentrated in the northern region, and there are small wildfires in the southeast region. In terms of temporal distribution, the minimum burned area occurs in autumn (March, 6.94 million ha). From March onwards, the burned area gradually increases until October. The burned area hits the largest peak in October (spring, 152.33 million ha). Subsequently, the burned area begins to gradually decrease from October to March. According to Figure 7, the year with the largest burned area is 2011 (78.95 million ha), followed by 2012 (78.33 million ha), while the year with the smallest one is 2010 (18.33 million ha), followed by 2020 (18.40 million ha). According to the composite analysis (Section 3.2.3), the high burned years are 2004, 2006, 2011, and 2012, respectively, while the low burned years are 2003, 2005, 2008, 2010, 2013, 2016, 2019, and 2020 (Figure 8).

4.3. The Connection between Hydrometeorological Factors and Wildfires

4.3.1. On a Seasonal Scale

The seasonal variation in burned area, TWSC, and the seven hydrometeorological factors in Australia during 2003–2020 were calculated by subtracting the average of the corresponding months for all years from the data for each month (Figure 9). It shows that TWSC, burned area, and hydrometeorological factors have significant seasonal change. In Figure 9a, both the maximum of TWSC (3.4 cm) and the minimum of burned area (−3.08 × 10⁻⁶ ha) appear in March, while the minimum of TWSC occurs in November (−1.2 cm), which is only one month away from the month of maximum of burned area (4.99 × 10⁻⁶ ha, October). Additionally, TWSC shows a strong negative correlation (CC = −0.76) with burned area. The maximum and minimum of PPT appear in January (5.4 cm) and August (−2.4 cm), respectively (Figure 9b). It shows that PPT has a two-month delay with burned area (CC = −0.74). ET has a similar seasonal distribution to PPT (Figure 9c), that is, the maximum and minimum values occur in January (2.2 cm) and July (−1.5 cm), respectively. This is attributed to the fact that water consumed by ET mainly comes from PPT. The response of the burned area to ET lags by two months, and its CC is −0.70. In Figure 9d, when RH is positive, the burned area is negative, and vice versa. This suggests that low RH promotes the occurrence of wildfire (CC = −0.77). When comparing Figure 9e,f, VPD and MT have the same seasonal change, that is, the maximum and minimum values occur in January (summer) and June (winter), which is consistent with the seasonal distribution in Australia. However, VPD and MT have no significant relationship with burned area on a seasonal scale. When comparing Figure 10a,g, SM and TWSC share the same seasonal change, that is, the maximum and minimum values appear in March and November. However, SM has a stronger correlation (CC = −0.89) with burned area than TWSC. In Figure 9h, the maximum and minimum of PCW appear in January (2.93 × 10⁻³ cm) and September (−1.60 × 10⁻³ cm), respectively. In addition, the CC between PCW and burned area is −0.82, and the response of burned area to PCW is one month.

According to the above results, a complete interaction process between TWSC, hydrometeorological factors, and wildfire in Australia is demonstrated on a seasonal scale. From January (summer) onwards, PPT and MT begin to decrease or drop from their peaks (PPT and MT, CC = 0.83), which leads to the weakening of ET and the decrease in PCW (PPT and ET, CC = 0.95; PPT and PCW, CC = 0.98; and MT and ET, CC = 0.94). However, due to the replenishment of SM, the decline rate of PCW is not as fast as that of PPT (SM and PCW, CC = 0.77). Although ET weakens, it continues to inject water vapor into the air, which causes a decrease in VPD (ET and VPD, CC = 0.86). With PPT continuing to decline, so do TWSC and SM after peaking in March. At this time, the burned area reaches the minimum. From June onwards, with TM at its bottom and continuing to grow, RH and VPD reach the maximum and minimum, respectively, which indicates that air moisture reaches its peak (TM and RH, CC = −0.41; VPD and TM, CC = 0.98). This leads to a small increase in PCW. At this time, air moisture begins to decrease (a decrease in RH and an increase in VPD). However, SM begins to decrease and air moisture continues to decline (RH and SM, CC = 0.77). The difficulty of obtaining water for vegetation increases, and PCW shows a downward trend. From September onwards, TWSC, SM, PCW, and RH successively reach the minimum, and VPD also reaches the maximum. The reduction in SM and air moisture creates a powerful external environment for wildfire to occur (burned area and RH, CC = −0.77), and the reduction in PCW provides a large number of combustibles for wildfire (burned area and PCW, CC = −0.65). Therefore, from August onwards, the burned area begins to increase. In October, the burned area reaches its peak. With the increase in PPT from September onwards, the external environment begins to change from dry to wet, and PCW increases, which causes a decrease in the burned area. This continues until January of the following year when a new cycle begins.

4.3.2. On an Interannual Scale

Figure 10 shows the interannual variation of the GRACE-DSI and the anomaly of burned area and hydrometeorological factors in Australia during 2003~2020. Among the GRACE-DSI and the seven hydrometeorological factors, burned area is only strongly correlated with the GRACE-DSI (0.66), MT (−0.72), and SM (0.62). It shows that the GRACE-DSI has a closer relationship with burned area in Australia than PPT (0.22). In addition, SM has a strong correlation with burned area, which is due to the close connection between SM and RH, and between VPD and PCW. On an interannual scale, the GRACE-DSI is more sensitive to climate change before a wildfire in Australia.

4.3.3. Performance of GRACE-DSI and Hydrometeorological Factors before Wildfire

From the previous analysis, there are strong relationships between the GRACE-DSI, the hydrometeorological factors, and burned area, and these relationships have certain time delay. It shows that climate change begins to play a role before a wildfire, and when this effect accumulates to a certain extent, it leads to a wildfire. We averaged the monthly GRACE-DSI, burned area, and hydrometeorological factors for a high burned year and its previous year (blue line in Figure 11). Before the abnormal high burned year, the GRACE-DSI and hydrometeorological factors all experience abnormal change. In Figure 11i, the burned area shows an upward trend that is significantly larger than the average level (green line), starting from July and continuing until October. ET and PPT show significant downward trends starting in February and March, respectively, and these trends continued until October. The GRACE-DSI shows a decreasing trend in May, which continues until December. With PPT and ET changing significantly, there are also significant changes in SM, VPD, RH, and PCW since March. Among them, the trends of SM, RH, PCW (decreasing), and VPD (increasing) are just opposite to each other.

Table 2. The maximum correlation coefficients and lag months between different variables in high burned years and the average level shown in Figure 11 (The correlation coefficient results have passed the 95% confidence level).

Variables	Correlation Coefficients		Lag Months
	High Burned Year	Average Level	High Burned Year	Average Level
PPT vs. GRACE-DSI	0.52	0.50	4	1
ET vs. GRACE-DSI	0.90	0.52	4	2
GRACE-DSI vs. Burned Area	−0.50	−0.32	4	2
PPT vs. SM	0.67	0.50	1	1
ET vs. SM	0.87	0.81	2	4
GRACE-DSI vs. SM	0.86	0.71	1	2
ET vs. VPD	−0.92	−0.51	0	4
ET vs. RH	0.82	0.79	0	4
SM vs. PCW	0.67	0.62	0	1
GRACE-DSI vs. VPD	−0.83	−0.46	4	2
GRACE-DSI vs. RH	0.65	0.49	4	1
GRACE-DSI vs. PCW	0.50	0.49	4	2
GRACE-DSI vs. MT	−0.61	−0.36	0	2
VPD vs. Burned Area	0.59	0.74	1	3
RH vs. Burned Area	−0.62	−0.51	1	3
PCW vs. Burned Area	−0.59	−0.58	1	3
MT vs. Burned Area	0.73	0.10	4	3

shows the maximum correlation coefficients and lag months between the GRACE-DSI, burned area, and hydrometeorological factors during the high burned years and the average level. The impact of climate change on wildfire in Australia shows a clear path. A reduction in PPT has a greater impact on SM (CC = 0.67), and the continuous accumulation of this impact leads to terrestrial water deficit (TWD), which leads to drought (PPT and GRACE-DSI, CC = 0.52; SM and GRACE-DSI, CC = 0.86). Drought may reduce ET (GRACE-DSI and ET, CC = 0.90), thereby reducing atmospheric water vapor (ET and VPD, CC = −0.92; ET and RH, CC = 0.82). Figure 11e,f show that there are clear decrease and increase in RH and VPD during the high burned years. Insufficient moisture in the air is conducive to the occurrence and spread of wildfire (burned area and RH, CC = −0.62; burned area and VPD, CC = 0.59) [17]. On the other hand, a SM deficit changes the PCW [69], causing the vegetation to become abnormally dry and, thus, preparing sufficient combustible materials for the occurrence and spread of wildfire (SM and PCW, CC = 0.67; PCW and burned area, CC = −0.59). A high MT speeds up the above process (MT and burned area, CC = 0.73).

When comparing the performance of the GRACE-DSI and the hydrometeorological factors during the high burned years and the average level (Figure 11 and Table 1), the burned area has a stronger correlation with the GRACE-DSI, RH, PCW, and MT (−0.50, −0.62, −0.59, and 0.73) than the ones (−0.32, −0.51, −0.58, and 0.10) in the average level. In particular, the correlations between the GRACE-DSI, MT, and burned area greatly improve during the high burned years. It indicates that the GRACE-DSI and MT are more sensitive to abnormal changes in the burned area. Additionally, the GRACE-DSI has a closer relationship with SM, RH, VPD, and PCW (0.86, −0.83, 0.65, and 0.50) during the high burned years. This means that the GRACE-DSI has a positive response to the dry and wet changes in the study region before wildfires. Therefore, it is feasible to use the GRACE-DSI to detect the impact of climate change on wildfires, and even this index can be used to achieve an early warning of wildfires in Australia.

5. Discussion

With the intensification of global warming in recent decades, extreme climates have occurred frequently around the world. Australia is located between the Pacific and Indian Oceans; therefore, its climate change is more susceptible to extreme climates associated with both oceans. When we discuss the impact of climate change on wildfires in Australia, it is important to consider the role of extreme climates (ENSO, PDO, and IOD). When combining Figure 12a–c, an El Niño event is prone to less PPT and higher MT, leading to drought, while a La Niña event is just the opposite. The burned area anomalies all appear positive during or after an El Niño event. It is due to the fact that an El Niño event creates a fire-prone climate. However, there is no significant connection between the IOD, PDO, droughts, and wildfires (Figure 12). From

Table 3. The maximum correlation coefficients and lag months between extreme climate and PPT, MT, burned area, GRACE-DSI (The correlation coefficient results have passed the 95% confidence level).

Variables	Correlation Coefficients	Lag Months
ENSO vs. PPT	−0.56	2
ENSO vs. MT	0.57	4
ENSO vs. GRACE-DSI	−0.33	3
ENSO vs. Burned area	0.13	5
DMI vs. PPT	−0.42	2
DMI vs. MT	0.47	6
DMI vs. GRACE-DSI	−0.23	5
DMI vs. Burned area	0.25	1
POD vs. PPT	−0.45	5
POD vs. MT	0.35	4
POD vs. GRACE-DSI	−0.31	2
POD vs. Burned area	0.18	5

, the ENSO has a stronger correlation with PPT and MT (−0.56 and 0.57) in Australia than the DMI and the PDO (−0.42 and 0.47; −0.45 and 0.35). It explains that the ENSO has a greater impact on Australia’s climate, and the IOD and POD have the same impact on PPT and MT as the ENSO, that is, a positive event leads to drought. The correlation coefficients results (ENSO and GRACE-DSI, −0.33; DMI and GRACE-DSI, −0.23; and PDO and GRACE-DSI, −0.31) verify the above conclusion (Table 3). The correlation coefficients between the burned area and the ENSO, the DMI, and the PDO (0.13, 0.25, and 0.18) show that the three extreme climates have a positive correlation with the burned area. Previous studies indicate that anomalous change in sea surface temperature is positively correlated with wildfires in adjacent continental regions [70,71]. An abnormal change in sea surface temperature is a basic element of extreme climate. Since an extreme climate event does not directly affect wildfire, and it indirectly affects wildfire by changing the local climate environment (PPT and burned area, CC = −0.61; MT and burned area, CC = 0.73), the correlation between an extreme climate event and the burned area is not strong.

Previous studies have shown the strong connection between drought caused by extreme climates and extreme fire phenomena on the local scale [17,72]. Chen et al. [34] indicates that regional TWSC is closely related to the number of wildfires in an active period, and the correlation between the two is a negative correlation. An extreme climate event leads to less PPT and higher MT, which causes SM deficit, and drought occurs (GRACE-DSI) [73]; then, the drought (GRACE-DSI) causes a decrease in ET, thereby reducing RH and increasing VPD (Figure 11c,e,f). A lower atmospheric humidity level would create a climate condition that is more favorable for wildfire occurrence and spread [17]. In addition, there is another way that drought (GRACE-DSI) can affect wildfire behavior, that is, drought (GRACE-DSI) leads to SM scarcity, and the water source that the vegetation can absorb is reduced (Figure 11c,d). Due to the impact of the reduction in PCW, the leaf area of vegetation is greatly reduced, and a large number of leaves fall off, which prepares sufficient combustible materials for wildfire [21]. Therefore, drought (GRACE-DSI) is a hydrometeorological factor that causes wildfire occurrence and spread. The GRACE-DSI is an important indicator to measure the occurrence and severity of droughts, and it has the characteristics of being simple and easy to understand [37]. When compared to the wildfire early warning system based on the traditional hydrometeorological factors (PPT, RH, MT, etc.), the GRACE-DSI can better reflect the overall dry and wet changes in the study region. Therefore, the GRACE-DSI has certain advantages in monitoring local climate changes before wildfires and establishing an early warning of wildfires.

6. Conclusions

In our study, we used six GRACE/GRACE-FO solutions, burned area, extreme climate indices, and seven hydrometeorological factors to investigate the connection between extreme climates and wildfires in Australia during the period between 2003 and 2020. The following conclusions can be drawn:

(1)

In terms of spatial distribution, Australia’s wildfires are mainly concentrated in the north, with sporadic wildfires in the southeast. In terms of temporal distribution, Australia’s wildfires are mainly concentrated in October and November. In 2011 and 2012, two of Australia’s worst wildfires occurred during the 18-year study period.

(2)

TWSC and the seven hydrometeorological factors are strongly correlated with burned area on a seasonal scale. Before the occurrence of wildfires, the regional climate generally changes abnormally, especially during a high burned year.

(3)

Droughts often lead to an increased chance of wildfires, which not only provides an external environment that is easy for wildfires to occur, but also provides an accumulation of combustibles for the occurrence and spread of wildfires.

(4)

An extreme climate event (ENSO, IOD, and PDO) is an important reason for the abnormal changes in regional climate, which has a strong influence on PPT and MT in Australia. An extreme climate event can lead to less PPT and higher MT, causing severe droughts.

(5)

The GRACE-DSI is a scientifically valid, easy-to-understand indicator of the occurrence and severity of droughts. Therefore, it can be used to evaluate the risk of wildfire occurrence.

However, limited by the mission period of the GRACE/GRACE-FO satellites, our study only used 18-year data for research. With the accumulation of follow-up data, we can obtain a longer-time observation series to explore the law of wildfire evolution in future research. Our results provide a new idea for establishing an early warning of regional wildfires and also have a great scientific significance for the impact of regional climate change on local ecological environment.