The Influence of Climate Change on Forest Fires in Yunnan Province, Southwest China Detected by GRACE Satellites

Abstract

Yunnan province in China has rich forest resources but high forest fire frequency. Therefore, a better understanding of the relationship between climate change and forest fires in this region is important for forest fire prevention. This study used the Gravity Recovery and Climate Experiment (GRACE) terrestrial water storage change (TWSC) data to analyze the influence of climate change on forest fires in the region during 2003–2016. To improve the accuracy and reliability of GRACE TWSC data, we used the generalized three-cornered hat (GTCH) and the least square method to fuse TWSC data from six GRACE solutions. The spatiotemporal variation of forest fires during 2003–2016 was investigated using burned area data. Then, the relationship between burned area and hydrological and climatic factors was analyzed. The results indicate that more than 90% of burned areas are located in northwestern and southern Yunnan (NW and S). On the seasonal scale, forest fires are mainly concentrated in January–April (dry season) and the burned area is negatively correlated with precipitation (correlation coefficient r = −0.83 (NW) and −0.51 (S)), relative humidity (r = −0.79 (NW) and −0.92 (S)), GRACE TWSC (r = −0.57 (NW) and −0.73 (S)) and evapotranspiration (r = −0.90 (NW) and −0.35 (S)). However, the burned area has no significant correlations with the above four factors on the interannual scale. The composite analysis suggests that the extreme climate affects precipitation, evapotranspiration and TWSC in this region, thereby changing water storage of the air in this region, leading to the formation of an environment prone to forest fires. Such conditions have led to an increase in the burned area in the above region. We also found that the difference between TWSC in high- and low-fire years is much greater than the precipitation in the same period. The above results show that GRACE satellites can detect the influence of climate change on forest fires in Yunnan province.

1. Introduction

Forest fires are serious natural disasters, which have caused enormous damage to human society, ecological environments and economic development [1]. According to statistics, there were more than 139,000 forest fires in North America each year, with a total area of 4.2 million hm² [2]. In 2019, the Amazon rainforest fire continued to burn for more than 21 days, and the burned area exceeded 1 million hm². It has caused serious damage to the local ecological environment and may be a greater influence on the local and global climate [3]. During 2019–2020, the wildfire in Australia killed at least 33 people and about 1 billion animals, and the affected area was more than 10 million hm². The heavy smoke caused by the wildfire lasted for several weeks and blocked many cities including Sydney, and the local ecological environment was severely damaged [4]. Therefore, it is of great significance to study the driving factors of forest fires for effectively controlling the occurrence of forest fires.

Forest fire occurrence is driven by many factors including weather condition [5], human activities [6], fuel characteristics [7], fire management [8], land use changes [9] and climate changes [10,11]. Among these, climate change is a key factor contributing to forest fires [12,13], which mainly affects forest fires in two ways. On the one hand, climate change affects the occurrence of forest fires by changing the dryness of combustibles through temperature, precipitation (PPT), evapotranspiration (ET), etc. [14,15,16]. On the other hand, climate change can change the composition, properties and biomass of combustibles by affecting the structure and growth of vegetation, which affects the behavior of forest fires [17,18,19]. Swetnam et al. [20] and Jones et al. [21] indicate that the large-scale climate model of Pacific Decadal-scale oscillations and strong El Niño/La Niña events is closely related to the frequency and severity of forest fires. The severity of forest fires has a certain correlation with abnormal ocean surface temperature and Palmer Drought Severity Index [5,22,23]. Gillett et al. [24] indicate that the increase in the area of burned areas in Canada from 1920 to 1990 is related to the trend of regional warming over the same period. Heyerdahl et al. [25] explain that the frequency of forest fires in inland forests is synchronized with the changes in the dry and wet climates in spring and summer. Westering et al. [26] analyzed the influence of climate change in the western United States on forest fires. The results show that the increase in average temperature and the extension of dry periods have increased the number of forest fires in the western United States by nearly three times, and the burned area has increased by nearly six times. Pausas et al. [27] studied the relationship between climate change and forest fires in the Iberian Peninsula from 1950 to 2000. The results show that the summer temperature and PPT are significantly related to the frequency of forest fires and the fire area in the same period.

In 2002, the Gravity Recovery and Climate Experiment (GRACE) was implemented to detect the changes of Earth’s gravity field on medium and long scales [28]. These changes mainly come from terrestrial water, so it is generally considered that Earth’s surface mass changes detected by GRACE are equivalent to terrestrial water storage change (TWSC) [29]. As an important part of the regional water cycle, TWSC can not only reflect the water balance relationship between water inputs (such as, PPT, snowfall, etc.) and water outputs (such as runoff, ET, etc.) but also was used to evaluate the influence of climate change on regional hydrological changes [30,31]. Therefore, it is feasible to use GRACE TWSC to study the influence of climate change on forest fires.

Yunnan province is one of the major forest regions and the most frequent and hardest-hit region of forest fires in China [32]. There is high forest coverage rate and rich forest resources. Forest fires in this region are characterized by a long fire-prevention period, complex terrain and difficulty in fighting the fires. During 1996–2008, there were 135.8 thousand forest fires that occurred, 0.1358 million hm² of forest burned, 67 people died from fighting fires and 142 people were burned [33]. Yunnan is located on a plateau with complex topography; there are significant regional differences and vertical changes in climate, and regional climate changes are inconsistent [32]. Therefore, studying the relationship between forest fires and climate change will help to scientifically understand the laws of temporal and spatial changes of forest fires in this region, and provide scientific data support for government departments’ forest fire prevention and macro decision making.

However, there are few studies about using GRACE TWSC to detect the influence of climate change on forest fires [34,35]. In particular, research about Yunnan province has not been reported in the literature to the best of our knowledge. Due to the differences between different GRACE solutions, it may lead to unreliable results [36]. Therefore, we first used the generalized three-cornered hat method (GTCH) to evaluate the uncertainly of six GRACE solutions in this paper, then the least squares method was used to integrate the TWSC results from the above GRACE solutions. Moreover, the fusion result was used to investigate the climatic characteristics of Yunnan province. Finally, the connection between climate change and forest fires was analyzed. This paper is organized as follows. In Section 2, Section 3 and Section 4, the study area, data and study methods are introduced separately. Section 5 presents the spatial and temporal distribution of fire burned area, the climate change characteristics in Yunnan and the analysis results of the influence of climate change on forest fires. The discussion and conclusion are provided in Section 6 and Section 7, respectively.

2. Study Area

Yunnan is located in Southwest China, approximately at 21°–30°N and 97°–106°E, with a total area of about 390,000 km² (see Figure 1). This region is a mountainous plateau region: most of the region is mountainous, and the overall terrain tends to be higher in the west and lower in the east [32]. Warm in winter and cool in summer, all seasons are like spring. There are also complex vegetation types and rich forest resources in this region. Because it is located at the intersection of frigid and tropical zone plant regions, plants of the frigid, warm and hot zones have grown. In addition, the distribution of vegetation in this region has obvious horizontal zonality. From north to south, the vegetation types are cold-temperate coniferous forest, subtropical semi-humid evergreen broad-leaved forest, subtropical monsoon evergreen broad-leaved forest, tropical monsoon forest and tropical rainforest [37].

3. Data

3.1. GRACE Data

To obtain monthly 1° × 1° gridded TWSC data, we used four GRACE RL06 spherical harmonic (SH) solutions provided by the Center for Space Research at the University of Texas at Austin (CSR), Helmholtz–Centre Potsdam–German Research Centre for Geosciences (GFZ), Jet Propulsion Laboratory (JPL) and Institute of Geodesy at Graz University of Technology (ITSG), respectively. Before calculating the gridded TWSC data, we needed to perform several types of preprocessing on the four SH solutions. Firstly, C₂₀ of SH coefficients have been replaced with those of satellite laser ranging [38]. The degree-1 coefficients have been replaced using the results from Swenson et al. [39]. To minimize the effects of high-frequency and correlated errors, the combined filtering (300 km Fan filter and de-striping method P3M6) has been applied for smoothing the SH solutions [40]. To reduce the signal attenuation effect caused by order truncation and filter smoothing, the GRACE TWSC data was corrected by using the scale factor approach [41].

In this study, we also obtained the monthly 1° × 1° global gridded TWSC data derived from two Mascon solutions provided by CSR and JPL. The JPL Mascon is Release 06 Version 1.0 with the coastal resolution improvement filter [42], and the CSR Mascon is Release 06 Version 2.0 [43]. Unlike SH solutions, Mascon solutions do not require additional processing and can be directly read and used. A number of measures have been taken to improve the accuracy of Mascon solutions. Firstly, the C₂₀ coefficients have been replaced with the solutions from SLR. Secondly, the degree-1 coefficients have been estimated using the results of TN-13a [44] and the approach of Sun et al. [45] and Swenson et al. [39]. Thirdly, the ice rebound correction has adopted ICE6G-D model [46].

Therefore, the four SH solutions and two Mascon solutions were used to obtain monthly GRACE TWSC gridded data in Yunnan from January 2003 to December 2016 in this study. For convenience, the four SH solutions and two Mascon solutions are termed CSR-SH, GFZ-SH, JPL-SH, ITSG-SH, CSR-M and JPL-M.

3.2. Burned Area Data

The monthly global gridded burned area data from January 2003 to December 2016 was provided by the Global Fire Emission Database Version 4.1s (GFED4.1s) [47,48] in our study, which was used to estimate the fire-affected area during forest fires in Yunnan. The GFED4.1s product is a collection of multiple satellite detection data, including the Moderate-resolution Imaging Spectroradiometer (MODIS), the Tropical Rainfall Measuring Mission (TRMM) Visible and Infrared Scanner (VIRS) and the Along-Track Scanning Radiometer (ATSR) family of sensors [49]. This database with spatial resolution 0.25° × 0.25° provides four components: burned area, carbon emissions, biosphere fluxes and other ancillary data sets.

3.3. In Situ PPT Data

The monthly gridded PPT data with spatial resolution 0.5° × 0.5° from January 2003 to December 2016 was provided by the China National Meteorological Science Data Center in our study. The datasets come from the monthly PPT data at national-level stations nationwide from 1961 to the present, which were collected and compiled by the China National Meteorological Information Center.

3.4. Evapotranspiration (ET) Data

We used the ET gridded data with spatial resolution 0.25° × 0.25° from January 2003 to December 2016 derived from the Global Land Evaporation Amsterdam Model (GLEAM) 3.5a [50,51] in our study. The model is a set of algorithms that separately estimate the different components of land evaporation: transpiration, bared-soil evaporation, interception loss, open-water evaporation and sublimation. Additionally, GLEAM provides surface and root-zone soil moisture, potential evaporation and evaporative stress conditions.

3.5. Relative Humidity Data

To further investigate the relationship between humidity changes in air and forest fires, we used the relative humidity gridded data derived from the National Centers for Environmental Prediction-Department of Energy (NCEP-DOE) Reanalysis 2 project. The data are provided by Physical Sciences Laboratory at the National Oceanic and Atmospheric Administration (NOAA). The NCEP-DOE Reanalysis 2 project is using a state-of-the-art analysis/forecast system to perform data assimilation using past data from 1979 through the previous year [52], and the temporal coverage of the data are 4-times daily, daily and monthly. In our study, we used the monthly relative humidity gridded data with spatial resolution 2.5° × 2.5° from January 2003 to December 2016.

3.6. Atmospheric Column Water Vapor

To investigate the role of water storage in regulating air and fuel moisture, the global atmospheric column water vapor monthly gridded observations with spatial solutions 1° × 1° were used in our study, which are derived from the Atmospheric Infrared Sounder (AIRS) [53]. In our study, the monthly atmospheric column water vapor gridded data were used from January 2003 to December 2016.

3.7. Extreme Climate Index

Since PPT and ET of Yunnan are affected by the El Niño-Southern Oscillation (ENSO) and the Indian Ocean Dipole [40,54], we included the above two extreme climates into the analysis scope when studying the influence of climate change on forest fires in Yunnan. The monthly Niño 3.4 index data indicate the magnitude of ENSO, which is provided by the National Oceanic and Atmospheric Administration (NOAA). The Indian Ocean Dipole Model Index (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 monthly DMI data were also provided by NOAA. The above index data in our study are from January 2003 to December 2016.

4. Method

4.1. Fusing Different Datasets

Due to discrepancies in different GRACE TWSC datasets, inconsistent results may appear, which may confuse the analysis results. Hereafter, to minimize the discrepancies caused by different datasets and improve the reliability of the results, we used the GTCH method to estimate the relative uncertainties of the six GRACE TWSC datasets, and fused them according to their uncertainties using the least square method.

4.1.1. The GTCH Method

This method can estimate relative covariances of different datasets without a priori knowledge as long as there are at least three datasets [56,57]. Suppose there are several different observation series, denoted by $x_i$, $i$ correspond to different data. Each observation series can be expressed as:

$$x_i=\widehat{x}+{\mathit{\epsilon }}_i,i=1,2,\cdots ,N$$

(1)

where $\widehat{x}$ is the real signal, ${\mathit{\epsilon }}_i$ is the error of the $i$ observation series (0 means white noise) and N is the number of observation series. Since the true signal cannot be obtained, any observation series was selected as the reference signal, and the choice of the reference series does not influence the final results [58]. Then, the difference series between the remaining observation series and the reference series was calculated [59].

$$y_i=x_i-x_R={\mathit{\epsilon }}_i-{\mathit{\epsilon }}_R,i=1,2,\cdots ,N-1$$

(2)

where $x_R$ is the reference observation series and ${\mathit{\epsilon }}_R$ is the errors of reference observation series. The TWSC series derived from the JPL Mascon solution was selected as the reference series. N − 1 difference series were stored in the following matrix.

$$Y=\left[\begin{array}{cccc}{y}_{11}& y_{12}& \cdots & y_{1\left(N-1\right)}\\ y_{21}& y_{22}& \cdots & y_{2\left(N-1\right)}\\ \cdots & \cdots & \ddots & \cdots \\ y_{M1}& y_{M2}& \cdots & y_{M\left(N-1\right)}\end{array}\right]$$

(3)

where M is the number of observations in the series. The covariance matrix of difference series is:

$$S=\mathrm{cov}\left(y_{iN},y_{jN}\right)=\left[\begin{array}{cccc}{s}_{11}& s_{12}& \cdots & s_{1\left(N-1\right)}\\ s_{21}& s_{22}& \cdots & s_{2\left(N-1\right)}\\ \cdots & \cdots & \ddots & \cdots \\ s_{\left(N-1\right)1}& s_{\left(N-1\right)2}& \cdots & s_{\left(N-1\right)\left(N-1\right)}\end{array}\right]\left(i,j=1,2,\cdots ,N-1\right)$$

(4)

where $\mathrm{cov}\left(\circ \right)$ is the covariance operator. $s_{i,j}$ is the variance estimate $\left(i\ne j\right)$ or covariance estimate $\left(i=j\right)$, that is, the variance or covariance estimate of difference series between the TWSC series derived from different GRACE solutions and the reference series. The unknown noise covariance matrix R was introduced, and its relationship with S is [60]:

$$S=J•R•J^T$$

(5)

where the matrix J is:

$$J_{N-1,N}=\left[\begin{array}{lllll}1\hfill & 0\hfill & \cdots \hfill & 0\hfill & -1\hfill \\ 0\hfill & 1\hfill & \cdots \hfill & 0\hfill & -1\hfill \\ ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill \\ 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & -1\hfill \end{array}\right]$$

(6)

The matrix R is:

$$R_{N,N}=\left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1N}\\ r_{12}& r_{22}& \cdots & r_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ r_{1N}& r_{2N}& \cdots & r_{NN}\end{array}\right]$$

(7)

From Equation (5), the following relationship can be obtained:

$$r_{ij}=s_{ij}-r_{NN}+r_{iN}+r_{jN}$$

(8)

Both R and S are real symmetric matrices by definition. There are $N\times \left(N+1\right)/2$ unknown parameters for R, but there are only $N\times \left(N-1\right)/2$ equations for S. Therefore, Equation (5) is underdetermined. As a results, the remaining $N$ free parameters need a reasonable way to obtain a unique solution [60].

To guarantee the positive definiteness of R, Galindo and Palacio [61] proposed an important constraint on the solution space for free parameters based on the Kuhn–Tucker theorem.

$$H_2\left(r_{1N},\cdots r_{NN}\right)=\frac{H_1\left(r_{1N},\cdots r_{NN}\right)}{K}>0$$

(9)

where $K=\sqrt[N-1]{\left|S\right|}$, which is introduced for a better numerical solution. $H_1\left(r_{1N},\cdots r_{NN}\right)$ is given by [62]:

$$H_1\left(r_{1N},\cdots r_{NN}\right)=\frac{\left|R\right|}{\left|S\right|}=r_{NN}-\left[r_{1N}-r_{NN},\cdots ,r_{\left(N-1\right)N}-r_{NN}\right]S^{-1}{\left[r_{1N}-r_{NN},\cdots ,r_{\left(N-1\right)N}-r_{NN}\right]}^T$$

(10)

This condition constrains the free parameters in the solution domain, but it is still not enough to obtain the unique solution of the free parameters [59]. It is necessary to give the optimal selection criteria to obtain the unique parameter solution. Tavella and Premoli [58] proposed that the determination of the free parameters is based on the minimum “global correlation” of all observation series and the positive definiteness of R. Therefore, the following objective function was defined and minimized to determine the free parameters [61].

$$F_1\left(r_{1N},\cdots r_{NN}\right)=\frac{1}{K^2}{\displaystyle \sum _{i<j}^N{r}_{ij}^2}$$

(11)

In order to make the initial value within the constraints, the initial value of the iterative calculation was set to [58]

$$\left\{\begin{array}{l}{r}_{iN}^0=0,i<N\hfill \\ r_{NN}^0=\frac{1}{2s^{\ast }},s^{\ast }=\left[1,\cdots ,1\right]S^{-1}{\left[1,\cdots ,1\right]}^T\hfill \end{array}\right.$$

(12)

Minimize the objective function (Equation (11)) to obtain a set of free parameter solutions under constraints (Equation (9)), that is, the variance of the uncertainty of different observation series. Other unknown elements in R can be obtained by Equation (8).

4.1.2. Data Fusion

We applied the GTCH method to estimate the relative covariance of TWSC derived from six GRACE datasets. Then we fused the six TWSC results by taking a weighted average of them [44].

$$\mathrm{TWSC}={\displaystyle \sum _{i=1}^N{p}_i•TWS{C}_i}$$

(13)

where $TWS{C}_i$ and $p_i$ are the TWSC results derived from the individual GRACE solution and its corresponding weight, respectively. The weights were determined based on estimated variances.

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

(14)

where $r_{ii}\left(i=1,2,\cdots ,6\right)$ is the variance of the $i$th TWSC time series estimated by the GTCH method. The above process was performed grid by grid until we fused the six datasets on all the grid nodes.

4.2. Composite Analysis

In the climate change research, composite analysis was used to explore the salient characteristics of certain special years [63]. In our study, we first standardized the time series of burned area. Then, the special years were selected according to the rule of greater than 0.5 times or less than 0.5 times the standard deviation based on the standardized time series of burned area, which were marked as high- and low-burned area years, respectively. Finally, the average of meteorological and hydrological elements of these special years were compared with the ones of all years to analyze the difference of these elements during high- and low-burned area years.

4.3. Pearson Correlation Analysis

In our study, we used the Pearson correlation analysis to evaluate the relationship between two sets of data. In statistics, the covariance of two sets of data is generally divided by the standard deviation of the corresponding data series to obtain the Pearson correlation coefficient [4]. For $x_i\left(i=1,2,\cdots ,n\right)$ and $y_j\left(j=1,2,\cdots ,n\right)$ two sets of time series, the Pearson correlation coefficient $r$ can be estimated by the following expression:

$$r=\frac{{\displaystyle \sum _{i=1,j=1}^n\left(x_i-\overline{x}\right)\left(y_j-\overline{y}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}}\sqrt{{\displaystyle \sum _{j=1}^n{\left(y_j-\overline{y}\right)}^2}}}$$

(15)

where $n$ is the number of observations in the time series $x_i$ and $y_j$. $\overline{x}$ and $\overline{y}$ are the average values of the time series $x_i$ and $y_j$ during the study period. The value range of $r$ is between −1 and 1.

4.4. Time Series Analysis

The time series of observation signals contains the long-term change trend, acceleration, seasonal change (annual change and semi-annual change) and residual signal. The corresponding signals can be extracted from the linear fitting model. The decomposition is expressed as follows [64,65]:

$$TWSC\left(t\right)=a_0+a_1\left(t-t_0\right)+a_2{\left(t-t_0\right)}^2+a_3\cos \left(2\pi t\right)+a_4\sin \left(2\pi t\right)+a_5\cos \left(4\pi t\right)+a_6\sin \left(4\pi t\right)+\epsilon$$

(16)

where $TWSC\left(t\right)$ is the time series of TWSC; $t$ is the time; $t_0$ is the midpoint of the entire research period; $\epsilon$ is the error and other signal; and $a_0$, $a_1$, $a_2$, $a_3$, $a_4$, $a_5$ and $a_6$ are the unknown parameters. $a_0$ is constant term, $a_1$ is long-term trend change, $a_2$ is the acceleration, $a_3$ and $a_4$ are annual terms and $a_5$ and $a_6$ are semi-annual terms.

5. Results

5.1. GRACE Solutions Fusion

Figure 2 shows the spatial distribution of uncertainties of six GRACE-based TWSC results estimated by the GTCH method in Yunnan. It can be seen that the uncertainties of TWSC results derived from the four SH solutions are smaller than those of two Mascon solutions. Specifically, the uncertainties of TWSC results from the four GRACE SH solutions are smaller than 4 cm, while TWSC results from the two GRACE Mascon solutions show uncertainties greater than 4 cm in most regions. Except for GFZ-SH, the regions with greater uncertainty are mainly concentrated in western Yunnan. Specifically, the uncertainty of TWSC results from two GRACE Mascon solutions are higher than 6 cm in the above region.

To evaluate the overall uncertainty of six GRACE solutions, we sorted the grid values of certainty in the study region in ascending order and took the median one to represent the overall uncertainty in this region. The results are presented in

Table 1. Uncertainties of TWSC results derived from six GRACE solution and fusion result estimated by the GTCH method.

GRACE Solution	CSR-SH	GFZ-SH	JPL-SH	ITSG-SH	CSR-M	JPL-M	Fusion Result
Median(cm)	2.66	2.63	3.40	3.02	5.11	4.62	0.85

. The six GRACE solutions were sorted in ascending order of the uncertainty of TWSC results, and their arrangement is GFZ-SH (2.63 cm), CSR-SH (2.66 cm), ITSG-SH (3.02 cm), JPL-SH (3.40 cm), JPL-M (4.62 cm) and CSR-M (5.11 cm). It can be seen that the results of GFZ-SH and CSR-SH are very close. In our study, we estimated only the relative accuracy. If you want absolute accuracy, you need ground truth data.

Based on the variances of six GRACE-based TWSC results estimated by the GTCH method, we fused these six results to yield the fused results. From Table 1, we found that the uncertainty of the fusion result is much smaller than the ones of six single solutions. Figure 3 shows the spatial distribution of uncertainties of the fusion result. Comparing Figure 2 and Figure 4, we found that all the uncertainties of fusion result are lower than 1.5 cm, while the ones of six single solutions are all higher than 2 cm. This suggests that the accuracy of the fusion result has been significantly improved.

Figure 4 shows the time series of TWSC results derived from six GRACE solutions and the fusion result in the study region. It can be seen that TWSC time series from six GRACE solutions and the fusion result have the similar magnitude. However, the magnitudes of TWSC time series from Mascon solutions are smaller than the ones from SH solutions. The fusion result is close to the TWSC results from SH solutions, because the uncertainties of the TWSC results from SH solutions are smaller than the ones of Mascon solutions, so the weights of SH solutions are larger than the ones of Mascon solutions in the fusion process. To further evaluate the fusion result, we calculated the correlation coefficients between the fusion result and the TWSC results from six single solutions (see

Table 2. The correlation coefficients between fusion results and TWSC results from six GRACE solutions.

GRACE Solution	CSR-SH	GFZ-SH	JPL-SH	ITSG-SH	CSR-M	JPL-M
Correlation coefficient for fusion results	0.9968	0.9949	0.9943	0.9889	0.9754	0.9741

). From Table 2, among six solutions, the correlation coefficient between the fusion result and CSR-SH is the highest (0.9968), followed by GFZ-SH, JPL-SH, ITSG-SH and CSR-M. The smallest one is JPL-M (0.9741). Therefore, the similar magnitude and high correlation suggest that the fusion result has good consistency with TWSC results from six GRACE solutions.

5.2. Spatial and Temporal Distribution of the Burned Area

Figure 5 shows the spatial distribution of monthly average burned area from 2003 to 2016 across Yunnan. From the figure, we found that the forest fires were mainly concentrated in northwest (NW), southern (S) and central Yunnan. Among the three regions, the region with the longest forest fire duration and the largest accumulation of burned area is region NW, followed by region S, and finally the central region. In region NW (25°–30°N and 97°–101°E), the high-fire months are February, March, October, November and December and the burned areas in most regions are larger than 18 ha, while the low-fire months are August and September and the burned areas in this region are smaller than 4 ha. In region S (21°–24°N and 98°–104°E), the high-fire months are February, March, April and December and the burned areas in most regions are larger than 4 ha, while the low-fire months are May, June, Auguest and October and the burned areas in most regions are smaller than 4 ha. There were no fires in June or October. In the central region, the forest fires only appeared in March, April and May. The largest burned area occurred in April. Except for the above three regions, there are fewer forest fires in other regions, and forest fires are mainly concentrated in region NW and S. The burned area in region NW is larger than the one in region S. Therefore, we focused on the above two regions in the following research.

Figure 6 shows the monthly average burned area during 2003–2016 in Yunnan, regions NW and S, and the percentages of burned area of regions NW and S divided by the ones in Yunnan. The maximum monthly average burned area (85,446 ha) and the minimum one (4110 ha) in Yunnan appeared in March and September, respectively; forest fires are mainly concentrated in winter and spring (January, February, March and April). The monthly burned areas in the high-fire months are more than 50,000 ha, which are much higher than those of other months. The maximum and minimum monthly burned areas in region NW appear in January and September, respectively; forest fires are mainly concentrated in autumn and winter (November, December, January and February). The monthly burned areas in these high-fire months are higher than 20,000 ha, while the ones in other months are lower than 14,000 ha. The largest proportion of burned area (82.06%) occurs in November, and the smallest one (13.15%) occurs in April. The average is 48.03%. In region S, the maximum monthly burned area appears in March, while the minimum one appears in September and October; forest fires are mainly concentrated in winter and spring (February, March and April). The monthly burned areas in these high-fire months are higher than 11,000 ha, while the ones in other months are lower than 5000 ha. The largest proportion of burned areas in Yunnan (26.25%) occurred in March, and the smallest one (3.16%) occurred in November. The average value is 12.62%.

According to the above numerical results, the low-fire months in the three regions are June, July, August and September. However, there is a significant difference between region NW and other two regions in the high fire months. Except for region NW, the temporal distribution of burned areas in Yunnan and region S are similar on the seasonal scale. This may be because region NW belongs to part of the Qinghai–Tibet Plateau. Therefore, its climatic characteristics are quite different from other parts of Yunnan [33]. According to the results of burned area percentage, it can be seen that regions NW and S are the regions that are more prone to forest fires in Yunnan, which is consistent with the results of Figure 5.

Figure 7 shows the temporal evolution of annual burned area for the period 2003–2016 in Yunnan, regions NW and S. According to composite analysis, the high and low burned years in three regions were determined. In Yunnan, the high-fire years are 2004, 2005, 2006, 2007 and 2009, while the low ones are 2008, 2012, 2014 and 2016. In region NW, 2003, 2006, 2007, 2009 and 2010 are the high-fire years, and 2004, 2008, 2011, 2012, 2013 and 2016 are the low fire years. In region S, the high-fire years are 2004, 2006, 2009, 2011 and 2015, while the low ones are 2007, 2008, 2010, 2012, 2014 and 2016.

Comparing Figure 7a–c the common high-fire years of three regions are 2006 and 2009, while the common low-fire years are 2012 and 2016. It can be seen that the interannual variation of forest fires in Yunnan is different from those in regions NW and S. Comparing Figure 7b,c, the common high-fire years of region NW and S are 2006 and 2009, while the common low-fire years are 2008, 2012 and 2016. 2004, 2011 and 2015 are the low-fire years in region NW, but they are the high-fire years in region S. 2010 is the high-fire year in region NW, but it is the low-fire year in region S. This indicates that the interannual variation of forest fires in region NW is different from that in regions S. This may be related to differences in climate and vegetation types in these regions. There are semi-arid and semi-humid climates in region NW, and the vegetation types in the region are mainly subtropical evergreen broadleaf and coniferous forest and boreal coniferous forests, while region S has a humid climate and its vegetation type is tropical seasonal rainforest [32,33].

5.3. The Relationship between Meteorological Factors and Burned Area

5.3.1. Long-Term Trend Change

To analyze the influence of climate change on forest fires in Yunnan, we estimated the spatial distribution of the long-term trend changes of GRACE TWSC, PPT, ET and burned area (Figure 8). We found that there is no significant trend in TWSC in most of Yunnan, and only regions NW and S show significant trend changes. Among them, TWSC in region NW shows a decreasing trend, while the one in region S shows an increasing trend. The long-term trend changes of PPT are consistent with the ones of GRACE TWSC in most regions, but there are significant differences in region S and western Yunnan. In regions NW and S only, there are significant growth changes in ET. Combining Figure 8a–c, it can be seen that due to the influence of climate change in region NW, no change in PPT and ET increased, so TWSC shows a decreasing trend, which creates a natural environment prone to forest fires. In addition, the increase in ET may reduce the water content of the vegetation itself, making the vegetation easier to burn. In region S, there is also a decrease in PPT and an increase in ET. However, TWSC shows an increasing trend. This may be related to the main stream of the Mekong River flowing through this region. The annual runoff of this river ranks ninth in the world, and runoff is an important part of TWSC [40]. From Figure 8d, we found that, similar to TWSC, PPT and ET, only regions NW and S have significant change trends in burned areas. The burned area shows a decreasing trend in most of these two regions. This is closely related to the active forest fire prevention measures adopted by the local government [66,67]. However, the burned area shows an increasing trend in some local regions, such as parts of western, northwest and southern Yunnan.

5.3.2. Seasonal Variations

Figure 9 shows the seasonal variations of TWSC, PPT, ET, relative humidity and burned area in region NW (a–d) and S (e–h) from 2003 to 2016. In the figure, PPT, ET, relative humidity and burned area data are the related change values (with the average value deducted). In Figure 9a, the monthly average TWSC in region NW is 1.61 cm, the maximum one (470.13 cm) appears in August and the minimum one (−359.99 cm) occurs in March. The monthly average burned area in region NW is approximately 0 ha, the maximum one (20,463.01 ha) occurs in January and the minimum one (−10,836.32 ha) appears in September. We found that the burned area change is positive in the months with the negative TWSC (from July to October), while the it is negative in the months with the positive TWSC (from December to March). The correlation coefficient between burned area and PPT on the seasonal scale was estimated and its value is −0.57 (see

Table 3. The correlation coefficients between burned area (BA), PPT, ET, relative humidity (RH) and TWSC on the seasonal scale in region NW.

Correlation Coefficient	BA	TWSC	PPT	ET	RH
BA	1	−0.57	0.83	−0.90	−0.79
TWSC	−0.57	1	0.67	0.36	0.89
PPT	−0.83	0.67	1	0.89	0.81
ET	−0.90	0.36	0.88	1	0.62
RH	−0.79	0.89	0.81	0.62	1

). This explains that burned area and TWSC have a strong negative correlation on the seasonal scale.

The season with the most PPT is summer, followed by autumn and spring, and the least is winter (Figure 9b). The maximum and minimum PPT are 291.72 cm (July) and −161.51 cm (December), respectively. When the monthly PPT is negative, the monthly burned area is positive, vice versa. This suggests that there is a strong negative correlation between PPT and burned area (−0.83, see Table 3). In Figure 9c,d, the maximum ET and humidity are 262.76 cm (July) and 11.37% (August), respectively, while the minimum ones are 270.42 cm (December) and −13.72% (March), respectively. It can be seen that there are positive burned areas in the negative ET and relative humidity months (from November to March), while there are negative ones in the positive ET and relative humidity months (from April to September, ET; from June to October, relative humidity). The correlation coefficient results also support the above conclusion (burned area and ET, −0.90; burned area and humidity relatively, −0.79; Table 3).

We also compared the monthly burned area with monthly TWSC, PPT, ET and relative humidity in region S (Figure 9e–h). In this region, the maximum burned area (16,750.48 ha) appears in September, while the minimum one (−5202.83 ha) occurs in March. We can see that there are similar relationships between burned area and TWSC, PPT, ET and relative humidity. In other words, like in region NW, the positive burned areas generally occur in the months with negative TWSC, PPT, ET and relative humidity, while the negative ones appear in the months with positive TWSC, PPT, ET and relative humidity. From

Table 4. The correlation coefficients between burned area (BA), PPT, ET, relative humidity (RH) and TWSC on the seasonal scale in region S.

Correlation Coefficient	BA	TWSC	PPT	ET	RH
BA	1	−0.73	−0.51	0.35	−0.92
TWSC	−0.73	1	0.37	0.24	0.82
PPT	−0.51	0.37	1	0.96	0.67
ET	−0.35	0.24	0.96	1	0.52
RH	−0.92	0.82	0.67	0.52	1

, the correlation coefficients between burned area and TWSC, PPT, ET and relative humidity are −0.73, −0.51, −0.35 and −0.92, respectively. However, we found that the correlation between burned area and ET is weaker than the one in region NW. Relative to region NW, the correlations between burned area and TWSC and relative humidity are stronger, while the one between burned area and PPT are weaker. The three above correlation coefficients are smaller than −0.5.

According to Figure 9, Table 3 and Table 4, we analyzed the relationship between forest fires and climate change in regions NW and S. The climate of Yunnan belongs to the Subtropical Plateau Monsoon Climate, and the dry and rainy seasons in Yunnan are relatively distinct. According to Figure 9b,f, the rainy season is from May to September and the dry season is from October to April in the two regions. From May, PPT begins to become positive and increase gradually. During the same period, ET also begins to increase as the temperature rises. The increase in PPT and ET lead to an increase in water vapor in the air, so relative humidity in the two regions rises in June (see Figure 9d,h). Similarly, TWSC also starts to increase in this month (see Figure 9a,e). The environment at this time is not conducive to the occurrence of fire. Therefore, the burned area is at a negative value during the rainy season.

The turning point appears in October. From October, PPT and ET begin to become negative. The relative humidity also becomes negative from November (NW) and January (S), and TWSC begins to become negative in December (NW) and January (S). The burned area begins to become positive in November (NW) and February (S). We found that there are significant relationships between climate change and forest fires in the two regions. Firstly, PPT and ET decrease due to climate change (from rainy season to dry season); 1–3 months later, the relative humidity and TWSC begin to fall and become negative. Immediately after, the burned area begins to increase. The correlation coefficients results (Table 3 and Table 4) show that the increase in PPT and ET causes an increase in relative humidity, and then the increase in relative humidity is reflected in TWSC. TWSC and relative humidity have a stronger negative correlation with burned area.

5.3.3. Interannual Variance

Figure 10 shows the interannual variations of TWSC, PPT, ET, relative humidity and burned area in region NW (a–d) and S (e–h) from 2003 to 2016. From these figures, the burned area has no significant correlation with TWSC, PPT, ET and relative humidity. In Figure 10a, the years with the opposite change direction in PPT and burned area are only 2003, 2006, 2008, 2009, 2014 and 2016, and there is no significant rule from the change range of the two. The same situation also appears between burned area and ET (see Figure 10b,c). From Figure 10d, it can be seen that, except for 2004 and 2010, the change direction of burned area and relative humidity are opposite. This shows that there is a significant negative correlation (−0.64,

Table 5. The correlation coefficients between burned area (BA), PPT, ET, relative humidity (RH) and TWSC on the interannual scale in region NW and S.

Correlation Coefficient	Region	TWSC	PPT	ET	RH
BA	NW	−0.15	−0.08	−0.27	−0.64
	S	−0.44	0.03	0.02	−0.11

) between the two in the interannual scale. From Figure 10e, we found that the years with the opposite change direction in PPT and burned area in region S are greater than the ones in region NW, which are 2004, 2005, 2006, 2008, 2009, 2012, 2013, 2014 and 2016. This suggests that the negative correlation between TWSC and burned area in region S is stronger than that in region NW. The results of Table 5 also verify this point (−0.15 in region NW and −0.44 in region S). In Figure 10f–h, the burned area has no significant correlation with PPT, ET and relative humidity.

Overall, the correlations between burned area and TWSC, PPT, ET and relative humidity in the seasonal scale are much stronger than that in the interannual scale. Except for burned area and ET in region S, the burned area has a strong correlation with TWSC, PPT, ET and relative humidity in the seasonal scale. This indicates that climate change has a strong influence on forest fires in Yunnan on a seasonal scale.

5.3.4. Climate Change before and during Extreme Forest Fire Months

To further study the relationship between climate change and forest fires in Yunnan, we compared the DMI, ENSO, PPT, ET, TWSC, water vapor column (Vcol), relative humidity (RH) and burned area in high- and low-fire years. Figure 11 shows the differences in monthly average of DMI, ENSO, PPT, ET, TWSC, Vcol, RH and burned area in high- and low-fire years.

From Figure 11, we found that the change trends of the above seven variances have considerable differences before and during the high-fire months. The differences between DMI and ENSO in high- and low-fire years mainly concentrate in June–October in region NW, while the ones in region S concentrate in July–January. The low PPT level appears from June to December in region NW, while it occurs from July to February in region S. The increase in ET starts in July in region NW, while the ET has no significant change in region S. Correspondingly, the water storage deficit in region NW and S appear from July to January and from September to March, respectively, which is basically the same as the time period with low PPT level. We also found that the water vapor column begins to decrease with water storage deficit. In region NW and S, the low water vapor column starts in September and lasts until February. As the water vapor column decreases, the relative humidity also decreases. The low relative humidity level in region NW and S appears from December to March and from January to March, respectively. This time period basically coincides with the high-fire months in the two regions. The maximum correlation coefficients and lag months between the different variable are shown in

Table 6. The maximum correlation coefficients and lag months between different variables shown in Figure 11.

Variables	Correlation Coefficients		Lag Months
	NW	S	NW	S
DMI vs. PPT	0.65	0.65	0	1
ENSO vs. PPT	0.84	0.49	4	1
PPT vs. TWSC	0.63	0.50	0	1
ET vs. TWSC	−0.30	−0.56	0	0
TWSC vs. Burned area	−0.70	−0.87	5	6
TWSC vs. RH	0.49	0.45	4	5
TWSC vs. Vcol	−0.61	−0.57	2	5
RH vs. Burned area	−0.73	−0.42	4	2
Vcol vs. Burned area	−0.61	−0.38	3	4

.

Combining Figure 11 and Table 6, we plotted the figures on a series of climatic and hydrological changes before forest fires (Figure 12). In region NW, the correlation coefficients between DMI and PPT and ENSO and PPT are 0.65 and 0.84, respectively. This shows that DMI and ENSO have an important influence on PPT. When DMI is negative, as the sea temperature of East Indian Ocean increased, the water vapor transport from the Bay of Bengal weakened, resulting in a decrease in PPT in Southwest China (including Yunnan) [68]. Similarly, when the La Ninã event happened (ENSO is negative), due to the strengthening of East Asian monsoon, the PPT belt moved northward, resulting in less PPT in Southwest China [40,69]. When ENSO events and Indian Ocean Dipoles appeared together, low PPT was strengthened [70]. However, the influence of ENSO on PPT in region NW is much greater than that in region S (correlation coefficient r = 0.84 and 0.49). With low PPT and high ET, TWSC experienced an abnormal decrease, because PPT and ET are the important input and output of the regional TWSC. In region NW, PPT has a stronger correlation (0.63) with TWSC than ET (−0.30). However, the situation is opposite in region S and the correlation coefficients between PPT and ET and TWSC are 0.50 and −0.56, respectively. Four to five months later, the atmospheric water vapor and relative humidity began to decrease, which caused the climatic condition in the forest to trend to be dry, and the water content in the vegetation also decreased [71,72]. This provides powerful conditions for the occurrence and spread of forest fires. The correlation coefficients between relative humidity and atmospheric water vapor and burned area are −0.73 and −0.61 in region NW, respectively. They are −0.42 and −0.38 in region S, respectively. Before the high-fire months, GRACE TWSC has an abnormal change. This explains that TWSC has a strong connection with forest fires (correlation coefficients between TWSC and burned area are −0.70 (NW) and −0.87 (S), lag months are 5 (NW) and 6 (S)). Overall, GRACE TWSC may be useful to reflect the influence of climate change on forest fires in Yunnan.

To further verify the response of GRACE TWSC to climate change before forest fires, we calculated the monthly average of GRACE TWSC and PPT in region NW and S for high- and low-fire years, respectively. From Figure 13 and Figure 14, we found that there is no significant difference between PPT in high- and low-fire years in region NW and S. However, comparing Figure 13a,b, GRACE TWSC shows the different performance in region NW in high- and low-fire years, which in high-fire years is significantly lower than that in low-fire years. This situation is particularly noticeable in the southern region. GRACE TWSC is lower than −9 cm in this region in high-fire years, while it is higher than −6 cm in low-fire years, and this significant difference is more pronounced in region S. Comparing Figure 14a,b, GRACE TWSC is lower than −4 cm in most of region S in high-fire years, and higher than −3.5 cm in low-fire years. Therefore, it can be seen that GRACE TWSC reflects climate change before forest fires in the two regions better than PPT. This is because climate change affects not only PPT, but also other hydrometeorological elements, GRACE TWSC includes soil moisture, runoff, snow water, vegetation canopy water, groundwater, etc. Among them, soil moisture and vegetation canopy water are more sensitive to changes in relative humidity and water vapor in the air [35,73]. Overall, GRACE satellites can detect the influence of climate change on forest fires in Yunnan.

6. Discussion

The changes of forest fires in Yunnan have significant seasonality (Figure 6), and forest fires are mainly affected by wet and dry climates [32]. In Yunnan, the dry season is from November to April of the following year, and the rainy season is from July to September [33]. The forest fires in this region are mainly concentrated in the dry season. The high-fire months in region NW and S are November–February and February–April, respectively. In the dry season, the decrease in PPT causes the ET of forest vegetation to decrease, thus atmospheric water vapor drops to a lower level (Figure 11) [73]. When the atmospheric water vapor is at a low level, the relative humidity is also relatively small (Figure 11). The lower relative humidity level creates a climatic condition more favorable for forest fires [74].

From a geographical point of view, most of the burned areas in Yunnan are mainly concentrated in region NW and S (Figure 5). By analyzing the relationship between climate change and forest fires in the two regions, we found that there is a difference in the influence of climate change on forest fires in the different regions. From Table 3 and Table 4, in region NW, the strongest correlation with burned area is ET (−0.90), followed by PPT (−0.83), relative humidity (−0.79) and GRACE TWSC (−0.57). However, in region S, according to the order of correlation with burned area, the order is relative humidity (−0.92), GRACE TWSC (−0.73), PPT (−0.51) and ET (−0.35). The reason for the above difference between the two regions may be the regional climate types. Region NW belongs to semi-dry and semi-wet climate, and region S belongs to a wet climate [74]. The biggest feature of semi-dry and semi-wet climate is that ET is relatively stronger with less PPT. The huge ET and lower PPT levels cause the water storage of local vegetation to be lower [15,16], resulting in a higher number of combustibles, which in turn increases the possibility of forest fires [25,32]. The wet climate is mainly characterized by humid air and abundant PPT, and the vegetation in region S is dominated by tropical rainforests. Therefore, relative humidity in the air is the main factor affecting forest fires [35]. ET is an important output item of TWSC, while relative humidity has an important influence on vegetation canopy water and soil moisture [71].

Our study found that less GRACE TWSC appeared before forest fires in Yunnan, which is related to climate change. Therefore, GRACE TWSC data may be used to predict the risk of forest fires [35]. From Figure 13 and Figure 14, GRACE TWSC can better reflect local climate change before forest fires than PPT. It can provide a new way of thinking and methods for forest fire prevention.

7. Conclusions

To detect the influence of climate change on forest fires in Yunnan, we used GCTH and least square method to fused TWSC results from six GRACE solutions to improve the accuracy and reliability of TWSC results, and analyzed the relationship between climate factors, GRACE TWSC and burned area. The results show that the temporal distribution of burned area has a significant seasonality, and climate factors and GRACE TWSC have a strong correlation with burned area on the seasonal scale. GRACE TWSC changed significantly in the first few months of forest fires, and it can reflect local climate change. Therefore, GRACE TWSC can well reflect local climate change before forest fires occur.

Our research results will provide a new research method for the influence of climate change on forest fires in Yunnan, and provide strong data support for local governments to formulate forest fire prevention measures.