Earlier peak photosynthesis timing potentially escalates global wildfires

ABSTRACT More intense fire weather due to climate change is implicated as a key driver of recent extreme wildfire events. As fuel stock, the role of vegetation and its phenology changes in wildfire dynamics, however is not fully appreciated. Using long-term satellite-based burned areas and photosynthesis observations, we reveal that an earlier peak photosynthesis timing (PPT) potentially acts to escalate subsequent wildfires, with an increase in the global average burned fraction of 0.021% (∼2.20 Mha) for every additional day of PPT advancement. Satellite observations and Earth System modeling consistently show that this fire escalation is likely due to intensified drought conditions and increased fuel availability associated with the climate feedback arising from earlier PPT. Current fire-enabled dynamic global vegetation models can reproduce the observed negative correlation between PPT and burned area but underestimate the strength of the relationship notably. Given the continued PPT advancement owing to climate change, the bioclimatic effects of vegetation phenology change suggest a potentially pervasive upward pressure on future wildfires.


Satellite-based global vegetation photosynthesis
The recent emergence of solar-induced chlorophyll fluorescence (SIF) serves as a direct indicator of vegetation photosynthesis activity [1][2][3] because it is a by-product of photosynthesis [4].Recent advances in the retrieval of SIF from satellite observation allow for monitoring the seasonal variations of vegetation photosynthesis from regional to global scales [5,6].Thus, we used SIF as a metric to estimate vegetation photosynthesis phenology at different spatiotemporal scales [7][8][9].To explore the robust trends in peak photosynthesis timing (PPT) and its impacts on subsequent wildfire activity, we employed three global reconstructed SIF products from 2001-2018 with diverse spatiotemporal resolutions and methodologies of processing and retrieval: (ⅰ) the clear-sky SIF from the contiguous SIF (CSIF) dataset with 0.05° spatial-and 4-day temporalresolution [10]; (ⅱ) the global 'OCO-2' SIF (GOSIF) dataset with 0.05° spatialand 8-day temporal-resolutions [11]; (ⅲ) and the global long-term contiguous SIF (LCSIF) dataset with a 0.05° spatial resolution and a bi-monthly temporal resolution [12].The three reconstructed datasets are not completely independent SIF products, because all of them were generated from the Orbiting Carbon Observatory-2 (OCO-2) SIF and MODIS reflectance observations and/or meteorological reanalysis data by using machine learning or neural network.Therefore, we further utilized a monthly temporally corrected long-term satellite SIF (LT_SIFc) dataset from 2001 to 2018 with a spatial resolution of 0.05° to track the seasonal variations of photosynthesis.This dataset integrates three satellite-based SIF products, including Global Ozone Monitoring Experiment (GOME), Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY), and GOME-2, correcting the temporal inconsistency of these three products [13].
The order of magnitude of current PPT trends are days, so we selected CSIF with the highest temporal resolution (4 days) as the primary product, with other three SIF (GOSIF, LCSIF, and LT_SIFc) products being used for validation.

Extraction of global peak photosynthesis timing
Before fitting the photosynthesis phenology curves, we first identified the start and end day of year (DOY) for the whole year: From January 1 st to December 31 st for the Northern Hemisphere, and July 1 st in the previous year and June 30 th in current year for the Southern Hemisphere.We used the method of combining Savitzky-Golay (SG) filter and cubic spline [14] to extract PPT and corresponding maximum photosynthesis.This method is flexible to fit a wide range of seasonal variations of vegetation [15], and is effective in capturing the peak photosynthesis timing at large scale [8,14].For each pixel per year, the first step was using the SG filter to smooth the SIF time series.The smoothing window size was set as 7 for CSIF, GOSIF, and LCSIF to get the resultant value positioned at the center of the window.However, given the coarse monthly resolution of the LT_SIFc data, the window size was set as 5 for it.Then, the day of the year with maximum SIF value (SIFmax) was considered as PPT based on the daily SIF reconstructed by cubic spline.Accordingly, we estimated global yearly PPT and maximum photosynthesis from 2002 to 2018 at the native grids of each product.Then, the estimated PPT and maximum photosynthesis were upscaled to 0.25° using the nearest-neighbor and bilinear interpolation methods, respectively.
In addition, to mitigate the influence of abnormally high values in SIFmax, we used the 95th quantile to remove outliers.
Despite being widely used in fire research, the spatial scale of 0.25° may result in the mixed pixel problem that contain different vegetation types and confound the detection of PPT trend and its impact on subsequent wildfires.Therefore, we performed an additional analysis at a finer spatial resolution (0.05°).Consistently, we found an overall advanced PPT at a global scale (Fig. S26a), and a predominantly negative correlation between PPT and subsequent burned area (Fig. S26b).
We calculated long-term trends in PPT at the grid level using the Theil-Sen slope estimator in combination with the Mann-Kendall test.In addition, to quantify the PPT trends for the globe, 4 climate zones, and 11 biomes, we computed the regional median PPT values for each year and used 5-year moving average approach to remove high-frequency variability and potential outliers [16], and estimated PPT trends and corresponding 95% confidence intervals via linear regression.
All SIF observations confirm an overall advancement in PPT at a global scale, with global PPT trends from -0.76±0.40days decade -1 for LT_SIFc to -1.65±0.54days decade -1 for GOSIF (Fig. 2b).The trends of PPT in northern ecosystems (>30°N) were from -0.90±0.44 days decade -1 for LCSIF to -2.32±0.80days decade -1 for GOSIF (Fig. S27), which were consistent with previous studies that reported the trend of PPT in northern ecosystems at -0.83 ~ -1.66 days per decade [8,16,17].The slight discrepancy may arise from the differences in study period, retrieved methodology and data sources.Generally, it demonstrates that our method used to calculate PPT is feasible and reasonable.

Fire perimeters
We utilized fire perimeter products from Canada (National Burned Area Composite, NBAC) and the US (Monitoring Trends in Burn Severity, MTBS) to extract information on burned area (BA) occurring after PPT for the period of 2002-2018.NBAC, developed by the Canada Centre for Mapping and Earth Observation and the Canadian Forest Service, provides yearly burned areas on a national scale since 1986.This dataset is derived from three burned area data sources, including the Canadian National Fire Database (CNFDB) agency data, finer-resolution burned area product from Landsat imagery, and coarse-resolution burned area product from SPOT VEGETATION/Proba-V imagery [18].NBAC provides two layers of information on fire start date: one is based on the first detected hotspot (SDATE) and the other is reported by the fire management agency (AFSDATE).Given no unified burned area mapping method and data source for all fire management agencies across Canada, which may create uncertainties in burned area estimates from AFSDATE [18], we took SDATE as the preferred information, and AFSDATE as the potential source if SDATE is not available.Based on this criterion, the yearly gridded burned area occurring after PPT was derived by three steps: (ⅰ) the fire perimeters were initially rasterized at a spatial resolution of 0.00025°, ensuring good correspondence between fire pixels and perimeters; (ⅱ) the values of SDATE and AFSDATE were allocated to each finer-resolution pixel, respectively; (ⅲ) the finer-resolution pixels with fires occurring after PPT were identified based on the above criterion, and then used to calculate total burned area (expressed as burned fraction with a unit of %) in the 0.25° grid cell.MTBS, developed by the U.S. Geological Survey Center for Earth Resources Observation and Science (EROS) and the USDA Forest Service Geospatial Technology and Applications Center (GTAC), provides national burn severity and extent of large fires since 1984 (https://www.mtbs.gov/).We used the national burned area boundaries dataset (i.e.fire perimeters) from 2002 to 2018.Given the much smaller extent of Hawaii and Puerto Rico, we only considered the continental U.S. and Alaska.This dataset is generated by primarily leveraging Landsat imagery at a 30-m resolution, and Sentinel-2 as a complemental data source since 2015 when there were no available high-quality Landsat observations [19].The dataset provides information on fire incident type (wildfire or prescribed fire) and ignition date.We only considered wildfires, and followed the same processing steps as NBAC to calculate yearly gridded burned areas that occur after PPT.
Accordingly, we obtained yearly gridded burned areas occurring after PPT with a 0.25° resolution in Canada and the US from 2002 to 2018.Note that though NBAC and MTBS datasets have been widely used in assessing fire dynamics and their attributions [20][21][22][23], it should be cautious when using these national-scale fire perimeters data.Because they are generated through completely different approaches and from different data sources, and even within one product some discrepancies could create uncertainties in the estimated burned area boundaries [18,19].
For example, NBAC is derived from three data sources, and there is no unified method and data source for all fire management agency data [18].Thus, we chose the detected hotspot as the preferred data source to partly mitigate the uncertainties from fire management agency data.
Moreover, in the case of MTBS, the minimum burned area mapping units are different in the western and eastern U.S., being 1000 acres and 500 acres, respectively [19].Therefore, we opted to rasterize these different fire perimeter data and conduct our time-series analyses at the grid level.This approach can mitigate potential uncertainties arising from combining these data with large discrepancies.It should be noted that the estimate of the impact of PPT advancement on the following burned area is conservative when using the MTBS because of its omission of small fires.

MODIS global burned area
We also used the MODIS global monthly burned area product (MCD64A1 V6) with 500-m resolution [24] over 2001-2018 to examine the impacts of PPT on subsequent burned areas at the global scale.The product provides information on burn date (expressed as DOY).We calculated the gridded burned area after PPT with a 0.25° resolution for each year through the following steps: (ⅰ) the start and end days for each year were identified, aligning with the global PPT calculation; (ⅱ) any burned pixels located in croplands and non-vegetated areas were masked out based on the IGBP classification scheme from MCD12Q1 V6 land cover product; (ⅲ) the 500-m pixels were allocated into the 0.25° grid cell, with total burned areas occurring after PPT calculated.Note that the pixels repeatedly marked as burned in one year were only recorded once in the burned area calculation.

VPD, CWD, and BUI
To elucidate the mechanistic linkages between PPT and subsequent wildfires, we considered three variables as the mediators: vapor pressure deficit (VPD), representing atmospheric aridity and fire weather [25]; climatic water deficit (CWD), indicating plant water stress [26]; and built up index (BUI), reflecting the total amount of fuel available for combustion [27].Monthly VPD was calculated based on air temperature (T) and dewpoint temperature (Td) at 2 m above the surface following the ref.[28], derived from ERA5-Land monthly data provided at a 0.1° resolution [29].Monthly CWD, calculated as the difference between potential and actual evapotranspiration, was obtained from the TerraClimate dataset with a 1/24° resolution [30].BUI was derived from the fire danger indices historical data produced by the Copernicus Emergency Management Service for the European Forest Fire Information System (EFFIS).This dataset provides daily 0.25° fire danger indices for three different models developed in Canada, the US, and Australia [27].Note that BUI, a key component of the Fire Weather Index System, represents the potential availability of accumulated fuel, thereby also serving as an effective drought indicator.
Increased BUI typically signifies enhanced drought conditions [31].Daily BUI was arithmetically averaged for each month and aggregated over the fire season after PPT.These three variables were upscaled to 0.25° using bilinear interpolation.

Other climatic and auxiliary data
To eliminate the compounding effects of other climatic factors on the post-PPT wildfires, we utilized the climatic variables derived from the ERA5-Land datasets, including 2-m T and Td, total precipitation (PRE), downwards surface solar radiation (RAD), and volumetric soil water content at the 0-7 cm layer (representing surface soil moisture content, SM).Besides, relative humidity (RH) has been recognized as an important variable controlling wildfire dynamics by increasing fire weather [25,[32][33][34], which was calculated with the T and Td based on the ref.[28].
We also used fire radiative power (FRP) to represent fire intensity [35], extracted from the MCD14ML V6 active fire product for the period 2001-2018 [36].This product provides information on hotspot detections, including coordinates, FRP, acquisition time, detection confidence, and fire type.For each year, we allocated the post-PPT vegetation fire hotspots with a detection confidence larger than 50% into 0.25° × 0.25° grid cells with the WGS84 coordinate system, and then averaged FRP to represent FRP per detection after PPT.
To test the relative importance of PPT and biomass available for burning in driving variations of burned area, we used the dataset of global simulated daily net primary productivity (NPP) from Boreal Ecosystem Productivity Simulator (BEPS) over 2001-2018 [37].The BEPS model is driven by remotely sensed leaf area index, clumping index, and land cover type, as well as meteorological and soil data, to simulate daily NPP with a 0.072727° resolution [37].NPP encompasses multiple processes associated with plant biomass accumulation, including carbon sequestration via gross primary productivity (GPP) and carbon loss via autotrophic respiration.Moreover, NPP has been proved to be critical for simulating fire occurrence and behavior [38,39], and is a key output variable representing vegetation productivity, biomass accumulation and subsequent fuel development in state-of-the-art fire-enable vegetation models [40,41].Therefore, pre-PPT accumulated NPP can serve as a proxy for the biomass available for burning.NPP was resampled to 0.25° by using the bilinear interpolation method.Given that SIFmax shows a strong spatiotemporal consistency with pre-PPT accumulated NPP (Fig. S22), it can also be used as an indicator of the biomass available for burning.
Our results were explored for different climate zones (tropical, arid, temperate, and cold), derived from a Köppen-Geiger climate classification map (Fig. S19) [42], and for different biomes derived from the Terrestrial Ecoregions of the World (Fig. S20) [43].We also used the PKU GIMMS NDVI product from combining AVHRR and MODIS observations [44], which provides quality control (QC) layer to check the quality of NDVI value.Thus, we excluded the pixels contaminated by cloud, snow and ice.The pixels with a mean annual NDVI low than 0.1 were regarded as areas with no or sparse vegetation cover, and were excluded from the analysis.

Characterization of climate conditions of potential fire season after PPT
For each year, our analyses were focused on the period from the month of PPT to November (generally from summer to autumn), representing the potential fire season after PPT that can affect the expansion of burned areas (Fig. 1).Therefore, we averaged the above-mentioned climatic variables from the month of PPT to November to characterize the climate conditions of potential fire season, except for precipitation which was accumulated within this period.Therefore, it should be noted that when exploring the causal mechanisms linking PPT and subsequent burned area, we focused on the feedbacks of PPT on climate conditions after PPT.

Analyses
We first calculated year-to-year variations (that is, the difference between two consecutive years) of PPT, subsequent burned area, maximum photosynthesis, and the climatic variables of potential fire season after PPT (including VPD, CWD, BUI, T, PRE, RAD, SM, and RH).This treatment can disentangle the resulting signal from possible long-term dependencies on covariates [45].For each grid cell at a 0.25° resolution, the time series of observations were constructed by combining these year-to-year variations of variables.Additionally, we only considered the grid cell in which fire changes in two consecutive years when constructing the time series.Thus, there are three situations: firefire; fireno fire; no firefire.The analyses were performed based on these constructed time series.However, the structural equation model (SEM) was an exception.
To ensure the model stability, SEM was performed based on the normalized year-to-year variations, calculated by dividing the year-to-year variations by the standard deviation.In addition, to enhance the reliability and better sample the statistical inferences by increasing the sample size, we applied a 9×9 spatial moving window (2.25°×2.25°) to represent the central grid cell during all analyses based on observations.We applied partial correlation to investigate the relationship between PPT and subsequent burned area after removing the effects of SIFmax, T, PRE, RAD, and RH.The impacts of other factors on burned area were also computed to compare with PPT.When calculating the partial correlation between each factor and burned area, we excluded the effects of the remaining factors.
We identified the most important factor for each grid cell based on the absolute value of partial correlation.Since PPT plays a crucial role in shaping vegetation productivity [16], as demonstrated in Fig. S4 where advancing PPT corresponds with increasing SIFmax across large areas of the globe, we also compared the relative importance of PPT and SIFmax (or pre-PPT accumulated NPP) in controlling variations in burned area.Both SIFmax and NPP can represent biomass available for burning (see Other climatic and auxiliary data for details).When comparing PPT with NPP, SIFmax was replaced by NPP in the partial correlation.
Considering the biophysical and biogeochemical feedbacks of vegetation phenology to the climate system and landscape, we utilized VPD, CWD, and BUI as the explanatory factors to explore the causal mechanisms linking PPT and subsequent burned area.First, we used partial correlation to elucidate the effects of PPT on these explanatory variables, and subsequently, the impacts of these variables on burned areas.Because soil moisture plays an importance role in the exchanges of water and energy between land and atmosphere [46], we excluded the effects of temperature, total precipitation, soil moisture and maximum photosynthesis when calculating these partial correlations.
To further explore the underlying mechanisms for the linkage between PPT and subsequent burned area, the SEM was employed.Three pathways were designed: PPT→VPD→BA, PPT→CWD→BA, and PPT→BUI→BA (Fig. S13).When constructing the SEM, we also considered temperature, precipitation, maximum photosynthesis, and soil moisture as the driving factors of VPD, CWD, BUI and burned area.Note that the PPT precedes these climatic conditions (see Characterization of climate conditions of potential fire season after PPT), thus the path is unidirectional.That is, PPT points to VPD, CWD and BUI, representing the feedbacks of PPT to the subsequent climate, but there is no reverse that represents the climate impact on PPT.The SEM was performed at grid level based on the normalized year-to-year variations of time series, and estimated the standard path coefficient for each path to represent the magnitude and direction of the impact.Following the ref.[47], five metrics were selected to assess the goodness of fit of the model, including the Goodness-of-fit Index (GFI ≥ 0.95), Comparative Fit Index (CFI ≥ 0.90), Root Mean Square Error of Approximation (RMSEA < 0.1), Non-Normed Fit Index (NNFI ≥ 0.92), and Standardized Root Mean Square Residual (SRMR < 0.08).The model for each grid cell was deemed reliable when three out of five criteria were met [47].The regional mean standard path coefficient was calculated for the globe, 4 climate zones, and 11 biomes, under the consideration of the goodness of fit of the model (Fig. S28) and the significance level of the path coefficient (pvalue < 0.05).To further quantify the global sensitivity of burned area to PPT, we also used the random forest and explainable machine learning (Shapley Additive Explanation, SHAP) method to strengthen the robustness of the impact of PPT on subsequent burned area than traditional statistical approach (see Random forest and SHAP for details).
We noted that the antecedent climate conditions before PPT may persistently influence those after due to the temporal autocorrelation of climate variables [48,49].This effect may confound the identification of the climate feedback arising from advance PPT.Therefore, we further verified the earlier PPT-induced drought conditions in two aspects.On one hand, we further excluded the effects of pre-PPT (from March to the month of PPT) temperature and precipitation when calculating the partial correlations between PPT and VPD, CWD, and BUI.On the other hand, for each grid cell, we used random forest to model the relationships between year-to-year variations of post-PPT VPD (CWD/BUI) and year-to-year variations of pre-PPT temperature and precipitation.Then, we calculated the VPD * (CWD * /BUI * ) by subtracting the effects of pre-PPT temperature and precipitation estimated by the trained model from the original factors.Finally, the estimated VPD * , CWD * and BUI * were used to develop a SEM with the same structure as Fig. S13.To increase the training sample size, we used a 9×9 spatial moving window to train the model representing the central grid cell.

Model configuration and experimental design
To strengthen the causal mechanisms regarding the climatic feedback of summer vegetation phenology, we further used the Community Earth System Model version 2.2, which is the latest version of the coupled Earth system model developed at the National Center for Atmospheric Research in collaboration with universities and other research institutions [50].We ran the model for AMIP-type simulations with two-way coupled atmosphere (Community Atmosphere Model version 6, CAM6) and land (Community Land Model version 5, CLM5) components, and prescribed sea surface temperature and sea ice concentration of present-day climatology (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005), with biogeochemistry model being inactive.The CAM6 and CLM5 use a nominal 2° (1.9° in latitude and 2.5° in longitude) horizontal resolution with 32 vertical levels and a model top pressure at 3.64 hPa.Satellite-observed vegetation phenology, monthly leaf area index (LAI), stem area index (SAI), and other vegetation features around the year 2000 are prescribed for each plant functional type (PFT) in CLM5.Monthly PFT LAI values were produced based on the 1-km MODIS-derived monthly grid cell average LAI [51,52].The SAI was calculated from the monthly PFT LAI using the method proposed by ref. [53].
To spin up the Earth System model, we first ran it for 50 simulation years starting from the prescribed initial conditions.Then we employed a consistent set of final restart files for the initial conditions in the following 100-year control and sensitivity simulations.We kept all the vegetation parameters unchanged in the control simulation, but modified the prescribed global LAI and SAI in the sensitivity simulation by shifting their summertime (June, July, and August for the Northern Hemisphere; December, January, and February for the Southern Hemisphere) growing phases earlier by 10 days.Specifically, we obtained summertime daily LAI and SAI series by linearly interpolating between monthly values and replaced the original indices with those 10 days later.
Note that the values of LAI and SAI in other seasons remained unchanged.Considering that the annual maximum photosynthesis for the majority of the globe occurs in summer (Fig. S18), the simulated differences of summer-autumn climatic conditions between the sensitivity and control experiments largely represent the climatic feedback of advanced PPT.

FireMIP models
To investigate whether the state-of-the-art fire-vegetation models were able to reproduce the observed impacts of PPT on subsequent burned area, we employed the outputs of FireMIP models to compare with satellite results over 2001-2012 [54].Note that this is the overlapping period (i.e.2001-2012) of satellite observations and FireMIP simulations.The FireMIP aims to evaluate the advanced global fire models and promote projections of global fire dynamics and impacts on ecosystems and human societies under climate change [41].These offline models employ empirical and process-based approaches and incorporate the impacts of climate, vegetation, and human activities on fire dynamics, but no local land-atmosphere feedback.We utilized the seven out of nine FireMIP models (including CLM, JSBACH-SPITFIRE, LPJ-GUESS-SPITFIRE, ORCHIDEE-SPITFIRE, CTEM, JULES-INFERNO, and LPJ-GUESS-SIMFIRE-BLAZE) that output monthly gross primary productivity (GPP) and burned area (Table S2), which allows for the calculation of PPT based on the SG-cubic spline method and subsequent burned area.
We examined the abilities of these models compared with satellite observations from 2001 to 2012 in two aspects: (ⅰ) using partial correlation to qualitatively compare the correlation between PPT and subsequent burned area; ⅱ) utilizing random forest and explainable machine learning (SHAP) method to quantitatively estimate the sensitivity of burned area to PPT (see Random forest and SHAP for details).Both eliminate the compounding effects of maximum photosynthesis, temperature, total precipitation, radiation, and relative humidity.To mitigate the implication of divergent spatial resolutions among FireMIP models and satellite observations, we calculated the area-weighted mean and 95% confidence interval of sensitivities across the globe.

Random forest and SHAP
We used random forest and explainable machine learning (SHAP) methods to estimate the sensitivity of burned area to PPT for the satellite observations and FireMIP models.In contrast to traditional statistical techniques, the methods integrate the advantages of random forest, encompassing bootstrap aggregation and non-distribution assumption, with the strengths of SHAP, ensuring coherence between global interpretations and local explanations [55,56].This method strengthens the robustness of the results from traditional statistics approaches [57].
We conducted this analysis following the ref.[57] through several steps.(ⅰ) For each grid cell, we used the year-to-year variations of burned areas as the dependent variables, and the year-toyear variations of PPT, maximum photosynthesis, temperature, total precipitation, radiation, and relative humidity as the predictors, to train the random forest model.The hyperparameter settings for random forest were number of estimators=100; maximum features=30%; random state=42.(ⅱ) SHAP was used to interpret the trained model to calculate the marginal contributions of PPT on burned areas.(ⅲ) The sensitivity of burned area to PPT was estimated as the slope of the relationship between SHAP-derived marginal contributions of PPT and its year-to-year variations by using the Theil-Sen regression.

Uncertainty from anthropogenic fires
Human activity can modify fire dynamics via ignition, active fire suppression, and manipulating the timing and fuel conditions of fires [58,59].For example, in the United States, human-caused fire can occur throughout the year with a much longer fire season than natural fires, and occur when climate and fuel conditions are not favorable for burning [60], thereby weakening the coupling between climate and fire activity [61].Thus, anthropogenic fires may confound the identification of PPT-fire interactions.However, differentiating between anthropogenic and natural fires remains a challenging task.A recent study compiled a fire-cause reference dataset and applied machine learning technique to predict global patterns of the fractions of fires and burned area related to natural and anthropogenic ignition sources at a large spatial scale.However, due to the limited fire-cause records globally, the fire-cause reference data primarily originated from North America, leading to large uncertainties in other regions.Actually, they did not identify whether a fire is nature-or human-ignited [61].
To partly mitigate the influence of human fire use, we excluded croplands when calculating burned area after PPT within a 0.25° grid cell, as croplands are human-dominated ecosystems (see MODIS global burned area).We further examined the impact of anthropogenic fires from three aspects: (ⅰ) Human activity often leads to land-use and land-cover changes (LULCC), such as landcover conversion for agricultural expansion in African and South American tropical savannas and Asian semi-arid grasslands [59], and deforestation in tropical forests [62], etc.Therefore, we used long-term MCD12Q1 land cover type products to identify the 500-m pixels where LULCC occurred during the study period.These pixels were masked when calculating burned area after PPT (Fig. S23).(ⅱ) We applied the global areas of low human impact dataset with 1 km resolution [63] to distinguish human-and nature-dominated regions.We calculated the fraction of low impact area (LIA) within a 0.25° grid cell.Larger LIA indicates low-impact land by human, where fire activity may be mainly determined by natural factors.Lower LIA indicates areas with larger human impact, where fire activity may be largely influenced by human activity (Fig. S24).(ⅲ) We also employed an intact forest landscape data for the year 2000, which can represent the forest and naturally treeless ecosystems with no remotely detected signs of human activity [64].We compared the effects of PPT on subsequent fire activity among the globe, global intact forest, global areas outside of intact forest, tropical intact forest (25°S-25°N), and extratropical intact forest (Fig. S25).