Importance of Uncertainties in the Spatial Distribution of Preindustrial Wildfires for Estimating Aerosol Radiative Forcing

Uncertainty in preindustrial aerosol emissions, including fires, is one of the largest sources of uncertainty in estimating anthropogenic radiative forcing. Here, we quantify the range in aerosol forcing associated with uncertainty in the location and magnitude of preindustrial fire emissions in a climate model based on four emission estimates. With varied emission location and magnitude among the fire estimates, we find the change in aerosol forcing from present‐day to preindustrial is between −0.4 and 0.3 W/m2 for direct radiative forcing and between −1.8 and 0.6 W/m2 for cloud albedo forcing. Altering only the spatial distribution of preindustrial fires for a fixed magnitude adds a previously unaccounted 25% uncertainty to the total aerosol radiative forcing range. Future studies must account for the uncertainty in the spatial distribution of fire and other aerosol emissions as regional differences contribute substantial additional uncertainty to anthropogenic radiative forcing estimates and the resultant climate sensitivity.


Introduction
Open fires influence the global energy budget by altering surface albedo, via changes in vegetation coverage (Randerson et al., 2006;Scott et al., 2018), and increasing the greenhouse gas burden (Randerson et al., 2006;Ward et al., 2012). Fires also emit aerosols and precursor gases that directly interact with solar radiation (Bowman et al., 2009;Rap et al., 2013) and alter cloud properties and formation processes (Andreae & Rosenfeld, 2008;Tosca et al., 2010). Despite understanding of key fire modeling processes advancing (Hantson et al., 2016) and an expanding network of satellites monitoring burned area, present-day (PD) estimates of fire emissions remain poorly constrained (Carter et al., 2020;Pan et al., 2020;Reddington et al., 2016Reddington et al., , 2019. Given such a high uncertainty in PD fire emission inventories (globally up to a factor of 4; Pan et al., 2020), we can expect an equal or greater uncertainty in preindustrial (PI) estimates (Li et al., 2019), which propagate into high uncertainties in PI-to-PD aerosol radiative forcing (RF) (Hamilton et al., 2018). Understanding and reducing the uncertainty in anthropogenic climate forcing thus requires the development of a more robust PI baseline atmospheric state (Carslaw et al., 2017).
A major component of the variance in PI natural aerosol burden is due to uncertainties in the spatial distribution and total magnitude of fire emissions. A lack of understanding about human influence on fires over the historical period results in a large range of possible land use and land cover change (LULCC) scenarios associated with PI fire estimates (Andela et al., 2017;Hamilton et al., 2018;Kloster et al., 2012;Knorr et al., 2014;Li et al., 2019;Marlon et al., 2008;Ward et al., 2014). Both the AeroCom (Dentener et al., 2006) and Coupled Model Intercomparison Project phase 6 (CMIP6) (Van Marle et al., 2017) estimates Abstract Uncertainty in preindustrial aerosol emissions, including fires, is one of the largest sources of uncertainty in estimating anthropogenic radiative forcing. Here, we quantify the range in aerosol forcing associated with uncertainty in the location and magnitude of preindustrial fire emissions in a climate model based on four emission estimates. With varied emission location and magnitude among the fire estimates, we find the change in aerosol forcing from present-day to preindustrial is between −0.4 and 0.3 W/m 2 for direct radiative forcing and between −1.8 and 0.6 W/m 2 for cloud albedo forcing. Altering only the spatial distribution of preindustrial fires for a fixed magnitude adds a previously unaccounted 25% uncertainty to the total aerosol radiative forcing range. Future studies must account for the uncertainty in the spatial distribution of fire and other aerosol emissions as regional differences contribute substantial additional uncertainty to anthropogenic radiative forcing estimates and the resultant climate sensitivity.

Plain Language Summary
Fires are a key driver of changes to physical landscapes, biogeochemical cycles, and climate on a variety of spatial and temporal scales. In particular, the aerosol particles emitted in fire smoke influence regional and global climate by altering the atmospheric energy budget via interactions with solar radiation and by modifying clouds. In this study, we quantify the aerosol radiative impacts on climate from changing the location and strength of preindustrial fires for four different emission representations. We find that changing only the location of preindustrial fires contributes an additional 25% uncertainty to the estimated range of radiative impacts calculated. Thus, it is important to consider not only changes in the magnitude of fires, but in which regions past fires are occurring when estimating human-induced changes to climate and the Earth System. WAN ET AL. assume a more pristine PI atmospheric state by emphasizing a positive relationship between human population density and fires over the Industrial Era. In contrast, state-of-the-science fire models, including SIMFIRE-BLAZE (Knorr et al., 2014(Knorr et al., , 2016 and LMfire (Pfeiffer et al., 2013), incorporate information from a variety of drivers, including agricultural expansion (Andela et al., 2017), landscape fragmentation, and active wildfire suppression (Kloster et al., 2010), which have historically decreased global burned area (Arora & Melton, 2018) despite population increases. The result of such anthropogenic suppression in global burned area is that PI fire emissions, and the resultant aerosol and gas atmospheric burden, is likely to be similar to, or higher than, PD levels (Hamilton et al., 2018;Rowlinson et al., 2020). Furthermore, when compared to paleoevidence from tree-ring scars, ice cores, and charcoal records (Chellman et al., 2017;Marlon et al., 2016;McConnell et al., 2007;Rubino et al., 2016;Swetnam et al., 2016;Thevenon et al., 2009;Wang et al., 2010;Zennaro et al., 2014), modeling incorporating a higher PI baseline aerosol burden from increased fire activity reflects this paleoenvironmental evidence more robustly than estimates which have a lower PI burden, such as those in CMIP6 (Hamilton et al., 2018;Rowlinson et al., 2020). Modeling including LMfire emission estimates, which explicitly represents human-fire interactions and fire processes in the PI, has been shown to align closely with changes in black carbon (BC) contained within a North American ice core (Hamilton et al., 2018), suggesting that a systematic PD anthropogenic bias can occur in PI emission estimates which are overly constrained by modern human-fire relationships.
Uncertainties in estimated RF from emission inventories and aerosol-cloud interactions (Carter et al., 2020;Hamilton et al., 2018;Rowlinson et al., 2020;Ward et al., 2012) have been found to be more significant than uncertainties due to natural variability in emissions (Hamilton et al., 2018), fire episodicity (Clark et al., 2015), and plume height (Veira et al., 2015). Here, we explore a previously untested aspect of fire emission uncertainty: the importance of altering the fire location for a fixed emission magnitude. We evaluate changes in cloud condensation nuclei (CCN) concentration, BC burden, and the PI-to-PD aerosol direct and indirect RF using the latest released version of the Community Earth System Model (version 2; CESM2). To compare to, and verify the results of, Hamilton et al. (2018), we retain the four PI fire estimates (LPJ-GUESS-SIMFIRE-BLAZE (SIMFIRE-BLAZE) model, LPJ-LMfire (LMfire) model, CMIP6 inventory, and AeroCom inventory) of that study. To isolate the impact of fire location, LMfire, CMIP6, and AeroCom emissions, each representing a plausible PI fire emission map driven by different assumptions of human activity over time is normalized to obtain the globally average total magnitude from SIMFIRE-BLAZE, as this represents the lower modeled PI fire emissions estimate that is consistent with paleofire evidence (Hamilton et al., 2018;Rowlinson et al., 2020).

Fire Emissions and Cases
In this study, we define fire emissions as the sum of natural and anthropogenic vegetation fires. To investigate the PI-to-PD radiative effects of changing the magnitude and spatial distribution of PI fires, we used four climatological (monthly mean) PI fire emission estimates including SIMFIRE-BLAZE (1750-1770 CE), LMfire (1770 CE atmosphere and land use; 150-year simulation), CMIP6 (1845-1855 CE), and Aero-Com (1750 CE) compared against the scientifically validated F2000climo COMPSET, as recommended for CESM2 by the National Center for Atmospheric Research (NCAR), which includes monthly climatological mean emissions for the PD from the CMIP6 data set (1995( -2005. Decadal variability between 1995-2005 and 2006-2014 in CMIP6 latitudinal mean emissions (Figures S1 and S2) is small, so the earlier decade was chosen as it represents the period with higher overall interannual variability (IAV). Since CMIP6 aerosol emissions have low decadal variability and IAV between 1750 and 1900 (Bellouin et al., 2020;Van Marle et al., 2017), 1850 CMIP6 emissions are used in this study to align with the defined PI baseline for the Intergovernmental Panel on Climate Change Sixth Assessment Report (IPCC AR6). We note that small PI IAV in the CMIP6 emissions is potentially unrealistic compared to the PD range, while modeled estimates of PI fire emissions (e.g., LMfire) can produce a PI IAV range more consistent with the PD ( Figure S1). Each fire data set includes horizontally and vertically distributed aerosol number and mass concentration emissions for BC, particulate organic matter (POM; where POM = 1.4 × organic carbon [OC]), and sulfate (SO 4 ) and mass emission of sulfur dioxide (SO 2 ). While injection height depends on a variety of factors including fire intensity and atmospheric conditions, we chose to vertically distribute surface emissions by latitude following Dentener et al. (2006) (Table S2a), as this simple parameterization simulates similar plume heights between neighboring regions produced by more complex vertical methods ( Figure S3) and allows for comparison to similar studies (e.g., Hamilton et al., 2020;Thornhill et al., 2018). However, using a more detailed vertical emission distribution that specifies injection height by biomass burning regions would likely lead to small improvements in simulating smoke transport and estimated radiative impacts, particularly in high-latitude regions where deviations in modeled vertical BC profiles are the largest ( Figure S3), which is consistent with Clark et al. (2015) that suggest fire intermittency is also more sensitive in the high latitudes.
For this study, two sets of PI simulations were conducted: (a) base simulations for each PI data set and (b) normalized simulations where the magnitude of emissions (Equation 1) is equal to SIMFIRE-BLAZE but retains each data set's unique spatial distribution (Table S1 and Figure S4). We also conducted additional simulations with no fire emissions for the PI and PD to isolate the radiative effects from fire aerosols only (e.g., Clark et al., 2015) and a PD simulation with NCAR's CMIP6 vertical fire emissions which assign emission heights by satellite-derived plant unctional Type (PFT) (Val Martin et al., 2018 and Table S2b) to test the RF sensitivity to changing vertical emission scheme. Equation 1. Calculations for each data set's annual average scaling factor (SF; 1.1-1.3) and normalized emissions denoted by LMfire*, CMIP6*, and AeroCom* (1.4-1.6) with spatially varying emissions in longitude (i), latitude (j), and altitude (k).

Model Setup
To model the emission, transport, deposition, and radiative effects of fire aerosols, we used the default released CESM version 2.1 (CESM2.1.0) configuration (Danabasoglu et al., 2020;Liu et al., 2016;Tilmes et al., 2019). The Modal Aerosol Module version 4 (MAM4) in the Community Atmosphere Model version 6 (CAM6) treats aerosols as internal mixtures, separated into four log-normally distributed size modes including Aitken, accumulation, coarse, and primary carbon (Liu et al., 2016), and runs on a 0.9° × 1.25° latitude by longitude horizontal grid with 56 vertical levels. Meteorology was forced for the PI and PD simulations by Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2; Gelaro et al., 2017) from 2006 to 2011 (averaging the last 5 years for analysis). Matching PI and PD meteorology enables isolation of the impact of changes in emissions on aerosol RF using 5-year simulations. A more detailed description of CESM2.1.0 can be found in Table S3.
A variety of metrics were calculated to measure the radiative effects of fire aerosols (Box S1) following the methods described by Ghan et al. (2012) and implemented by Clark et al. (2015). In this study, RF is a term used to describe the direct RF (DRF) and first indirect RF (cloud albedo effect; CAE) at the top of the atmosphere (Myhre et al., 2013). We neglect the second indirect effects and semidirect effects, which tend to be smaller and more uncertain (Forster et al., 2007). Hamilton et al. (2018) used global aerosol microphysics model simulations with PI fire emissions from SIMFIRE-BLAZE and LMfire data sets to predict a 35% and 91% reduction in PI-to-PD CAE, respectively, compared to simulations using CMIP6. Here, we first test the robustness of their methodology (offline aerosol RF calculation) by performing the calculations within a climate model (CESM2). We then present new results from the SIMFIRE-BLAZE normalized cases (LMfire*, CMIP6*, and AeroCom*) to quantify the changes in anthropogenic aerosol RF owing to altering the spatial distribution of PI fire emissions, and thus, climatologically relevant mean PI atmospheric state.

Distribution of Fire Emissions
Pyrogenic aerosol emissions are composed of predominantly carbonaceous material (Andreae, 2019), so the spatial distribution of fire emission is well-characterized by the sum of BC and OC (Figure 1). To quantify regional differences among data sets, emissions within 10 fire regions are compared (Figure 1c). Overall, LMfire carbonaceous aerosol emissions are 1.8 times (67.0 Tg/a) larger than SIMFIRE-BLAZE emissions (36.9 Tg/a), with different regional patterns in the magnitude of emissions (Figure 1a). However, after normalizing the magnitude of LMfire emissions to SIMFIRE-BLAZE (LMfire*), we find a similar relative distribution in emissions between fire model estimates.
At the other end of the emission spectrum, AeroCom emissions (10.3 Tg/a) are less than a third of SIM-FIRE-BLAZE and also exhibit a significantly different regional distribution ( Figure 1a). As PI AeroCom emissions are derived by reducing PD emissions according to changes in population density over time (Dentener et al., 2006), AeroCom's spatial distribution remains similar to the PD CMIP6 pattern ( Figure S4). Aer-oCom emissions are characterized by high emission contribution from Africa (AFRI) and low contribution from TCNA and South America. CMIP6 PI emissions have a similar spatial distribution to AeroCom but a slightly higher total magnitude (15.4 Tg/a) due to CMIP6's inclusion of fire proxy records and modeling (Van Marle et al., 2017). Given AeroCom's poor alignment with paleofire evidence (Hamilton et al., 2018) and its strong dependence on PD observations to estimate PI fire emissions, we exclude AeroCom from calculations of the RF uncertainty range and use CMIP6 as the lower bound for all forcing estimates herein.
Since most anthropogenic emissions occur in the temperate Northern Hemisphere (NH) (Hoesly et al., 2018), we expect anthropogenic aerosol forcing estimates to be most sensitive to changes in the PI baseline within this region (Figure 1b). SIMFIRE-BLAZE (17%) and LMfire* (21%) have similarly high percentages of emissions in the temperate latitudes, which is consistent with ice core and charcoal data (Hamilton et al., 2018;Rowlinson et al., 2020). In contrast, CMIP6* (5%) and AeroCom* (1%) have relatively low emissions in this region owing to disproportionately high emissions from Africa for both fire data sets.

Atmospheric Aerosol State Differences
The CCN formed by pyrogenic aerosol plays an important role in cloud processes (Andreae & Rosenfeld, 2008). Altering the location of fires causes changes in global CAE by enhancing or suppressing cloud formation in different regions (Tosca et al., 2010). The weighted zonal mean CCN concentrations in the base cases ( Figure S5a) have greater variation than the concentrations in the normalized cases ( Figure S5b), suggesting CCN concentrations are less sensitive to spatial changes compared to changes in fire emission magnitude. Additionally, greater variation in CCN when changing absolute magnitude is supported by higher root-mean-squared error differences and statistically significant Kolmogorov-Smirnov (K-S) p-values (p < 0.05) between base estimates compared to the statistical values between the normalized cases (Table S4). Furthermore, CCN concentration distributions for the base cases vary more within the tropics than within the extratropics (Figures S5c and S5d). Variability in CCN concentration between data sets is dominated by modeled differences in CCN concentrations over the ocean while CCN concentrations over land are more similar between data sets (Table S4 and Figures S5a and S5b).
The relative (extra)tropical BC burden for the SIMFIRE-BLAZE (extratropics: 66%, tropics: 34%) and LMfire (extratropics: 67%, tropics: 33%) base cases is slightly higher in the extratropics than for the CMIP6 base case (extratropics: 63%, tropics: 37%) ( Figure S6), following the differing spatial distribution of fire emissions (Figure 1). Percent differences in absolute BC burden relative to SIMFIRE-BLAZE are substantial for LMfire (extratropics: 98%, tropics: 96%) and CMIP6 (extratropics: −58%, tropics: −51%), which are greater than the maximum regional BC burden differences found in Veira et al. (2015) due to varying fire emission height (30%-40%). It should be noted that model parameterizations of deep convection are likely more important than both emission height (Veira et al., 2015) and spatial distribution in realistically characterizing tropical fire BC column burden. Normalized histograms of the modeled BC column burden overlap substantially more than the base histograms, suggesting that the magnitude of BC burden scales proportionally to changes in the magnitude of fire emissions while maintaining similar spatial distributions ( Figure S6). The PD CMIP6, PI CMIP6, and PI CMIP6* cases have similar BC column burden distributions in the tropics and greater variation in the extratropics ( Figure S7).

Sensitivity of Aerosol RF to Preindustrial Fires
DRF measures an aerosol particle's direct interaction with radiation and is important for understanding changes in climate in response to perturbing the atmospheric aerosol burden. The PI-to-PD global average all-sky DRF for each base fire emission case was calculated as follows: −0.40 W/m 2 (LMfire); 0.06 W/m 2 (SIMFIRE-BLAZE); 0.29 W/m 2 (CMIP6); and 0.35 W/m 2 (AeroCom) (Figure 2). LMfire, with the highest PI fire emissions among the examined estimates, exhibits a strong positive (fire) aerosol forcing in the Arctic ( Figure S8), which leads to an overall negative PI-to-PD DRF ( Figure S9), unlike the CMIP6 and SIM-FIRE-BLAZE estimates which have positive global averages (Figure 2). While the CMIP6 and AeroCom cases produce similar spatial distributions of DRF to SIMFIRE-BLAZE, the magnitude of positive forcing over the NH desert regions is larger and no negative forcing is simulated over Greenland ( Figure S9). The DRF spatial distribution of the normalized cases, summarized in Figures 3a and 3b, is similar to that of the each of their respective base cases and dominated by large variations over the Arctic and Northern Africa. The overall PI-to-PD range in DRF (excluding AeroCom) due to uncertainty in PI fire emissions (magnitude and spatial distribution) is 0.70 W/m 2 . Uncertainty in the spatial distribution of PI fires alone contributes an additional 0.12 W/m 2 to the DRF uncertainty range.
CAE is another key quantity in calculating changes in climate that measures the net radiative effect from clouds due to changes in aerosol concentrations. The PI-to-PD aerosol CAE for each base case was calculated as follows: 0.62 W/m 2 (LMfire); −0.50 W/m 2 (SIMFIRE-BLAZE); −1.76 W/m 2 (CMIP6); and −2.17 W/m 2 (AeroCom) (Figure 2). All estimates have strong negative forcing over Southeast Asia due to high anthropogenic emissions in this region ( Figure S10). The SIMFIRE-BLAZE and LMfire cases have similar CAE patterns, dominated by negative forcing in the NH and net positive forcing over the equatorial Pacific Ocean and extending into the Southern Hemisphere ( Figure S10). Conversely, the CMIP6 and AeroCom cases have a different regional CAE pattern dominated by strong cooling in both hemispheres, especially over Southeast Asia, the equatorial Pacific, Brazil, and South Africa ( Figure S10). The CAE spatial distributions for the normalized cases (Figures 3c and 3d) are similar to the base patterns characterized by large variations over the equatorial Pacific and strong cooling in Southeast Asia. For global mean PI-to-PD CAE, the total uncertainty range (excluding AeroCom) from changing emission magnitude and location is 2.38 W/m 2 . Uncertainty in the spatial distribution alone adds an extra 0.65 W/m 2 to the CAE uncertainty range.
When isolating the CAE from fire aerosols only, the standard deviation from fire aerosols (Figure 3f) is identical to the all aerosol plot (Figure 3d) since the variance in both is due to changes only in fire aerosols. The mean CAE for fire aerosols (Figure 3e) also has strong signals over the equatorial Pacific but has weaker forcing over Southeast Asia, stronger forcing over the boreal latitudes, and predominantly positive forcing from PI-to-PD.
Lastly, we find vertically distributing emissions by PFT results in higher PI-to-PD RF magnitudes than vertically distributing by latitude (DRF: +0.01 W/m 2 , CAE: −0.19 W/m 2 ; Figure 2). Compared to the range in DRF and CAE from horizontal spatial uncertainties, RF uncertainty from vertical emissions is smaller, WAN ET AL.
10.1029/2020GL089758 6 of 11  Figure 14-20 from Ramaswamy et al., 2018) between PI BC emission magnitude and the PI-to-PD CAE (blue) and DRF (red). Each point represents one of the four emission data sets (symbol in legend). Dashed lines and transparent symbols show the resultant PI-to-PD CAE and DRF ranges using the PD simulation with emissions vertically distributed by PFT (Table S2b). Black error bars show the additional uncertainty from changing the spatial distribution (fixed SIMFIRE-BLAZE emission magnitude). The bars (right) show the total uncertainty range between LMfire (max) and CMIP6 (min), the black error bar shows the expanded range when AeroCom is included, and the lighter bar shows the range due to spatial distribution. PI, preindustrial; BC, black carbon; PD, present day; CAE, cloud albedo effect; DRF, direct radiative forcing; CMIP6, Coupled Model Intercomparison Project phase 6.
suggesting that horizontal (latitude x longitude) emission uncertainties are likely dominating uncertainty in RF estimates from fires.

Discussion
The absolute value of DRF and CAE uncertainty between our four examined PI fire emission estimates (varying magnitude and spatial distribution) is 0.70 and 2.38 W/m 2 , respectively, scaling almost linearly with the magnitude of emissions ( Figure 2). Moreover, from the linear relationship of forcing to emissions, we can assume that normalizing to a different total magnitude will result in a similar spread around the fitted line. The range in CAE found here is over double the range estimated in Hamilton et al. (2018)   predicted 35% (SIMFIRE-BLAZE) to 91% (LMfire) reductions in CAE compared to estimates from CMIP6, results from this study suggest even larger reductions in CAE between 71% (SIMFIRE-BLAZE) to 135% (LMfire) relative to CMIP6. The high CAE sensitivity to aerosol found here for CESM2 is consistent with recently improved cloud processes in CAM6, compared to CESM1-CAM5, particularly turbulence and ice nucleation parameterizations Danabasoglu et al., 2020;Gettelman et al., 2019Gettelman et al., , 2020; highlighting the importance of repeating experiments with multiple models to assess forcing uncertainty.
Altering only the PI horizontal spatial distribution of fires contributes a nonnegligible amount of additional uncertainty in PI-to-PD estimates of aerosol RF: +17% (0.12 W/m 2 ) to the uncertainty in DRF, +27% (0.65 W/m 2 ) to the uncertainty in CAE, and +25% (0.77 W/m 2 ) to the uncertainty in total RF. Compared to other parameterizations, such as plume height and fire episodicity, the location of emissions constitutes a larger uncertainty in RF. We find that changes in vertical emission distribution contribute +1% (0.01 W/ m 2 ) uncertainty to DRF, +8% (0.19 W/m 2 ) to CAE, and +6% (0.18 W/m 2 ) to total RF. Veira et al. (2015) concluded that plume height parameterization contributes a 0.04 W/m 2 change in total RF, which is even less than the vertical emission uncertainty presented here, while Clark et al. (2015) found episodicity results in a 0.2 W/m 2 change in CAE, around 9% of the overall range in CAE presented in this study.
While the spatial distribution of fires plays an important role in the global energy budget, we must also consider regional impacts, as changes in PI-to-PD forcing between normalized simulations are not uniformly distributed. Spatial uncertainties in fire estimates used in future climate and air pollution projections may not fully capture certain regions' susceptibility to localized climate and human health impacts. For example, arid regions with high fire emissions are more likely to experience drought stress Tosca et al., 2010) and lower income populations are disproportionately impacted by smoke-related mortalities (Johnston et al., 2012;Van Der Werf et al., 2010). Understanding the impacts of spatial heterogeneity in fires on climate and society highlights the importance of better constrained emission data sets.
Beyond fire aerosols, these results have implications for improving current understanding of the spatial uncertainty for all aerosols and on estimates of climate sensitivity. Figure S11 shows a schematic of the RF probability density functions for aerosol cooling and greenhouse gas warming, based on Figure 8.16 from Myhre et al. (2013). Given the narrow curve for greenhouse gas forcing, most of the uncertainty in total anthropogenic RF can be attributed to large uncertainty in aerosol RF (Myhre et al., 2013). Adding the aerosol forcing due to uncertainty in the spatial distribution of PI fires (median absolute change between DRF, CAE, and total RF; ±0.32 W/m 2 per tail) results in a broader and more uncertain distribution in the total forcing over the historical period ( Figure S11). A study examining multiple streams of aerosol RF evidence including modeling, theory, and observations concluded that a PI-to-PD aerosol effective RF more positive than 0.4 W/m 2 is unlikely (Bellouin et al., 2020), which is exceeded here when accounting for the additional forcing uncertainty due to location of PI fire emissions. While Figure S11 is only a crude schematic based on one model, it highlights the need for further research on the impacts of spatial distribution of emissions, in addition to changes in magnitude. Future studies should aim to further constrain the relative importance of spatial distribution in radiative uncertainty by including sensitivity studies with different models, normalized emission magnitudes, and improved vertical fire emission parameterizations. Lastly, multiple studies have shown increased cloud feedbacks in CESM2 result in a stronger surface temperature change and higher climate sensitivity than produced by CESM1 Danabasoglu et al., 2020;Gettelman et al., 2019Gettelman et al., , 2020. Therefore, as models increase CAE sensitivity to aerosol loading , quantifying the PI atmospheric baseline becomes increasingly more important.

Conclusion
Using the spatial pattern within one PI emission estimate neglects the large uncertainty in LULCC over the historical period. We show that estimates of PI-to-PD changes in aerosol RF are likely missing an additional 25% uncertainty in total RF (17% in DRF and 27% in CAE) due to an as yet unconstrained horizontal spatial distribution of PI fire aerosol emissions and 6% uncertainty in total RF (1% in DRF and 8% in CAE) from vertical fire emission uncertainty.
The PI aerosol burden is one of the largest uncertainties in calculating anthropogenic aerosol RF since the Industrial Revolution. Characterizing and potentially reducing the uncertainty in RF from fire emissions could improve our estimates of climate feedbacks between current and PI climates and climate sensitivity. However, more research is needed to better characterize the relationship between land use and population density changes to improve our understanding of the human impact on changing fire regimes with time. Collaboration with experts in the fire ecology community could also further reduce uncertainty in PI emission distributions due to their extensive knowledge about past fire regimes and pyrogeography. Additionally, climate feedbacks within the fire-climate-human system were not quantified in this study. We suggest that future modeling advances could aim to couple fires with climate and human activity to characterize the historical radiative impacts of a different subset of fire feedbacks than presented here. Lastly, we have not considered the radiative impacts from the spatial distribution in other PI natural aerosol species as we assume the uncertainty in fires to be the largest. The results of this study are thus applicable to understanding the radiative impacts caused by spatial uncertainties in other aerosol species in addition to those emitted by fires.