Improved estimates of preindustrial biomass burning reduce the magnitude of aerosol climate forcing in the Southern Hemisphere

See allHide authors and affiliations

Science Advances  28 May 2021:
Vol. 7, no. 22, eabc1379
DOI: 10.1126/sciadv.abc1379


Fire plays a pivotal role in shaping terrestrial ecosystems and the chemical composition of the atmosphere and thus influences Earth’s climate. The trend and magnitude of fire activity over the past few centuries are controversial, which hinders understanding of preindustrial to present-day aerosol radiative forcing. Here, we present evidence from records of 14 Antarctic ice cores and 1 central Andean ice core, suggesting that historical fire activity in the Southern Hemisphere (SH) exceeded present-day levels. To understand this observation, we use a global fire model to show that overall SH fire emissions could have declined by 30% over the 20th century, possibly because of the rapid expansion of land use for agriculture and animal production in middle to high latitudes. Radiative forcing calculations suggest that the decreasing trend in SH fire emissions over the past century largely compensates for the cooling effect of increasing aerosols from fossil fuel and biofuel sources.


Both climate variability and human activity drive changes in wildfire frequency and magnitude. During the past millennium, the human imprint on wildfire has become increasingly important because of landscape fragmentation through land use and, more recently, through large-scale active fire suppression (1). For the time scale relevant to climate change in the industrial era, a key uncertainty is where and to what extent human activity has altered fire activity (2). As a major source of fire ignition, humans use fire for land clearance, thus introducing fire to areas that are unlikely to burn naturally, such as tropical rainforests or peatlands (3). In contrast, recent analyses have suggested that anthropogenic land cover change and landscape fragmentation significantly reduce fire in savannas by affecting fuel load and fire spread, and the fire activity over human-managed land is lower than that under natural conditions (46). For example, the global burned area observed by satellite decreased 24% over the past two decades, mainly driven by agricultural expansion and intensification (5). Before the satellite era, regional and global fire trends have been reconstructed using several types of proxy records, such as charcoal from lake sediments (7, 8), fire-scarred tree rings (9), and chemical impurities or trace gases preserved in ice cores (1014). However, large discrepancies remain among different records, and there is an especially large uncertainty in the trend of fire emissions over the past two centuries. On the global scale, the use of fossil fuel and biofuels has increased and become the major source of carbonaceous aerosols, methane, ethane, and carbon monoxide (CO) in the present day (PD), which may confound interpretation of fire activity from these proxies in ice cores (10, 12, 14). Chemical transport models considering these different sources are therefore needed for the interpretation of ice core records. Dynamic global vegetation models have been used to simulate historical fire emissions. However, different models demonstrate quite different trends of fire activity from the late preindustrial (PI) Holocene to the PD, mainly because of divergent assumptions regarding the response of fire to human demographic growth and to changes in land use and land cover (1517).

Large uncertainty in the PI aerosol loading thus results from uncertainty in PI fire emissions (18). Knowledge of PI aerosol loading is, however, a key for global climate assessments that consider aerosol forcing (19). This forcing typically quantifies the PD aerosol effect on radiative fluxes at the top of the atmosphere (TOA), relative to the aerosol effect in the PI (1750 or 1850 CE). Biomass burning emits both light-absorbing black carbon (BC) and light-scattering organic carbon (OC) aerosols, thus directly influencing the radiative balance via aerosol-radiation interactions (20). Physically and chemically aged smoke particles also serve as cloud condensation nuclei (CCN), consequently altering cloud albedo and indirectly affecting the radiative balance via aerosol-cloud interactions (20). In addition, the cloud albedo forcing of other emissions, such as fossil fuel and biofuel emissions, is highly nonlinear and largely depends on the CCN concentration of the PI baseline, which is, in turn, determined by biomass burning in the PI (19). A recent work by Hamilton et al. (18) suggests that a revised PI biomass burning emission scenario that is consistent with Northern Hemisphere ice core records can reduce the calculated mean global cloud albedo forcing magnitude by 35%, compared to the estimate using emissions prescribed in the Sixth Coupled Model Intercomparison (CMIP6).

In this study, we perform a comprehensive analysis of fire activity and its associated aerosol radiative forcing for the Southern Hemisphere (SH) over the past 250 years. We achieve this by combining an array of Antarctic ice core records of BC deposition, dynamic global vegetation and fire modeling, and atmospheric chemistry transport modeling. We show that BC deposition fluxes in most Antarctic ice cores were roughly constant from 1750 CE to the PD, despite the fact that other anthropogenic emissions—i.e., fossil fuel and biofuel emissions—increased markedly in the SH over the past century (21). This unexpected result can be explained by a large human-induced reduction in wildfire over the same period, as suggested by the fire modeling that we developed independently of the ice core records. The reduced biomass burning emissions largely compensated for the increase in BC emissions from fossil fuel and biofuels.


BC deposition fluxes

We use a collection of 14 new and published BC ice core records retrieved throughout Antarctica (table S1). Figure 1 shows the locations of the ice core records and the BC deposition fluxes from 1750 to the PD, calculated using measured BC concentrations and ice accumulation rates. A recently published Andean ice core record [Illimani, Bolivia, 6300 m above sea level, 16.7°S, 67.8°W] is also included to constrain the fire emissions from the mid to low latitudes of South America (Fig. 2) (22). Although BC is not well mixed on the hemispheric scale, different ice core records represent BC emitted from different source regions across the SH and transported long distances. The simulated contribution of Antarctic BC deposition from each source region is illustrated in fig. S1. The complete array of Antartic BC ice core records can thus provide constraints on the history of fire emissions from the middle to high latitudes across the SH. Compared to charcoal records that are primarily representative of local to regional fire activity (7, 8), the long-range transported BC in ice cores better integrates emissions over large spatial scales (14, 23). A confounding factor, however, is the increasing contribution of anthropogenic fossil fuel and biofuel BC emissions during the time frame (fig. S2) (21). In the present study, these emissions are accounted for in a state-of-the-art chemical transport model GEOS-Chem.

Fig. 1 PI-to-PD variations of BC deposition mass fluxes derived from an Antarctic-wide array of synchronized high-resolution ice core records.

Gray curves represent annual mean BC fluxes from different ice cores. Thick black curves represent the 30-year running average values. Different symbols represent the ice core locations. Age uncertainty is estimated to be less than ±1 to ±3 years.

Fig. 2 Measured and modeled trends of BC deposition fluxes.

(A) Antarctica and (B) Bolivia Andes. The trends are shown as the ratios of historical BC fluxes to PD values (centered around the year 2000). Thick curves in (A) represent the median values of 12 ice core records or colocated modeled results. Shaded areas in (A) represent the 25th to 75th percentile range across sites. The black curve in (B) represents the 30-year running average of BC ratios in an Andean ice core record from Illimani, Bolivia (22). The gray curve represents annually averaged values. Modeled values in both panels are from GEOS-Chem-TOMAS simulations using two different biomass burning emission inventories—biomass burning for CMIP6 (blue) (17) and LPJ-LMfire (orange). Fossil fuel and biofuel emissions are from the CEDS for CMIP6 (21). Shaded areas in (B) represent the SD across different years.

Measured BC deposition fluxes range from 100 to 102 μg m−2 year−1 for different ice core sites across Antarctica (table S1). The high-resolution ice core data show large year-to-year variations (Fig. 1) because fire is episodic in nature, and the peaks in ice cores provide a record of the plumes from large fires transported to Antarctica. However, factors other than fire emissions, such as the interannual variability of atmospheric transport and precipitation, may also contribute to the large variability in BC flux within and across ice cores. To delineate the possible roles of climate variability on the observed BC deposition fluxes, we calculate the correlation coefficients between the BC flux anomalies (i.e., z scores) and several climate indices, including the El Niño–Southern Oscillation NINO3, the Indian Ocean dipole mode index, and the Southern Annular Mode index (fig. S3 and table S2). Although these large-scale climate variations can influence fire weather, atmospheric transport, and precipitation (24, 25), most of the Antarctic ice core records do not significantly correlate with these indices, indicating that the variability of BC fluxes is not driven primarily by climate (table S2). The evaluation of sea salt sodium or water isotope records from the same ice core array indicates no significant long-term changes in atmospheric transport during this period, suggesting that the long-term variability of BC flux in the ice core records is driven principally by the change of emissions from their source regions, likely linked to human activity. For example, several ice core records indicate an increase in BC deposition flux during the early 19th century (Fig. 1), possibly related to land clearing during the early industrial era in the SH. Moreover, most ice core records (10 of 14) show increasing trends after 1970, consistent with the rapid increase of fossil fuel and biofuel emissions quantified in historical emission inventories (21). Even so, we do not fully rule out the possibility that other processes, such as changes in atmospheric circulation and precipitation, can influence the variability of individual ice core records on interannual and decadal time scales. Therefore, our analysis mainly focuses on the long-term historical trend of the ensemble mean of the array of Antarctic ice cores.

Given the complexity of BC sources, atmospheric transport, and deposition, we use a chemical transport model with detailed aerosol microphysics to aid the interpretation of the ice core records [GEOS-Chem-TOMAS (TwO Moment Aerosol Sectional); see Materials and Methods]. BC is short-lived in the atmosphere (~1 week), and the model suggests that its atmospheric concentration can vary several orders of magnitudes during long-range transport. To evaluate our model output, we simulate BC deposition fluxes for the PD (1998 to 2002) using satellite-constrained fire emissions [Global Fire Emissions Database v.4 (GFED4)] (26) along with a state-of-the-art fossil fuel and biofuel emission inventory [Community Emissions Data System (CEDS)] (21). The resulting simulated BC deposition fluxes over Antarctica agree well with the contemporary BC deposition fluxes derived from ice cores (1985 to 2010; fig. S4), indicating that our model can reasonably capture the magnitude and spatial heterogeneity of PD BC deposition over Antarctica. The model simulations further suggest that wet deposition involving precipitation and in-cloud scavenging is the dominant BC loss pathway, contributing 87 ± 8% (means ± 1 SD) of total BC deposition over the model grid cells corresponding to the ice core sites. Dry deposition accounts for the remaining modeled BC flux.

Despite the large temporal variabilities during the past 250 years, average BC deposition fluxes derived from observations for the late PI Holocene period (1750 to 1780) over Antarctica are not significantly lower than those for the PD (1985 to 2010) (Figs. 1 and 2A). The observed PI-to-PD ratio of BC deposition flux is 0.83 ± 0.30 (means ± 1 SD for the eight Antarctic ice cores with records spanning both PI and PD). This ratio is unexpectedly high, given the widely held view that both open biomass burning and fossil fuel/biofuel BC emissions markedly increased in the SH over the past century (17, 27, 28), although some of the increases might occur in the lower-latitude regions to which Antarctic ice cores are less sensitive. This discrepancy is evident in the comparison between ice core observations (Fig. 2A) and the modeled BC deposition trend obtained with standard CMIP6 emission datasets, i.e., BB4CMIP (17) open biomass burning emissions and CEDS (21) fossil fuel and biofuel emissions. Modeled BC deposition with CMIP6 emissions shows a sharp increase from the PI to PD, with a modeled PI-to-PD ratio of 0.23 ± 0.05 (means ± 1 SD) over grid cells colocated with the Antarctic ice cores. This value is less than one-third of the observed value. The modeled increase in BC deposition over Antarctica is caused by both increases in the open biomass burning (BB4CMIP) and fossil fuel/biofuel (CEDS) emissions (Fig. 3). For the PD, BB4CMIP uses exactly the same emissions as GFED4 for the PD, and the simulations with GFED4 and CEDS for PD agree with the ice core measurements (fig. S4). For the PI, emissions from fossil fuel and biofuel sources are negligible. Together, these results suggest that simulations using the CMIP6 fire emissions underestimate PI BC deposition over Antarctica.

Fig. 3 PI-to-PD changes of BC emissions from biomass burning in the SH.

(A) Maps show the annually averaged BC emission fluxes from biomass burning for PI (1750 to 1780) and PD (1998 to 2014). (B) Maps show the PI to PD changes of BC emission fluxes. Zonal mean BC emissions are shown for the PI period in (C) and for the PD in (D). (E) Annual mean BC emissions from biomass burning in the SH from 1750 to 2014. Biomass burning emissions are from the historical biomass burning emissions for CMIP6 (blue, BB4CMIP) (17) and a fire model coupled with a dynamic global vegetation model (orange, LPJ-LMfire) (30). The PD emissions from BB4CMIP are the same as the GFED4 dataset. Fossil fuel and biofuel emissions of BC from CEDS are shown in red in (C) to (E) for comparison. Maps in (A) and (B) are displayed for the SH using south pole stereographic projection to emphasize the spatial variation in emissions across biomass burning regions.

Among the Antarctic ice cores, only the record from James Ross Island at the northern tip of the Antarctic Peninsula shows a large increasing trend of BC flux (Fig. 1), with a measured PI-to-PD ratio of 0.40. In contrast, the PI-to-PD ratio modeled with CMIP6 emissions for this site is just 0.12. The observed BC increase in James Ross Island can be explained by increasing fossil fuel and biofuel emissions: Sensitivity simulations suggest that more than 70% of BC deposition over the island in the PD is from fossil fuel and biofuel sources, while these emissions in the PI are negligible (fig. S1).

We find a similar discrepancy between the model and an ice core record from the SH tropics. Figure 2B shows the trend of BC flux in the Illimani ice core from the central Andes (16.7°S). Model results suggest that the BC deposition in the Illimani core is primarily from open biomass burning and fossil fuel/biofuel emissions in mid- to low-latitude South America, which includes the Amazon rainforest. Emissions from this region are not well captured by the Antarctic ice cores (fig. S1), and thus, the Illimani record complements the Antarctic records. Between 1750 and 1850, the Illimani record shows that the ratio of BC deposition flux relative to PD increased from 0.25 to 0.6 (Fig. 2B). In comparison, the PI-to-PD ratio modeled with CMIP6 emissions is about 0.08 (Fig. 2B), much lower than observed. In addition, the charcoal flux in the Illimani ice core, which indicates local fire activity, exhibits a PI-to-PD ratio of 0.53 (29), in contrast to the 15-fold increase in BB4CMIP biomass burning emissions within a ~25-km radius of Illimani.

Improved fire modeling

Comparison between model simulations and SH ice core proxies suggests that the CMIP6 inventories might underestimate fire emissions in the PI period, especially from higher latitudes in the SH. The historical biomass burning emission inventory for CMIP6 (BB4CMIP) (17) was constructed on the basis of different data sources. For the most recent few decades, BB4CMIP is principally driven by observational datasets including GFED4 and visibility. During the PI era, BB4CMIP is mainly the ensemble mean of six global fire models. Most of these models assume that a fire activity is positively correlated with human demographics at low to moderate population densities. BB4CMIP consequently shows a slight increasing trend in fire emissions throughout the SH (Fig. 3, A and B), mainly in the low latitudes (Fig. 3, C and D), and hemispheric total BC emission increases from 0.6 TgC year−1 in the PI to 1.0 TgC year−1 in the PD (Fig. 3E). Although the Antarctic ice core records are less sensitive to emissions from lower latitudes, the inconsistency between ice core records and modeled trends of BC deposition using BB4CMIP combined with CEDS highlights the need for improvements in fire trends in middle to high latitudes.

To improve our understanding of historical fire emissions in the SH, here we perform a new set of simulations using a process model specifically designed to represent anthropogenic fire in PI time. We use the LPJ-LMfire global dynamic vegetation and fire model (30) to simulate fire emissions from 1750 to 2015 in a transient simulation. LPJ-LMfire represents the human influence on fire through explicit simulation of managed burning carried out by agricultural and pastoral societies, by accounting for changes in the size and geographic distribution of these populations over time and by representing the effects of land use on wildfire behavior through fire protection and landscape fragmentation. Even considering this passive suppression, LPJ-LMfire still overestimates burned area in many regions in the PD as compared with GFED4, such as Europe, Middle East, North America, and South America. To correct this systematic bias, we calculated the ratio between GFED4 burned area and that simulated by LPJ-LMfire for each of the 17 regions (fig. S5) and used these ratios to scale the modeled biomass burning emissions for the entire period from 1750 to 2015 in postprocessing (Materials and Methods). This region-specific correction brings the simulated PD emissions into agreement with GFED4 (Fig. 3, A and D), while keeping the historical fire trend consistent with the original LPJ-LMfire transient simulation. The scaled LPJ-LMfire emissions are used throughout this study.

Between 1750 and the early part of the 20th century, LPJ-LMfire simulates a relatively stable level of fire emissions in the SH, followed by a 30% decrease until about 1990 (Fig. 3E). SH fire emissions are stable during the first ca. 150 years of our simulations primarily because of slow demographic growth and anthropogenic land cover change; meteorology adds interannual variability to the simulation but is not responsible for a secular long-term change. During the 20th century, the simulated strong decline in fire emissions is caused by the rapid expansion of land use for agriculture and animal production in middle- to high-latitude regions, particularly in South Australia, South Africa, and Argentine Patagonia (Fig. 3, A and B, and fig. S6) (31). Anthropogenic landscape fragmentation as a result of the conversion of natural landscapes to agriculture and rangelands is a well-known form of passive fire suppression (4). Expansion of pastures and rangelands for animal production has a particularly strong influence on fire because high densities of livestock reduce the amount of biomass fuels and occur in combination with landscape fragmentation from roads, fence lines, paddocks, farmsteads, and other structures. In the late 20th century and up to the end of our study period (2010), climate variability combined with the influence of increasing atmospheric CO2 concentrations leads to a small increase in biomass in semiarid regions of the SH [e.g., (32, 33)] with concomitant fire emissions.

After the region-specific correction, LPJ-LMfire predicts that biomass burning emissions of BC over the SH decrease from 1.4 Tg year−1 in the PI to 1.0 Tg year−1 at the PD (Fig. 3E). Most of the reductions occur in higher latitudes from which emissions can be more efficiently transported to Antarctica (Fig. 3, C and D). During the same period, fossil fuel and biofuel BC emissions over the SH increase from <0.1 TgC year−1 in the PI to 0.9 TgC year−1 in the PD. Using LPJ-LMfire emissions as input, the modeled PI-to-PD ratio of BC deposition is 0.97 ± 0.18 (means ± 1 SD) over grid cells colocated with Antarctic ice cores, in much better agreement with measurements than that modeled with the CMIP6 emissions (Fig. 2A). The simulation with scaled LPJ-LMfire emissions also improves the comparison with the tropical Illimani ice core. The modeled PI-to-PD ratio is 0.25, consistent with measurements (Fig. 2B).

Implications for aerosol radiative forcing

Here, we calculate the aerosol radiative forcing from all sources (i.e., biomass burning + fossil fuel and biofuels) relative to 1750 at the TOA for the SH using two different biomass burning emission inventories (i.e., BB4CMIP and LPJ-LMfire). We include biomass burning in the forcing calculations because the change in biomass burning emissions is due largely to anthropogenic activities. For both simulations, we use CEDS for anthropogenic fossil fuel and biofuel emissions. The simulations are performed for six time slices between 1750 and 2000. Maps of direct radiative forcing and cloud albedo forcing simulated for the year 2000 relative to 1750 conditions are illustrated in Fig. 4. Overall, the change in total aerosol from the PI to the PD has a cooling effect mainly because of the increase in CEDS anthropogenic fossil fuel and biofuel emissions. The differing trends in the biomass burning emission inventories, however, change the sign of direct radiative forcing in some regions in the SH, such as Australia and the Amazon (Fig. 4, A and B).

Fig. 4 Aerosol radiative forcing calculated for the year 2000 relative to 1750 conditions.

Aerosol direct radiative forcing calculated by the RRTMG radiative transfer model using GEOS-Chem-TOMAS model output simulated with (A) CEDS fossil fuel and biofuel emissions and BB4CMIP biomass burning emissions and (B) CEDS and LPJ-LMfire biomass burning emissions. Cloud albedo forcing simulated with (C) CEDS and BB4CMIP emissions and (D) CEDS and LPJ-LMfire emissions.

Figure 5A shows the time series of the cloud albedo forcing since 1750 owing to changes in both anthropogenic and biomass burning emissions. The shaded area represents the estimated uncertainty considering the error propagation from the variabilities in the input of emissions and meteorology. Simulation with the LPJ-LMfire emissions yields a less negative cloud albedo forcing than that with the BB4CMIP emissions. For the year 2000, the simulation with LPJ-LMfire predicts a mean cloud albedo forcing of −0.33 W m−2 for the SH, compared with a value of −0.52 W m−2 using the BB4CMIP emissions. These results indicate that cloud albedo forcing in the SH is very sensitive to the change in biomass burning emissions. The difference between BB4CMIP and LPJ-LMfire biomass burning emissions can also influence estimates of the PI CCN number concentration in the SH (fig. S7), thus changing the baseline of climate assessments. Under the relatively clean conditions of the PI SH, changes in the CCN number concentration have a greater impact on cloud albedo forcing than they would under the more polluted conditions of the Northern Hemisphere (fig. S8).

Fig. 5 Simulated SH annual mean aerosol radiative forcing since the PI era (1750).

(A) Cloud albedo forcing, and (B) direct radiative forcing. The forcing calculations incorporate changes in biomass burning, fossil fuel, and biofuel emissions. Open markers represent the model calculated values using GEOS-Chem-TOMAS and RRTMG (see Materials and Methods). Solid curves represent interpolated values based on the anthropogenic SO2 emissions from CEDS (red) and the biomass burning OC emissions from LPJ-LMfire (orange) and BB4CMIP (blue). Shading areas represent the uncertainty considering the variability of emissions.

Figure 5B depicts the values of direct radiative forcing due to aerosol-radiation interactions calculated for total aerosol (i.e., including biomass burning, fossil fuel, and biofuel emissions) using different biomass burning emission inventories. To separate the contributions of fossil fuel and biofuel versus biomass burning aerosols, we also show the direct radiative forcing of anthropogenic emissions only (Fig. 5B). The increase of anthropogenic emissions alone from 1750 to 2000 has a direct radiative forcing of −0.05 W m−2. Over the same period, the increase in biomass burning emissions suggested by BB4CMIP has an additional negative forcing of −0.03 W m−2, and the total aerosol direct forcing is −0.08 W m−2. In contrast, the total aerosol direct radiative forcing calculated when using LPJ-LMfire emissions is just −0.02 W m−2, indicating that the positive forcing of decreasing biomass burning largely compensates the negative forcing of the increasing anthropogenic emissions. These results suggest that the difference in biomass burning emissions can dominate the magnitude of aerosol direct radiative forcing in the SH. Even so, the values of direct radiative forcing are generally one order of magnitude smaller than those of cloud albedo forcing, suggesting that the climate impact of biomass burning emissions is primarily caused by the cloud albedo effect.


The continuous high-resolution BC ice core records from Antarctica and the Andes are highly complementary to previously reported Holocene records of fire proxies, such as charcoal records (7, 34) as well as ice core records of CO (10), methane (11, 12), ethane (12), and levoglucosan (35). Compared with other ice core proxies, BC is chemically stable, thus less susceptible to postdeposition processes. A general consensus of these other records is that the global fire emissions may have been relatively high early in the past millennium (1000 to 1500 CE), with a decreasing trend from the Medieval warm period (~1000 CE) to the Little Ice Age (LIA; 1600 to 1800 CE) (7, 1012). However, there are large discrepancies in the reported trends of fire emissions over the past two centuries. Charcoal (7) and CO (10) records suggest that global fire emissions increased from the LIA minimum to a peak in the late 1800s to early 1900s, followed by a 70% decrease in the late 20th century, although the fire modeling performed by van der Werf et al. (2) suggests that such a large reduction seems unlikely. On the contrary, the methane and ethane records suggest a sharp increase in fire emissions over the past 150 years (11, 12). Our SH BC records and fire modeling support our contention that the PI (1750 to 1850) fire emission level is 20 to 30% higher than that of the PD (~2000). Hamilton et al. (18) also suggested greater PI fire activity compared with the PD, based on several BC records from the Northern Hemisphere. These results imply that extensive fire events may not be uncommon before 1900. Even so, our fire modeling suggests that fire emissions can potentially increase in the 21st century, driven largely by the warmer climate and the fertilizing effect of increasing CO2 on vegetation, especially over unmanaged land such as national forests and parks (36). A historical record of area burnt in Australia shows a declining trend in the 20th century (37), followed by a significant increase since 2000, reaching an unprecedented level in 2019 to 2020 (38, 39). These fire trends are roughly consistent with our model predictions.

These records indicate that the CMIP6 biomass burning emissions widely applied to climate models may underestimate SH fire emissions in the late PI era and further affect estimates of contemporary aerosol radiative forcing. With the improved biomass burning emissions presented here, PI-to-PD aerosol forcing (direct radiative forcing + cloud albedo forcing) in the SH changes from −0.61 to −0.35 W m−2, indicating that large uncertainties in aerosol radiative forcing may stem from uncertainties in the historical trend in biomass burning. Similarly, on the basis of ice core records from Greenland, Europe, and North America, Hamilton et al. (18) suggest that the reduction in biomass burning emissions may also occur in the Northern Hemisphere.

As a caveat, most of the ice core records used in this study are from Antarctica, and only one record is from a tropical/subtropical region. Considering the relatively short atmospheric lifetime of BC, we do not rule out the possibility that SH fire emissions may have increased in tropical/subtropical regions (as suggested by BB4CMIP), a trend to which Antarctic ice cores would be insensitive. More paleoenvironment records from these low-latitude regions are needed to better constrain fire emissions on hemispheric to global scales. Because the scaling factor is small (0.1 to 0.2) in SH South America, our approach of keeping the scaling factor uniform over time leads to an estimate of fire emissions that we consider to represent a lower bound. In late PI time, unscaled biomass burning estimates from LPJ-LMfire would likely have been closer to reality, which would imply an even greater reduction in fire emissions in high latitudes over the 20th century. Another priority for future research should be the compilation of additional paleoenvironmental records of fire activity over the 19th and 20th centuries, particularly from the high latitudes of the SH. The radiative forcing calculations we present here consider neither rapid adjustments of atmosphere or surface to aerosol perturbations nor effects on cloud lifetime (20). Another caveat is that the meteorological fields used as input for the fire model and chemical transport model do not consider the climate feedback of the updated fire emissions. These caveats could be addressed in future studies using fully coupled Earth system models.

Accurate estimates of aerosol radiative forcing are also crucial for better understanding the transient climate response (TCR) and equilibrium climate sensitivity (ECS) to increasing CO2 and more accurate projection of future climate change (40). The negative aerosol radiative forcing can, in part, cancel out the positive forcing of increasing greenhouse gases and contribute to the uncertainty of total radiative forcing. An overly large aerosol cooling implies that models might overestimate TCR and ECS to reproduce historical temperature response (41). A recent study using one of the latest-generation CMIP6 climate models (E3SM) suggested that reducing both the magnitudes of negative aerosol radiative forcing and climate sensitivity yields a better agreement with the observed historical record of the surface temperature (42). Ten in 27 of the CMIP6 climate models have an ECS higher than the upper end of the range (1.5° to 4.5°C) estimated by previous generation models. These high ECS values, however, are not supported by paleoclimate constraints (43). Modest aerosol forcing and climate sensitivity values have also been suggested by other observationally based studies (44, 45). Our improved fire emissions may help to bridge the gap between aerosol forcing estimates from current climate model simulations and the constraints from observations.


Ice core measurements

Details of ice cores included in the Desert Research Institute (DRI) Antarctic array are listed in table S1. All ice cores were analyzed using the well-established DRI method for continuous BC measurements (14, 46, 47). In brief, BC concentrations were measured in a continuous stream of ice core meltwater by coupling a nebulizer to a Single-Particle Soot Photometer (SP2; Droplet Measurement Technologies, Boulder, CO, USA). Annual depositional fluxes were determined using net water-equivalent snowfall rates derived from annual layer thickness (48) and measured density profiles following standard procedures. Note that modern ice core–derived accumulation rates over Antarctica (fig. S9) agree well with precipitation rates from the reanalysis dataset MERRA2 (Modern-Era Retrospective analysis for Research and Applications, version 2) used for model simulations. All of the ice core records were mapped from depth to year using either volcanically constrained annual layer counting or volcanic synchronization followed by linear interpolation (23), with non–sea salt sulfur in the West Antarctic Ice Sheet Divide record on the WD2014 age scale (49) serving as the reference.


The LPJ-LMfire dynamic global vegetation model (30) includes processes to simulate natural wildfires from lightning ignition and a unique representation of anthropogenic fire caused by PI societies, which distinguishes the different relationships between humans and fire among hunter-gatherers, pastoralists, and farmers. To run the model, we developed a transient climate scenario by combining climatological mean climate fields that provide a high–spatial resolution baseline with monthly climate anomalies from simulations of the Goddard Institute for Space Studies (GISS) model E2R Earth system model. The GISS meteorology used to drive LPJ-LMfire combined ensemble results from simulations of the past millennium (850 to 1850) (50), historical period (1850 to 2005) (51), and RCP4.5 (2006 to 2100) (52) prepared for CMIP5. We generated monthly anomalies for the period 1701 to 2015 relative to a 1961-to-1990 averaging period. These anomalies were added to the high-resolution climate fields described in Table 3 of (30). Further details of the protocol we used to prepare the input climate datasets were described in (18).

To capture the effects of land use, we merged the KK10 anthropogenic land cover change scenario (53), which ends in 1850, with the HYDE 3.2 dataset (54), following the protocol described in (55). LPJ-LMfire included a representation of managed fire on agricultural and rangelands in which 50% of the litter on 20% of managed land (cropland and rangeland) was burned annually. Wildfire did not otherwise occur on managed land in the model, and wildfire ignitions were caused only by lightning or foraging people [for further details, see (30)].

To perform a PD simulation with LMfire, we made a transient model simulation for the period 1701 to 2015 driven by the climate and land use time series described above and reconstructed atmospheric CO2 concentrations. We performed the LPJ-LMfire simulation continuously over our study period and archived model state variables with monthly resolution. The spatial resolution of the model output was 0.5° × 0.5°. The land cover change used in LPJ-LMfire model is shown in fig. S6. There was a large increase in pasture in the temperate SH (~30°S) starting in 1880s and intensifying through the 20th century. The cropland area also increased significantly in the 20th century. Modern agricultural technology and air quality legislation can largely suppress fire on human-used land. Therefore, cropland or pasture do not burn or burn much less than natural vegetation. This land cover change, along with associated landscape fragmentation, caused a reduction in fire activity for an increase of human-managed land. For example, the rapid expansion of land use for agriculture and animal production in the SH resulted in a strong decline in fire emissions over Australia, South Africa, and Argentine Patagonia in the late 20th century (Fig. 4).

To futher reduce systematic biases in the fire model output, we scaled the output for each region to the GFED4 burned area for the overlapping 1997-to-2015 period. This scaling method was similar to that used in BB4CMIP for the FireMIP model output (17). However, we used burned areas instead of carbon emissions for the scaling, as GFED4 burned areas were from satellite observations and so more accurate. We divided the globe into 17 fire regions, consistent with the definitions in BB4CMIP. To derive fire emissions, we combined the regional scaled dry matter burned for different vegetation types from LMfire with the emission factors from Akagi et al. (56) and Andreae and Merlet (57). These emission factors are consistent with those used in BB4CMIP (17) and GFED4 (26).

GEOS-Chem-TOMAS simulations

We implemented the historical fire [i.e., BB4CMIP (17) and LPJ-LMfire] and anthropogenic [CEDS (21)] emissions into the three-dimensional (3D), global chemical transport model GEOS-Chem version 12 ( The version of GEOS-Chem used in this study had 47 vertical layers and 4° × 5° horizontal resolution and was driven by meteorology from the MERRA2 from the NASA Global Modeling and Assimilation Office. To simulate aerosol number, mass, and size distributions in the atmosphere, we coupled GEOS-Chem with the TOMAS (58) microphysics module. The TOMAS module used 15 size bins ranging from 3 nm to 10 μm for aerosol components of sulfate, organic aerosols, BC, sea salt, and dust (59). TOMAS simulated particle nucleation, coagulation, growth, and wet and dry deposition. In addition to biomass burning and anthropogenic emissions, GEOS-Chem included emissions for biogenic secondary organic aerosols [SOA; Model of Emissions of Gases and Aerosols from Nature (MEGAN) 2.1 (60) with a simple scheme for SOA formation (61)], sea salt (62), volcanic aerosol (Aerocom), and wind-blown dust [Dust Entrainment and Deposition Model (DEAD)] (63). GEOS-Chem-TOMAS was validated against global measurements of aerosol number concentration and size distribution (64).

We performed simulations for six time slices between 1750 and 2000 at 50-year intervals. For each time slice, the simulation was carried out for 5.5 years using 1998 to 2002 meteorology but changing emissions. Applying PD meteorology allowed us to isolate the effect of changing emissions. Results from the first 6 months in each time slice were discarded, and the remaining results were used for analysis. The output of size-resolved aerosol compositions in the atmosphere and deposition fluxes were archived monthly. Modeled results were shown as 5-year averages, and year-to-year variations were used to estimate the uncertainty in natural variabilities in meteorology and emissions. Using the monthly output of GEOS-Chem-TOMAS, we derived the hygroscopicity parameter κ as a volume-weighted mean of each individual species. We used the κ value for each size bin and particle number size distribution to calculate the CCN number concentration at a supersaturation of 0.2%.

Aerosol radiative forcing

The all-sky direct radiative effect and cloud albedo effect were calculated at the TOA based on the monthly size-resolved aerosol output from GEOS-Chem-TOMAS. Details of these calculations are in (65). To derive the direct radiative effect, we first calculated aerosol optical properties, including aerosol optical depth, single-scattering albedo, and asymmetry factor using Mie code (66). We assumed an external mixture of BC with other species in the Mie calculations. To take into account the possible lensing effect of core-shell mixing, we multiplied the calculated absorption coefficient by a constant absorption enhancement factor of 1.5 (67). The calculated optical properties were then used as input to the offline version of the Rapid Radiative Transfer Model for Global climate models (RRTMG). Using the TOA shortwave and longwave fluxes, we calculated the aerosol direct radiative effect as referenced to an aerosol-free condition.

For the cloud albedo effect, we used the sectional representation of cloud activation from Abdul-Razzak and Ghan (68) to calculate cloud droplet number concentrations (CDNCs) for an updraft velocity of 0.3 m s−1. The cloud activation scheme took monthly data of particle number size distribution and size-resolved κ values as input. As a control, a uniform cloud droplet radius of 10 μm was assumed for the PD, consistent with the satellite measurements of the International Satellite Cloud Climatology Project. The effective radii of low- and mid-level clouds (up to 600 hPa) were perturbed by taking the ratio of CDNC between historical and PD to the one-third power. The implied assumption is that aerosol did not change the liquid water path as the definition of the cloud albedo aerosol indirect effect (i.e., the first indirect effect or Twomey effect). The change in radiative fluxes at the TOA due to the change in cloud effective radius was calculated by offline RRTMG. We did not consider the aerosol effect on cloud lifetime. Both direct radiative forcing and cloud albedo forcing were defined as the differences in the PD aerosol effect relative to the PI baseline (1750), considering changes in all aerosols including biomass burning, fossil fuel, and biofuel emissions.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the students and staff of the DRI Ice Core Laboratory, as well as the U.S. and other field teams. The B53 and B40 ice core samples were provided by the Alfred-Wegener-Institut, the James Ross Island and Gomez samples by the UKRI-NERC-British Antarctic Survey, and the Aurora Basin and W10 samples by the Australian Antarctic Division. Funding: This research primarily was funded by the Geosciences Directorate of the NSF under grants AGS-1702814 and 1702830, with additional support from 0538416, 0538427, and 0839093. P.L. acknowledges start-up support from Georgia Institute of Technology. Author contributions: P.L., L.J.M., J.O.K., and J.R.M. designed the study. N.J.C., M.M.A., M.S., and J.R.M. carried out BC measurements using Antarctic ice core samples provided by J.F., R.M., and M.A.J.C. M.S. provided Illimani ice core dataset. J.O.K., P.L., and Y.L. performed LPJ-LMfire simulations and prepared fire emission data. P.L. performed GEOS-Chem-TOMAS simulations. P.L., J.K.K., and J.R.P. performed radiative forcing calculations. P.L., J.O.K., L.J.M., and J.R.M. led the manuscript writing with specific comments and edits from all other coauthors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Datasets of ice core measurements and LPJ-LMfire emissions presented in this study are available at Figshare repository: For code availability, the global 3D chemical transport model GEOS-Chem (including the TOMAS microphysics module) is available at The LPJ-LMfire Dynamic Global Vegetation model is available at The Python codes used in this paper to analyze ice core data, prepare emission inventory, perform radiative forcing calculations, and produce the figures are available at Figshare repository: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article