Variability of fire emissions on interannual to multi-decadal timescales in two Earth System models

Connections between wildfires and modes of variability in climate are sought as a means for predicting fire activity on interannual to multi-decadal timescales. Several fire drivers, such as temperature and local drought index, have been shown to vary on these timescales, and analysis of tree-ring data suggests covariance between fires and climate oscillation indices in some regions. However, the shortness of the satellite record of global fire events limits investigations on larger spatial scales. Here we explore the interplay between climate variability and wildfire emissions with the preindustrial long control numerical experiments and historical ensembles of CESM1 and the NOAA/GFDL ESM2Mb. We find that interannual variability in fires is underpredicted in both Earth System models (ESMs) compared to present day fire emission inventories. Modeled fire emissions respond to the El Niño/southern oscillation (ENSO) and Pacific decadal oscillation (PDO) with increases in southeast Asia and boreal North America emissions, and decreases in southern North America and Sahel emissions, during the ENSO warm phase in both ESMs, and the PDO warm phase in CESM1. Additionally, CESM1 produces decreases in boreal northern hemisphere fire emissions for the warm phase of the Atlantic Meridional Oscillation. Through analysis of the long control simulations, we show that the 20th century trends in both ESMs are statistically significant, meaning that the signal of anthropogenic activity on fire emissions over this time period is detectable above the annual to decadal timescale noise. However, the trends simulated by the two ESMs are of opposite sign (CESM1 decreasing, ESM2Mb increasing), highlighting the need for improved understanding, proxy observations, and modeling to resolve this discrepancy.


Introduction
Extreme fire emissions from Indonesia in 2015 were in the global news and were linked to major degradation of regional air quality (Chisholm et al 2016), which has been suggested to contribute to increased mortality in Southeast Asia (Marlier et al 2012). These fires were brought on by dry conditions associated with the strong 2015-2016 El Niño and, preliminarily, are thought to be the most extreme episode of fires in this region since the 1997-1998 El Niño and could bring even higher economic costs (Chisholm et al 2016).
Factors underlying fires in this region, especially exposure of large amounts of peat following forest clearing, exacerbate the potential severity of the burning and subsequent emissions (Marlier et al 2014, but it is natural climate variability that drives the timing and scale of these events (Chen et al 2011(Chen et al , 2016. Fire variability in other regions has also been connected to the El Niño/southern oscillation (ENSO) (e.g.  Tosca et al 2015), and also play an important role in the variable annual growth rate of CO 2 (Nevison et al 2008).
Despite the clear connection between Pacific sea surface temperature (SST) anomalies and fires in equatorial Asia, our understanding of these relationships on a global scale and for regions without long observational records of fires (including much of the tropics and subtropics) is limited by the length of the satellite record in which active fires can be sensed remotely, roughly 18 years (Giglio et al 2013). Moreover, there are reasons to suspect that global fires exhibit decadal and multidecadal variability because important fire drivers have been shown to vary on these timescales, including precipitation and soil water (Ault et al 2012, Dai 2012, Chikamoto et al 2015, temperature (Meehl 2015), and ENSO itself, which may undergo periods of extreme behavior on a wide range of timescales (Wittenberg 2009, Emile-Geay et al 2013, McGregor et al 2013, Wittenberg et al 2014, Capotondi et al 2015. While the satellite record of fires is too short to fully characterize variability on decadal timescales, charcoal sediment records, used as a proxy for fire emissions, cover much longer time periods but typically with century-averaged values that do not provide information on multidecadal aspects (e.g. Marlon et al 2008, Daniau et al 2012. An alternate proxy for fire activity, burn scars on tree rings, has been used to suggest a connection between decadal climate oscillations and fires in western North America (Kitzberger et al 2007, Taylor et al 2008, Trouet et al 2010. Ice core records of black carbon deposition and trace gas tracers for fire activity (e.g. Zennaro et al 2014) have not been utilized for this purpose to our knowledge, but could provide higher time resolution than charcoal sediments (Bauer et al 2013).
The gap in our knowledge of interactions between climate variability and fires on a global scale spans timescales that are important for near-term predictions, and limits our ability to address questions of detection and attribution. With the lack of available observational data of sufficient length, Earth System models (ESMs), which account for interactions among multiple fire drivers, may be used to provide information about how fires respond regionally and globally to variations in climate across timescales. Numerical experiments with CO 2 -concentration driven ESMs allow a separation between forced signal (e.g. historical warming or land use and land cover change) and internal climate variability. In this study we investigate fires in long model integrations with preindustrial forcing, as well as ensemble simulations of the historical time period using the Community Earth System Model (CESM1) (Hurrell et al 2013, Kay et al 2015 and the Geophysical Fluid Dynamics Laboratory (GFDL) ESM2Mb (Dunne et al 2012, Malyshev et al 2015). We aim to characterize the role of internal climate variability (Baede et al 2001) in driving fires in order to inform future modeling investigations of fire regime changes, extreme events and trend detection.

Methods
The CESM1 Large Ensemble Community Project (LENS) comprises historical and future simulations with many ensemble members, with an aim of improving our ability to distinguish between climate change and climate variability. Kay et al (2015) provide a detailed description of the LENS simulation protocol. In this project, the CESM1 land, ocean, atmosphere and sea-ice components have approximately 1°r esolution and it was forced with 1850 solar and radiative forcing and 1850 land cover for the preindustrial control. We use the last 1800 years of the CESM1 control simulation and also examine the 40 members of the CESM1 historical ensemble from years 1920 to 2005. The ensemble members are initialized with slightly perturbed initial conditions to provide a sampling of internal climate variability under historical climate forcing. Land cover change in the 20th century ensemble is represented by plant functional type (PFT) transitions on a yearly basis and follows Hurtt et al (2011), adjusted to match the CESM1 land model PFT scheme by Lawrence et al (2012). The CESM1 fire model is based on the Thonicke et al (2001) scheme, which simulates fires on a daily basis. In this scheme, probability of fire occurrence is parameterized as a function of soil moisture in the top 0.5 m of the soil, with separate 'moisture of extinction' values for woody and herbaceous PFTs above which fires will not occur. Fire occurrence in this model requires dead litter availability above 100 gC m −2 and also a ground temperature above zero Celsius (Thonicke et al 2001). Fire season length and annual area burned are computed from an empirically-derived fuel moisture function (Thonicke et al 2001). Carbon (C) emissions from fires are determined by applying PFT-specific combustion completeness and mortality factors to the available biomass within the burned area.
ESM2Mb is based on ESM2M (Dunne et al 2012, 2013) with updated parameter settings for the land model LM3 (Malyshev et al 2015), approximately 2°horizontal resolution for the atmosphere and land, and roughly 1°horizontal resolution for seaice and ocean components. The control run was continued for 6000 model years after reaching quasi-equilibrium with solar and radiative forcing representative of year 1860. This simulation used potential vegetation (i.e. undisturbed by human land use) instead of preindustrial vegetation cover and includes dynamic vegetation, important distinctions from the CESM1 control run. In addition to the preindustrial control, we analyze the corresponding 10-member ensemble of historical forcing simulations running from 1861 to 2005 with 20th century land cover changes following Hurtt et al (2011). The ESM2Mb transitions from potential vegetation in the preindustrial control run to 1861 land cover via a 160 year 'bridge' experiment with 1860 radiative forcing and 1700-1860 land cover historical reconstructions (Hurtt et al 2011). The ESM2Mb fire model, described by Shevliakova et al (2009), predicts fires as a function of soil moisture but fires occur annually. Under this scheme, individual months are designated as drought or nondrought, based on the soil moisture deficit. The yearly sum of drought months and the available above ground fuel are used to modulate a PFT-specific climatological fire return-interval factor that yields annual fire emissions of C at each grid point (Shevliakova et al 2009). In contrast to CESM1, this scheme is bound more to the seasonal cycle of soil moisture and fuel availability and has no direct air temperaturedependence.
The fire models used in this study do not capture subtle ecosystem shifts that impact fire dynamics, such as different species make-up within similar PFTs Lightning, however, may be a better predictor of fire on a regional basis (Abatzoglou et al 2016). Moreover, lightning has been shown to vary substantially with tropical Pacific SST anomalies (Sátori et al 2009). We are unable to explore any connections between IAV in lightning and fires with the current model setup. Models that include lightning ignitions typically use a lightning climatology and are, therefore, also unable to explore these connections.
Our analysis of fire emissions from both ESMs is done on an annual basis. Fires in some arid regions respond to the robustness of the previous wet season and associated vegetation growth (Van der Werf et al 2008). We shift the start month of the annual averages for CESM1 emissions at each grid point to eight months prior to the model climatological peak in fire emissions, determined by harmonic analysis, to account for the different regional seasonality. Annual windows for locations with biannual fire cycles were not shifted. This shifting is applied in the ENSO composite analysis only.
We fire emission inventories (details in supplementary material) for comparison to present-day model fire emission variability for 14 regions (figure S1). In addition, we compiled charcoal sedimentation records from the Global Charcoal Database (GCD) for 11 sites in western North America (table S1)  as a proxy for past fire emissions in this region (details in supplementary material). We note that as a metric of fire activity, C emissions emphasize forest fires, while another commonly used measure, area burned, is dominated by savannah and grassland fires (van der Werf et al 2008).
To compute ENSO indices we use SSTs averaged over the NINO3 region (150W-90W, 5S-5N), as in Wittenberg (2009) S2). Precipitation teleconnection patterns are similar between the two models for North America and equatorial Asia but show substantial differences in both the sign and seasonal timing of the response in Africa (figure S3).
We quantify the interdecadal Pacific oscillation (IPO) index as the leading principal component of Pacific basin SSTs after applying a 13 year low-pass filter to the unforced long-control model run (Meehl et al 2009). We split the model timeseries into periods of 300 years before filtering and computing the empirical orthogonal functions. For CESM1, the first EOF explains 39% of the variance and exhibits a spatial pattern very similar to the IPO associated with observed SSTs (not shown) suggesting that the model captures this mode of interdecadal variability. For ESM2Mb, the first EOF explains less than 25% of the variance. Therefore, we do not use the IPO index computed from the ESM2Mb in the remainder of the study. Hereafter we refer to this index as the PDO, noting that the IPO and PDO are highly correlated and the terms are often used interchangeably (Meehl 2015).
The Atlantic Meridional Oscillation index is here defined as the SST anomalies in the North Atlantic Ocean (80W-0W, 0N-60N) minus global SST anomalies (60S-60N), following Trenberth and Shea (2006) as recommended by Phillips et al (2014). We also compute 10 year and 50 year low-pass filtered AMO index timeseries, as in Knight et al (2006), to remove variability related to ENSO and the PDO, respectively.

Results
We assess the IAV in fire emissions globally and by region using the coefficient of variation (CV; quotient of the standard deviation and the mean) (figure 1). Globally, CESM1 and ESM2Mb underpredict the CV in comparison to the GFED4s natural fires. The largest regional biases occur in boreal regions where fires are more episodic than in the tropics (Clark et al 2015). Precipitation is not underestimated in boreal regions (figure S4), suggesting the bias is caused by poor representation of other fire drivers, such as ignition and temperature, which are more important in boreal zones. Temperature, for example, drives boreal fire The power spectrum of CESM1 fire emissions shows substantial variance at periods of 3-7 years, which is significantly different from red noise at a 95% confidence level ( figure 2(a)). The variance of ESM2Mb fire emissions is not clearly distinguishable from red noise at any frequency, although the signal is strongest for periods of 3-7 years ( figure 2(b)). Individual regions show statistically significant variance at these timescales, including temperate North America in both ESMs, and all Asian regions in CESM1 (figures 2(a), (b)). In the CESM1 output, variance is distinguishable from red noise in southern hemisphere Africa and the Middle East on timescales of the PDO (10-30 years), and in Central America and boreal N America on AMO timescales (75-100 years) (figures 2(a), (d)).
The shortness of the GFED4s record of fire emissions makes it difficult to interpret its power spectrum ( figure 2(c)). There is a suggested connection to ENSO with higher variance at a period of 5 years but this periodicity cannot be confirmed statistically with the small sample size. GCD timeseries from western North American sites exhibit strong variability at longer periods of about 40 years and 70 years (figure 2(c)). However, given our knowledge of variability associated with ENSO, it is possible that aliasing of higher frequency oscillations is responsible for some of the variance. The potential for aliasing is difficult to state with certainty given that the GCD data are accumulation rates, each encompassing several years, and have different original timesteps depending on the site. The problematic nature of spectral analysis on the available datasets of fire emissions and proxies stands to highlight the gaps in our knowledge of fire/climate interactions, especially on decadal timescales.
We assess the variability of fire emissions with the aforementioned climate oscillation indices using composite analysis. We group the long control run annual fire emissions by positive or negative index, ENSO, PDO and AMO, and compute the difference in means between the data subsets at each model grid point. A simple statistical test of difference between means is then applied to identify locations where the mean fire emissions shift depending on the sign of the oscillation index. The fire emissions timeseries are standardized such that the differences in mean are plotted as fractions of the standard deviation characteristic of each location. We present analysis of fire emissions relationship to the decadal and multi-decadal climate oscillations for CESM only because ESM2Mb does not generate a reliable PDO and does not produce an AMO with variability on multidecadal timescales ( figure 2(e)).
The spatial pattern of the CESM1 response to ENSO ( figure 3(a)) compares well to the fire emission anomalies observed during the 1997-1998 El Niño with anomalously high emissions in boreal North America, the Amazon, southeast Asia, and Australia, and decreased emissions in northern hemisphere Africa (van der Werf et al 2004)) (figure S5). The ESM2Mb fire emissions exhibit a similar spatial pattern of response to ENSO but the signal is weak ( figure 3(b)), as suggested by the spectral analysis ( figure 2(b)).
The response of CESM1 fire emissions to the PDO index is weaker than the response to ENSO but follows a similar spatial pattern (figure 3(c)) except in boreal North America where fire emissions are negatively correlated with the PDO index. In this region, the average soil moisture increases for warm-phase PDO (not shown), contributing to a decrease in fires during this phase. CESM1 fire emissions also exhibit a statistically significant response to different phases of the AMO in several regions ( figure 3(d)). The abrupt change in sign of the response across the equator in South America, positive southward and negative northward ( figure 3(d)), is consistent with analysis of satellitederived fire counts from Chen et al (2011) who attribute the change to shifts in the ITCZ and precipitation climatology that are associated with the different phases of the AMO. In general, the CESM1 fire emissions response to the AMO is opposite in sign to the PDO response ( figure 1(c)), which is consistent with drought and pluvial relationships between the PDO and AMO (e.g. Findell and Delworth 2010). However, Figure 2. Power spectra of fire emissions (black line) with variance required for statistical difference from red noise with 95% confidence (black dashed line) for (a) the CESM1 long control run, (b) the ESM2Mb long control run and (c) the GFED4.1 s. The spectrum of charcoal influx data from 11 western North America sites interpolated to 20 year time resolution are also included in (c). Regional spectra are shown as gray lines for a) CESM1 and (b) ESM2Mb with frequencies for which the spectra are significantly different from red noise with 95% confidence highlighted by color. Spectra and 95% confidence level curves are plotted for the PDO (brown) and AMO (black) indices for (d) CESM1 and (e) ESM2Mb. There is no PDO plotted for ESM2Mb since the long control run did not produce a reliable PDO pattern.  The ESM2Mb ensemble mean fire emissions increase by about 400 TgC yr −1 between 1920 and 2005 while the CESM1 produces a small decrease (∼25 TgC yr −1 ) in 10 year mean over the same time period (figures 4(a), (c)). The divergence of the trends results from differences in the fire model response to environmental factors, such as temperature and precipitation, and differences in the ESM predictions of these environmental factors. Efforts to determine 20th century trends in fire activity from observations have also led to mixed results. Mouillot and Field (2005) compiled a detailed historical reconstruction of area burned by fires, combining satellite retrievals, field records, and tree-ring data. They conclude that global fires have increased since 1900, with a mid-century minimum in fire activity. However, it is difficult to compare their dataset with the model results reported here which do not account for deforestation fires or human fire ignition and suppression.
CESM1 exhibits a small global trend in fire emissions but there are major regional increases in the boreal northern hemisphere, equatorial Asia and much of the Amazon (figure 4(b)) that appear to be driven by land cover changes. Parts of southeast Asia underwent a transition from largely primary forest vegetation to herbaceous vegetation during the 20th century (Hurtt et al 2011), which reduces the moisture of extinction for fires in the CESM1 fire scheme and also shrinks fuel availability, leading to reduced fire. In boreal forests, increasing temperatures lengthen the fire season in CESM1, enhancing fire emissions despite higher soil moisture.
Historical trends in emissions in ESM2Mb are largely driven by soil moisture changes. Where regional trends are statistically significant from 1920 to 2005 they are uniformly positive ( figure 4(d)). The changes are related to increases in drought months over this time period, and enhanced in the western United States, central Asia, and southern South America by co-located increases in fuel availability. Here we are able to show that historical trends in fire emissions for CESM1 and ESM2Mb, while opposite in sign, are statistically significant at a 95% confidence level with a null hypothesis of zero trend. This is to say that a 20th century anthropogenic signal, including human impacts on climate, CO 2 , and land cover, is detectable in natural fires above the year-to-year noise in fire emissions in these two ESMs. The significance testing used here can be generalized to any size trend and number of ensemble members (table 1).

Discussion and conclusions
The magnitudes of significant trends differ substantially between the CESM1 and ESM2Mb due to  Li et al (2013), which represents peat and deforestation fires, as well as natural and human ignition and fire suppression, and is able to improve model reproduction of fire spatial patterns, total emissions, and IAV (Li et al 2013). Small-scale land cover effects on wildfires, such as ecosystem edge effects and local variations in surface hydrology, are still not well represented in global fire models . These effects could be especially important in regions of intense land conversion, including equatorial Asia where, as mentioned in the Introduction, deforestation can expose massive amounts of peat and increase the sensitivity of regional fire emissions to climate variability . Greater attention to these small-scale, often sub-grid-scale, processes could improve our understanding of anthropogenic impacts on fire variability. The spatial pattern of the fire emissions response to ENSO in both CESM1 and ESM2Mb is roughly consistent with observational records of fires in Canada (Skinner et al 2006) and Alaska (Hess et al 2001), and studies of tree-ring burn-scar synchrony. Greater synchrony, an indicator of increased fire activity (Falk et al 2011), occurs during the cold phase of ENSO in the southwest and interior west United States, and during the warm phase of ENSO in the Pacific northwest (Kitzberger et al 2007, Trouet et al 2010. Burn-scar synchrony studies of connections between fire activity and decadal oscillations, such as the PDO, are generally inconclusive due to uncertainties in the proxies of climate variability (Kipfmueller et al 2012), but possible connections would have major implications for decadal fire prediction on a global scale. For example, best estimates of global fire emissions, such as the GFED4s, are based on data collected almost entirely during a negative phase of the PDO, and there is evidence that a phase change may have occurred in 2014 (Meehl 2015).
Our results suggest that, with this phase change in the PDO, southeast Asia, Australia, the northern Amazon region, and eastern equatorial Africa will see a shift to higher fire emissions contributed by natural climate variability, and the southwestern United States, northern Mexico, and the southwestern Amazon will shift to lower fire emissions contributed by natural variability. Increasing confidence in long-term forecasts such as these will require further progress in our understanding of the natural variability in fires from additional observations and proxy data, as well as improved modeling of global fires.