High-latitude cooling associated with landscape changes from North American boreal forest fires

eScholarship provides open access, scholarly publishing services to the University of California and delivers a dynamic research platform to scholars worldwide. Abstract: Fires in the boreal forests of North America are generally stand-replacing, killing the majority of trees and initiating succession that may last over a century. Functional variation during succession can affect local surface energy budgets and, potentially, regional climate. Burn area across Alaska and Canada has increased in the last few decades and is projected to be substantially higher by the end of the 21st century because of a warmer climate with longer growing seasons. Here we simulated changes in forest composition due to altered burn area using a stochastic model of fire occurrence, historical fire data from national inventories, and succession trajectories derived from remote sensing. When coupled to an Earth system model, younger vegetation from increased burning cooled the high-latitude atmosphere, primarily in the winter and spring, with noticeable feedbacks from the ocean and sea ice. Results from multiple scenarios suggest that a doubling of burn area would cool the surface by 0.23 ± 0.09 °C across boreal North America during winter and spring months (December through May). This could provide a negative feedback to winter warming on the order of 3–5% for a doubling, and 14–23% for a quadrupling, of burn area. Maximum cooling occurs in the areas of greatest burning, and between February and April when albedo changes are largest and solar insolation is moderate. Further work is needed to integrate all the climate drivers from boreal forest fires, including aerosols and greenhouse gasses. Abstract. Fires in the boreal forests of North America are generally stand-replacing, killing the majority of trees and initiating succession that may last over a century. Functional variation during succession can affect local surface energy budgets and, potentially, regional climate. Burn area across Alaska and Canada has increased in the last few decades and is projected to be substantially higher by the end of the 21st century because of a warmer climate with longer growing seasons. Here we simulated changes in forest composition due to altered burn area using a stochastic model of ﬁre occurrence, historical ﬁre data from national inventories, and succession trajectories derived from remote sensing. When coupled to an Earth system model, younger vegetation from increased burning cooled the high-latitude atmosphere, primarily in the winter and spring, with noticeable feedbacks from the ocean and sea ice. Results from multiple scenarios suggest that a doubling of burn area would cool the surface by 0.23 ± 0.09 ◦ C across boreal North America during winter and spring months (December through May). This could provide a negative feedback to winter warming on the order of 3–5 % for a doubling, and 14–23 % for a quadrupling, of burn area. Maximum cooling occurs in the areas of greatest burning, and between February and April when albedo changes are largest and solar insolation is moderate. Further work is needed to integrate all the climate drivers from boreal forest ﬁres, including aerosols and greenhouse gasses.


Introduction
Boreal forests represent one of the key biomes for understanding feedbacks to climate change because of their large area, carbon stocks (McGuire et al., 2002), biophysical imprints Euskirchen et al., 2010), and sensitivities to disturbance regimes (Weber and Flannigan, 1997;Lloyd et al., 2006). These systems experience particularly strong climate changes: northern high latitudes have been warming approximately twice as fast as the global mean during the last century (IPCC, 2007), and are projected to warm by a further 3-8 • C by 2100 (Christensen et al., 2007). One of the primary feedbacks from boreal forests involves the effects of vegetation on surface energy budgets. Vegetation feedbacks are amplified in these high-latitude environments by seasonal snow cover and can influence global climate: on long timescales, vegetation changes may even be necessary for the initiation of ice ages (Gallimore and Kutzbach, 1996;de Noblet et al., 1996;Horton et al., 2010).
Wild fires are a major agent of change for forest composition in boreal North America (Viereck, 1973;Payette, 1992;Chapin III et al., 2006). These fires are typically crown fires that kill most overstory trees and initiate century-long vegetation succession that alters energy absorption and partitioning. Wildfire frequency and extent is controlled by longterm climate (Carcaillet et al., 2001) and short-term patterns of summer high-pressure systems (Macias Fauria and Johnson, 2008) that alter temperature, precipitation, humidity, and wind (Flannigan and Harrington, 1988). Such conditions are projected to result in increased fire weather and burn area during the 21st century. A comprehensive understanding of the potential impacts of altered fire regimes on vegetation distributions, energy budgets, and climate at the continental scale is thus critical for our understanding of high-latitude feedbacks to anthropogenic climate change.
Species and plant functional compositions of North American boreal forests are largely controlled by successional processes after disturbance, including fire, insects, and to a lesser degree, harvesting (Kurz and Apps, 1999). Mature forests consist primarily of flammable evergreen conifers, including black spruce (Picea mariana), white spruce (Picea glauca), and to some extent in Eastern Canada, balsam fir (Abies balsamea) and white cedar (Thuja occidentalis). Immediately after a fire, winter and spring albedos are increased dramatically (up to 0.8) because of higher snow exposure, while on-site deposition of black carbon and charred surfaces reduce summer albedo (Chambers and Chapin III, 2002;Amiro et al., 2006;Jin et al., 2012). Turbulent exchange is decreased because of lower surface roughness, causing an increase in ground temperature and emission of longwave radiation. Thus, despite more energy absorption from low summer albedos, net radiation and heat fluxes are suppressed (Chambers et al., 2005;Liu and Randerson, 2008). Although black carbon has a substantial impact on summer surface energy budgets, its persistence is limited to the first few years as it is washed away and degraded, and grasses, herbs, and shrubs colonize and/or re-sprout (Amiro et al., 1999;Yoshikawa et al., 2003). Re-growing landscapes during the first two decades are subsequently dominated by shrubs and tree saplings that exhibit higher summer albedos than pre-fire mature ecosystems. Winter and spring albedos continue to remain high, although they generally start to decline within a decade or two (Amiro et al., 2006;.
Approximately 20 yr after a fire, forest succession is dominated by trees. Two principal pathways are possible at the plant functional type (PFT) level: (1) self-replacement by evergreen conifers, or (2) relay floristics, in which deciduous broadleaf trees, primarily aspen (Populus tremuloides), birch (Betula spp.), or balsam poplar (Populus balsamifera), dominate and are gradually replaced by conifers after 60-150 yr. With the latter pathway, summer albedo and latent heat remain higher, and sensible heat lower, than in mature forests because deciduous trees have brighter leaves and partition more of the available energy into transpiration during the growing season Liu and Randerson, 2008). Winter and spring albedos gradually decline as the canopy closes and/or deciduous trees, with near-zero leaf area outside the summer months, are replaced by conifers (Amiro et al., 2006;Lyons et al., 2008). Relay floristics therefore extend post-fire trends in albedo and transpiration for many decades, whereas conifer self-replacement expedites the return to pre-fire conditions. While studies have demonstrated that pre-fire forest age and composition, fire severity, climate, topography, aspect, and initial tree recruitment influence tree successional pathway (Mann and Plug, 1999;De Grandpré et al., 2000;Fastie et al., 2003;Johnstone and Ka-sischke, 2006;Johnstone and Chapin III, 2006;Kurkowski et al., 2008;Johnstone et al., 2010b;, the prevalence of these across North America remains largely uncharacterized. Boreal forest fires are projected to increase in frequency and severity by the end of the 21st century, due primarily to longer growing seasons and exacerbated mid-summer drought. Studies using fire weather indices or statistical relationships between burn area and climate variables from climate models generally predict increases in burn area on the order of 40-150 % across Canada and Alaska by 2100 (Flannigan and Van Wagner, 1991;Stocks et al., 1998;Flannigan et al., 2005;Krawchuk et al., 2009;Amiro et al., 2009;Bergeron et al., 2010;Wotton et al., 2010). In an extreme case, Balshi et al. (2009) projected increases of 250-450 %. Indeed, climate change-mediated intensifications of fire regimes are already being observed (Gillett et al., 2004). Since the mid-20th century, burn area has been increasing in Canada (Podur et al., 2002;Stocks et al., 2003) and Alaska (Kasischke et al., 2010), with consequent increases in young deciduous stands (Bond-Lamberty et al., 2009). If fire regimes do follow such trajectories, we should expect a future landscape with younger stands and modified biophysics.
Recent research has identified the main drivers of firemediated biophysical impacts on high-latitude climate in North American boreal forests. Field studies of fire chronosequences have documented typical post-fire trajectories of energy budget constituents, including albedo, net radiation, and fluxes of sensible, latent, and ground heat (Chambers and Chapin III, 2002;Amiro et al., 2006;Liu and Randerson, 2008;McMillan and Goulden, 2008). Studies utilizing remote sensing have quantified trends of post-fire albedo (Lyons et al., 2008;Jin et al., 2012) and leaf habit (deciduous vs. evergreen,  in various regions of Alaska and Canada. Randerson et al. (2006) demonstrated that a set of fires in Alaska resulted in globally averaged radiative cooling. Euskirchen et al. (2009) used a dynamic vegetation and ecosystem model to project future fire regimes, PFT areas, and feedbacks to atmospheric heating in Alaska.
Here, we extend these analyses of boreal forest fire biophysical impacts to the continental scale and simulate spatially and temporally explicit climate impacts within a coupled Earth system model. Generalized succession curves were derived from remote sensing and the literature, and used to drive vegetation distributions and climate under varying burn area scenarios. Results suggest a cooling trend with increased burning, strongest in late winter and early spring, that becomes amplified due to ocean and sea ice feedbacks after an approximate doubling of burn area.

Introduction
This research involved two primary undertakings: (1) the construction of driving datasets, including boreal forest extent, fire frequency, and vegetation succession, and (2) model simulations of fire-driven vegetation changes, continental energy budgets, and high-latitude climate. In the following sections we describe our approaches to these main tasks and to sensitivity analyses used to quantify uncertainty.
CESM's land model, the Community Land Model version 4 (CLM, Oleson et al., 2010), grid cells are composed of urban, glacier, lake, wetland, and natural vegetation fractions derived from the Moderate Resolution Imaging Spectroradiometer (MODIS), the Advanced Very High Resolution Radiometer (AVHRR), and long-term climate and crop datasets (Lawrence and Chase, 2007). Within natural vegetation, there are 16 plant functional types (PFTs), four of which occur in boreal/arctic North America, and bare ground. We defined our domain on this grid to be North American grid cells that contained at least 50 % boreal/arctic vegetation cover within the vegetated fraction, at least 1 % tree cover from MODIS, and at least one fire within the Alaska Large Fire Database (Kasischke et al., 2002) or the Canadian National Fire Database (Amiro et al., 2001;Stocks et al., 2003;Parisien et al., 2006) (NFDBs) (Fig. 1). Based on analyses of remote sensing data (below), within each grid cell a fraction of boreal/arctic vegetation was defined as forest and the compliment as non-forest. The latter was assigned a composition of boreal shrubs and arctic grass based on MODIS land cover over non-forest pixels, and the PFT composition of the boreal forest fraction changed as a function of our vegetation simulations (Sect. 2.3.1).
A variety of land cover datasets can potentially distinguish boreal forests from non-forest at moderate pixel resolution (approximately 300 m-1 km). These datasets are derived from classification schemes using high-resolution training datasets and vegetation indices from satellite-based sensors, including AVHRR (e.g., Loveland et al., 2000;DeFries and Hansen, 2010;Ramankutty and Foley, 2010), MODIS (Friedl et al., 2002), the VEGETATION sensor (Bartholome and Belward, 2005), and the Medium Resolution Imaging Spectrometer (Bontemps et al., 2011). Initial analyses revealed discrepancies among and between datasets at field sites with known land cover, although MODIS performed best (data not shown). We therefore developed a pixel-level classification scheme based on MODIS land (MCD12Q1 year 2003) and tree cover (MOD44B) at 500 m resolution to separate forest from non-forest vegetation. Chronosequences of natural vegetation pixels over fire scar polygons from the NFDBs demonstrated that while MODIS land and tree cover displayed expected post-fire successional patterns, mean tree cover remained above 22 % (Fig. 2). 87 % of unburned grass and shrub pixels contained tree cover under 18 %, and 87 % of burned pixels or unburned forest types contained tree cover over 18 %. We therefore defined forest pixels as those with natural vegetation cover and tree cover greater than or equal to 18 % (Fig. 1).

Historical burn area
Mean historical boreal forest fire return intervals (FRIs) were quantified on the 2 • CLM grid using point versions of the NFDBs (Fig. 1). We defined FRI as the length of time required to burn an area equal to the boreal forest fraction  In (a), a regression line was fit to annual probabilities of re-burn for 0-60 yr from the Canadian National Fire Database and the Alaskan Large Fire Databases (NFDBs). The sensitivity analysis ("SA") line was forced through (0,0) and used in the "nonflammable young stands" sensitivity analysis. (b-d) were developed using MODIS vegetation continuous fields (MOD44B) and land cover (MCD12Q1) over fire scars from the NFDBs. The shaded area in (b) represents ±1 standard error. A logistic curve was fit to EN trees, and a Weibull curve to shrubs. DB = "deciduous broadleaf". The char type in (e) is included as part of the nonvegetated fraction. of a given grid cell. Because burn area prior to 1960 is likely underestimated due to incomplete statistics (Amiro et al., 2001), we used data from 1961-2010 and accounted for missing years in particular Canadian provinces. In relevant southern grid cells, temperate and boreal vegetation were assigned equal relative burn areas. Because forest pixels accounted for 85 % of burn scars for the domain as a whole, 15 % of the total burn area was assigned to non-forest. Boreal forest FRI within a given grid cell was then calculated as: FRI = boreal forest area km 2 mean annual boreal forest burn area km 2 yr −1 , and constrained to have a maximum value of 1000 yr. Empirical evidence suggests that younger boreal forests are less likely to burn due to lower levels of fine fuels and higher fractions of relatively nonflammable deciduous vegetation (Schimmel and Granström, 1997;Cumming, 2001;Lavoie, 2004). To account for this, we quantified agedependent burn probabilities from fire scar polygons in the NFDBs. Polygons were rasterized to 0.005 • (approximately 500 m), and post-fire stand ages were tracked for all burned pixels. Burn probabilities for each stand age were calculated based on the annual re-burning of pixels, and subsequently aggregated for the entire time period. A linear regression was fit to the annual probabilities for 1-60 yr after a fire. Because data coverage becomes increasingly sparse and many forests approach maturity at this time, we assumed burn probability to remain constant after 60 yr (Fig. 2). While this approach could not fully capture the age dependence of fire susceptibility, it qualitatively represented the expected pattern across the continent. We also tested the sensitivity of our domainwide results to this function (Sect. 2.5).

Succession trajectories
Post-fire succession, even in America's boreal forests of low species diversity and high fire severity, varies substantially across sites depending on environmental and disturbance characteristics (Bonan and Shugart, 1989). Revegetation involves complex interactions between fire severity, mortality, reproductive strategies, local seed sources, soil attributes, topography, and micro-climates (Viereck, 1973(Viereck, , 1983Chapin III et al., 2006;Johnstone et al., 2010b). Because of its long-term nature, there exists little quantitative information on succession patterns across large spatial domains. In this study, we attempted to characterize the mean patterns of succession at the PFT level across boreal North America. The ultimate goal was to represent surface energy budget anomalies via fire-mediated changes to vegetation distributions.
We used remote sensing analyses and expert opinion to construct post-fire vegetation succession trajectories (Fig. 2). Because they displayed clear patterns and directly translated into CLM PFTs, we quantified evergreen needleleaf (EN) tree and shrub succession from chronosequences of MODIS land cover over fire scars from the NFDBs. A logistic curve was fit to the EN data, and a Weibull curve to the shrub. Deciduous broadleaf trees were assumed to dominate during mid-succession (approximately 20-60 yr), resulting in a pattern that qualitatively mimics "relay floristics" succession. We prescribed a mixture of non-vegetated ground and herbaceous plants, represented by grasses, during the first few years following fire. To account for the effect of charred surfaces on summer energy budgets, we assigned part of the non-vegetated fraction to a new PFT entitled "char" (Sect. 2.3.2). Although these trajectories were necessary over-simplifications, the validity of applying them to post-fire landscapes across boreal North America was tested in several ways (Sect. 2.4).

Vegetation
A simple stochastic model of fire occurrence and post-fire succession was developed to quantify spatially explicit boreal forest vegetation as a function of fire regime. Within each 2 • CLM grid cell, the boreal forest fraction was divided into sub-grid cells to represent the size of an average fire. Using fires 200 ha or larger from the point version NFDBs, this was approximated as 7000 ha. Although these fires only represent 2-3 % of the total number of fires in Canada, they account for over 97 % of the burn area (Amiro et al., 2001;Stocks et al., 2003). Most grid cells contained a few hundred sub-grid cells.
This simple fire return time and succession model ran on an annual time step. For each sub-grid cell, the annual probability of burning (P burn) was computed from the product of a general "fire probability seed", which was uniform across the grid cell (and initially estimated from the grid cell FRI shown in Fig. 1), and stand age via the age-dependent fire probability function (Fig. 2, Sect. 2.2.2). A random number from 0-1 was selected and, if this was less than P burn, the sub-grid cell burned and its stand age was reset to zero. PFT composition then followed our succession curves until subsequent burning. When the grid cell mean stand age reached steady state, mean FRI was calculated and compared to the desired FRI. The fire probability seed was adjusted until the two fell within 0.2 % of each other.
For the simulation of present-day vegetation, the desired FRIs equaled those for the historical 1961-2010 period (BA × 1, Sect. 2.2.2). FRIs corresponding to multipliers of 0, 2, and 4 on historical burn area (BA × 0, BA × 2, and BA × 4) were used to quantify changes to vegetation distributions due to the absence of fire and with a doubling and quadrupling of burn area. While a transition from BA × 1 to BA × 2 is the most likely scenario given recent projections, BA × 0 and BA × 4 were included to better constrain fire-climate relationships and examine feedbacks. We effectively applied a multiplier of zero for the BA × 0 scenario by assigning PFTs to the maximum age composition from our PFT trajectories.

Land model
We used the CLM version 4 with satellite phenology (CLM-SP) to quantify the impacts of vegetation on surface energy budgets. This version of the CLM uses remote sensing-based estimates of vegetation phenology to prescribe monthly leaf area index (LAI) for a given PFT and simulate fluxes of energy, momentum, and water. To more realistically represent early post-fire water fluxes, we added a multiplicative scalar of 0.05 to bare soil evaporation in non-vegetated surfaces. This was applied only to the fraction of non-vegetated surface that resulted from burning. Recent parameterizations have addressed the CLM's high evaporation bias over bare soil Lawrence et al., 2011), although it remains an outstanding issue (P. Lawrence, personal communication, 2012). The bias was likely exacerbated in the CLM simulations described here because of the presence of organic matter near the surface of many freshly burnt stands, which is assumed to be absent in the CLM's non-vegetated PFT. Additionally, we added a new PFT, "char", to capture the effect of charred surfaces and deposited ash on summer albedo. Char was represented by a plant with zero leaf area, small stature (0.2 m), and low reflectance (visible and nearinfrared equal to 0.09 and 0.2, respectively, compared to 0.16 and 0.39 for tree and shrub stems, and 0.31 and 0.53 for grass stems).
Uncoupled CLM simulations were conducted for 40 yr, driven by our modified vegetation distributions and repeating year 2000 climate from Qian et al. (2006). Although certain aspects of the model take years to reach equilibrium, especially soil water, the CLM displays no internal variability other than that inherited from the external climate forcing datasets. We therefore extracted 40 yr for our analyses.

Climate model
The CESM includes interactive land, atmosphere, ocean, ice sheet, and sea ice models (Gent et al., 2011). Multiple configurations are possible in which particular components are inactive and replaced with observation-derived data. For our purposes, we ran CESM with active land (CLM-SP), atmosphere (CAM5), sea ice (CICE4), and a "slab ocean" model (Bailey et al., 2011). The slab ocean simulates sea surface temperatures by using the mixed layer heat capacity and prescribing horizontal and vertical heat transport consistent with a fully coupled simulation. This allows for feedbacks between the atmosphere, ocean surface, and sea ice without adding the computational expense of the full threedimensional ocean model.
CESM was spun-up with repeating year 2000 climate forcing datasets (top of atmosphere solar radiation, atmospheric CO 2 concentration, and aerosols from volcanic, anthropogenic, and biomass burning sources) for 35 yr using land cover from our historical (BA × 1) vegetation scenario. Four branch runs were then initiated with land surfaces corresponding to BA × 0, BA × 1, BA × 2, and BA × 4. These scenarios therefore varied only in their surface biophysics due to altered North American boreal forest compositions. Simulations ran for 120 yr after branching with the exception of BA × 4, which proceeded for only 80 yr due to limited computing resources and this scenario's stronger surface forcings and climate responses.
Impact analyses focused on mean anomalies to climate variables over the North American boreal forest domain (Sect. 2.2.1), although surrounding areas were assessed for feedbacks and broader effects. Unless stated, summer, fall, winter, and spring refer to JJA, SON, DJF, and MAM time periods, respectively.

Validation
We validated our generalized PFT curves, fire regimes, and land model fluxes in several ways. Firstly, our succession trajectories were input to the CLM to simulate post-fire energy budgets. We compared modeled albedo and energy fluxes during spring and summer to a number of field-and satellite-based observational datasets (Chambers and Chapin III, 2002;Liu et al., 2005;Amiro et al., 2006;Randerson et al., 2006;Liu and Randerson, 2008;Lyons et al., 2008;McMillan and Goulden, 2008). In all cases, the CLM was run uncoupled under repeating year 2000 climate for 40 yr, after which we changed vegetation annually (on 1 January of each year) to match our PFT succession trajectories ( Vegetation distributions in the boreal forest fraction of grid cells were functions of fire frequency and PFT succession from our vegetation model (Sect. 2.3.1). This approach effectively isolates landscape changes due to forest fires, but risks misrepresenting historical vegetation by assuming post-fire succession to be the sole driver of forest composition. Although fire is one of the major determinants, vegetation in boreal North America is controlled by a variety of factors, including climate, solar radiation, soil moisture and temperature, nutrient availability, permafrost, depth of forest floor organic layer, proximity to seed sources, reproductive legacies, and other disturbances (Heinselman, 1978;West et al., 1981;Dyrness et al., 1986;Bonan, 1989;Bonan and Shugart, 1989;Soja et al., 2007). To validate this approach, we compared our simulated present-day vegetation to continental-toglobal land cover datasets, including the CLM input PFT distributions (Lawrence and Chase, 2007), MCD12Q1, Glob-Cover 2009 (Bontemps et al., 2011), GLC2000 North America (Bartholome and Belward, 2005), University of Maryland Land Cover (DeFries and Hansen, 2010), ISLSCP II Potential Natural Vegetation Cover (Ramankutty and Foley, 2010), and the Global Land Cover Characteristics Data Base Version 2.0 (Loveland et al., 2000).

Sensitivity analyses
Four sensitivity studies were conducted to characterize key sources of uncertainty in surface fluxes caused by altered vegetation distributions. Due to the demanding computer resources of the coupled CESM, our sensitivity analyses focused on energy budgets in the uncoupled CLM. In each case, vegetation distributions from BA × 1 and BA × 2 were used, and the CLM was run for 40 yr with repeating year 2000 climate forcing datasets.
In our main simulations, burn area was assumed to increase spatially in direct proportion to historical distributions. This assumption may be problematic for frequently burning grid cells because younger stands are less likely to burn. Additionally, historical data are limited to 50 yr, and the spatial pattern of fires may change in the future. To test the influence of altered burn distributions, we implemented a scenario ("uniform BA increase") in which burn area doubled across the continent but was added uniformly to all boreal forests.
Analysis from the NFDB polygons revealed that nearly all young forest re-burns (< 10 yr) occurred near the perimeters of burn scars (data not shown), where mapping errors were suspected to be greatest. This suggested that the actual probability of re-burn in young stands may be smaller than our regression analysis indicated (Sect. 2.3.2). To address this, our second sensitivity analysis ("nonflammable young stands") used an age-dependent burn probability function that passed through the point (0,0) but was otherwise identical to the main regression.
Despite being indirectly validated in several ways, there remained a high level of uncertainty regarding our succession curves, particularly the contribution of deciduous broadleaf (DB) trees. In our main trajectories, the EN tree fraction was taken directly from a regression of MODIS EN tree pixels over fire scars. We assigned the non-EN tree fraction in mid-successional years exclusively to DB trees. However, in MODIS, this includes DB trees as well as mixed forests and savannas, which undoubtedly contain contributions from conifers. We therefore conducted simulations ("half deciduous tree") in which the DB fraction was halved, and the residual area assigned to EN trees.
In our main simulations, multipliers were applied to boreal forest burn area. There is, however, a negative feedback to changes in burn area due to altered vegetation and stand ages (Fig. 2). Although characterized by Krawchuk and Cumming (2011) for a region in Alberta, Canada, this biotic feedback has been neglected in projections of continent-wide burn area (e.g., Flannigan et al., 1998;Stocks et al., 1998;Amiro et al., 2009;Balshi et al., 2009). We quantified this feedback by conducting a simulation ("biotic BA feedback") in which the × 2 multiplier was applied to the historical fire probability seed, a proxy for the combination of historic fire weather and fuel dynamics, instead of burn area, for each grid cell. This sensitivity study provided an estimate of how much burn area and energy flux anomalies may be dampened by younger, less flammable vegetation under scenarios of increased burning.
The above sensitivity analyses were directed at quantifying uncertainty in surface fluxes caused by altered vegetation distributions from varying fire prescriptions. We also used a Monte Carlo approach to quantify uncertainty arising from low signal-to-noise ratios in our climate simulations (Sect. 3.5). There are, however, numerous other sources of uncertainty and error in our approach that we were unable to address. Prominent among these include the assumption of static forest distributions and succession patterns, biases in our driving vegetation and fire data, and biases in the CLM land and CESM climate models. A more detailed discussion of model caveats is provided in Sect. 4.3.

Validation
After minor parameterization changes and implementation of the "char" PFT, CLM-simulated post-fire energy budgets compared favorably with observations (Fig. 3). A higher level of uncertainty existed for observations and modelobservation comparisons during winter and spring. This was likely due to a greater dependence of winter and spring albedo on burn severity, snow cover, and time that dead boles remain standing. While comparisons with Amiro et al. (2006) implied a high bias in modeled post-fire winter/spring albedo, remote sensing observations from Lyons et al. (2008) suggested model estimates were high for the first decade and low thereafter. Comparisons of monthly albedo with Randerson et al. (2006) were favorable over the entire year for a mature and early burn site, as were latent and sensible heat (data not shown). Considering the variability and uncertainty inherent in post-fire succession and measurement techniques, we believe our approach captured many of the important processes regulating post-fire energy budgets.
Domain-wide vegetation fractions from remote-sensing based datasets displayed a high amount of variability (Table 1). This was to be expected for categorical data derived from different satellite sensors, spectral bands, training layers, analysis years, vegetation classes, classification techniques, and experts' interpretation. Nonetheless, our simulations of historical vegetation compared reasonably well. It must be noted that direct comparisons were not possible due to differences in vegetation classification: observational datasets included mixed forests, woodland/savanna, and/or tundra, all of which were not represented by distinct CLM PFTs. Although our simulated deciduous broadleaf tree coverage was higher than most datasets, deciduous stands are comparatively patchy on the landscape and likely comprise a substantial portion of the mixed forests and woodland/savanna classes present in many of the land cover products shown in Table 1. In our "half deciduous tree" sensitivity analysis, in which the deciduous tree fraction of succession was halved, domain-wide DB tree coverage dropped from 7.9 % to 3.9 %.

Vegetation simulations
In our vegetation model, North American boreal forests varied from nearly pure evergreen needleleaf (EN) stands at long fire return intervals (FRIs) to mixtures of EN trees, deciduous broadleaf (DB) trees, shrubs, grasses, and non-vegetated ground at short FRIs (Fig. 4). DB trees were more dominant than EN trees in stands with FRIs less than approximately 6 yr, and peaked at about 40 % coverage in grid cells with FRIs of about 45 yr. This agreed qualitatively with observations: across Canada, forest inventory data suggest that the   (Loveland et al., 2000). 10 Global Land Cover Database Version 2.0, USGS classification (Loveland et al., 2000).
deciduous fraction of forests varies from 1 % to 40 %, depending on ecozone (Amiro et al., 2001). Modified fire regimes, represented by multipliers of 0, 2, and 4 on historical burn area (BA × 0, BA × 2, and BA × 4), had a pronounced effect on simulated boreal forest vegetation across the domain (Fig. 4). EN tree fraction was approximately halved from 62 % in BA × 0 to 32 % in BA × 4, with DB trees and shrubs accounting for most of converted vegetation. Including tundra pixels, shrub coverage nearly equaled EN tree coverage under BA × 4 (28 % shrubs and 32 % EN trees).

Uncoupled land model simulations
The largest impacts on domain-wide CLM surface energy budgets occurred during late winter and spring (Fig. 5). Increased burning resulted in higher albedos caused by more snow exposure. With increasing daylight as winter progressed to spring, this acted to substantially reduce February-April net radiation (−7.4 % for BA × 2 and −17.3 % for BA × 4) and sensible heat (−7.9 % for BA × 2 and −18.1 % for BA × 4).
Anomalies during summer months were comparatively smaller (Fig. 5) and mainly controlled by differences in energy partitioning and turbulent exchange. With increased burning, domain-wide albedos were slightly elevated. Shorter roughness lengths decreased turbulent energy exchange between the biosphere and atmosphere, which tended to heat the ground surface. Consequently, although less solar energy was absorbed, outgoing longwave radiation increased. Despite this decrease in net radiation, latent heat was still elevated (+1.6 % and +3.6 % for BA × 2 and BA × 4, respectively, from June-August) because deciduous trees partition more available energy into transpiration during the growing season Fall and early winter displayed minimal changes, as solar insolation and snow accumulation were relatively low. On an annual basis, more burning resulted in higher albedos, less net radiation and sensible heat, and slightly more latent heat (Table 2). One would expect this to result in cooling over the domain, strongest in winter and spring, and atmospheric moistening during the summer. However, continental responses depend on both local surface forcings and remote feedbacks involving the atmosphere, ocean, and sea ice.
While surface flux anomalies displayed a dampened response to increasing burn area, temperature changes were linear, if not amplified (Fig. 8). A number of feedbacks were responsible for this. Prevailing westerlies relayed continental temperature anomalies to the North Atlantic and Greenland (Fig. 7), cooling the sea surface and increasing annual sea ice coverage (Fig. 9): mainly polar ice in the summer and fall, and sub-polar ice in the winter and spring. This produced further cooling in eastern Canada. Snow coverage within the domain also increased (+2.6 % and +8.1 % in April-May for BA × 2 and BA × 4, respectively), amplifying the positive land surface albedo anomalies (Fig. 5). Total column atmospheric water dropped slightly with increased burning (from +0.6 % in BA × 0 to −1.1 % in BA × 4 annual mean anomalies), reducing its greenhouse effect. Regressions on the relatively linear responses from BA × 0-BA × 2 for annual, winter-spring, and February-April time periods yielded slopes of −0.12 ± 0.05, −0.24 ± 0.09, and −0.43 ± 0.12 • C per unit increase in historical  burn area. This was equal to −0.06 ± 0.02, −0.11 ± 0.04, and −0.20 ± 0.06 • C/(Mha yr −1 ) ( Table 4). As a result of climate feedbacks and, potentially, model biases, domain-wide surface energy flux anomalies were somewhat dissimilar between coupled and uncoupled simulations ( Table 2). While land surface albedo changes were comparable, the ground surface cooled considerably more with increased burning in the coupled model, which tended to reduce longwave emission and dampen net radiation anomalies. Additionally, when coupled to cooler atmospheres with less evaporative demand, the land surface in BA × 2 and BA × 4 evaporated less water throughout the winter, spring, and early summer (Fig. 5). Although surface-specific humidity decreased during these months, so did temperature, so that spring vapor pressure deficits were reduced by about 3 % in BA × 2 and 7 % in BA × 4 ( Table 3). The CESM also displayed a high temperature bias over the domain during these months (Fig. 6), which may have exacerbated vapor pressure deficit anomalies due to the non-linear Clausius-Clapeyron relationship between saturation vapor pressure and temperature. Annual latent heat changes in BA × 2 and BA × 4 were therefore opposite in sign between the two sets of simulations (Table 2), and sensible heat anomalies were further dampened. As a result of this and other potential dynamical changes, precipitation displayed a weak negative trend with increased burning (Table 3).
The boreal atmosphere responded in other notable ways to fire-induced changes in vegetation composition. Because sensible heat decreased with more burning, planetary boundary layers were lowered throughout the year, with greatest reductions in spring (Table 3). Cooling also created high surface pressure and anti-cyclonic wind anomalies in BA × 2 and BA × 4; the opposite was true for BA × 0. Wind anomalies were greatest in winter, when anti-cyclonic flow dominated, with domain-wide values of −1.3 %, +1.5 %, and +3.4 % for BA × 0, BA × 2, and BA × 4, respectively. A number of atmospheric responses reached beyond the boreal domain into temperate North America and the North Atlantic (e.g., Fig. 7).    1 Values in parentheses represent percent changes from the main simulation differences. 2 Burn area was added evenly to the boreal forest fraction of all grid cells, instead of in proportion to historical burn area, as was done in the main simulations. 3 The probability of fire as a function of stand age was forced to pass through the point (0,0). 4 The original deciduous broadleaf tree fraction during succession was halved, and the remainder assigned to evergreen needleleaf trees. 5 Multipliers were applied to historical "fire probability seeds", instead of burn area. 6 W m −2 .  includes vegetation from tundra pixels. Not shown are contributions from temperate vegetation, glaciers, wetlands, lakes, and urban areas.

Uncertainties
Sources of uncertainty, bias, and fallible assumptions in a study such as this are numerous. In an effort to evaluate their impacts, we quantified two major sources of uncertainty: that associated with the construction of vegetation distributions, and that from the large variability in highlatitude climate and consequent low signal-to-noise ratios produced by our surface perturbations. For the former, we conducted four sensitivity analyses with CLM that modified varying components of the BA × 1 and BA × 2 simulations (Table 5): (1) burn area was added uniformly across the domain ("uniform BA increase"), (2) the burn probability vs. stand age function was forced through the point (0,0) ("nonflammable young stands"), (3) the DB tree composition throughout succession was halved ("half deciduous tree"), and (4) multipliers were applied to historical "fire probability seeds" instead of burn area ("biotic BA feedback"). Because net radiation and surface temperatures responded relatively linearly from BA × 0-BA × 2 (Fig. 8, Table 2), we used the variance in net radiation anomalies between BA × 1 and BA × 2 as an estimate of uncertainty resulting from our vegetation distributions (Table 4). This was done separately for each seasonal period. Of the four experiments, "uniform BA increase" impacted surface energy budgets to the greatest degree. This scenario amplified fire-induced annual net radiation and sensible heat changes by over 30 % (Table 5), suggesting that new patterns of burning may strengthen the land surface cooling feedback from fires. The "nonflammable young stands" scenario also amplified energy budget anomalies, but to a much smaller degree. Alternatively, "half deciduous tree" and "biotic BA feedback" dampened energy budget changes in a nearly identical manner, resulting in an approximate 17 % reduction in net radiation anomalies. In "biotic BA feedback", burn area was reduced by 21.7 % in BA × 2 and by 30.5 % in BA × 4.

BAx4
BAx0 Fig. 6. Mean monthly historical (BA × 1) and changes to 2-m air temperature. The shaded area in (a) delineates monthly standard error, and shaded regions in (b-d) represent monthly 95 % confidence intervals. The uncoupled time series is driven by Qian et al. (2006) reanalysis data and used to assess temperature bias in the coupled model.
The CESM coupled climate model displays a considerable amount of high-latitude inter-annual and interdecadal variability because of long-term oscillations in atmosphere-ocean dynamics (Jahn et al., 2012;Landrum et al., 2012). This is true even with the slab ocean model , and is especially prominent in winter: we found the standard error of linearly de-trended surface temperature for land north of 45 • N in CESM for BA × 1 to be 2.37 • C for December-February vs. 1.97 • C in Qian et al. (2006) re-analysis data. Because of limited computing resources, we were able to run our equilibrium climate simulations for a maximum of 120 yr each, with the exception of BA × 4, which we ran for 80 yr. This resulted in somewhat low signalto-noise ratios for many climate variables. To quantify the uncertainty in our relationship of domain-wide temperature vs. burn area due to these limited run times, we performed 1 × 10 6 Monte Carlo simulations of the estimated slope from BA × 0-BA × 2. For each seasonal period, a population of (temperature) was created from the annual values of seasonally averaged temperature. Simple linear regressions were fit to three points chosen from the normally distributed populations of (temperature) in BA × 0, BA × 1, BA × 2 (Fig. 8). In this way we were able to use 360 simulation years from three experiments to estimate the mean and variance of modeled temperature change associated with a doubling of burning.
These two major sources of uncertainty for surface temperature effects, i.e. vegetation distributions and climate model signal-to-noise ratios, were pooled to provide an overall estimate for our approach. We calculate uncertainties of 41, 37, and 28 % for the slopes of temperature vs. burn area during annual, winter-spring, and February-April time periods (Table 4). There are, however, other uncertainties due to fundamental assumptions and potential model biases that we were unable to address directly (Sect. 4.3).

Background and comparisons to previous work
Earth's climate-vegetation system is particularly sensitive to change at high latitudes. Temperature changes alter spatial patterns of plant establishment, survival, and competition. Once the potential for forest is realized, climate-dependent disturbance regimes place strong controls on plant functional composition. Snow exposure amplifies albedo effects from structural vegetation changes. The relatively shallow boreal atmosphere responds sharply to these radiative perturbations, and additional feedbacks from glaciers, sea ice, and sea surface temperatures amplify the responses.
Researchers have used models to explore the sensitivity of boreal climate to vegetation since the inception of credible land components within climate models. Otterman et al. (1984) calculated the high-latitude atmosphere to be 4.6 • C cooler without the masking of snow from vegetation. Thomas and Rowntree (1992) demonstrated that boreal deforestation could increase land albedo and snow depth, and decrease surface temperature, net radiation, heat fluxes, and precipitation during spring. Bonan et al. (1992) estimated this cooling to be upwards of 12 • C at 60 • N in April when amplified by ocean feedbacks. More recently, Betts (2000) calculated a +10 to +20 W m −2 surface shortwave forcing from afforestation in boreal regions. In a complemen- Our study builds upon this earlier work, and uses a data-informed modeling approach to characterize responses from a particularly influential and climate-sensitive land cover agent: boreal forest fires. Our surface forcings were in general agreement with recent modeling work in Alaska. Euskirchen et al. (2009) estimated an increase in summer albedo between 0.006 and 0.015, and winter albedo of 0.05, across the western arctic associated with burn area increases on the order of 100-200 % under future climates. For the same domain, we estimated albedo increases of 0.003 and 0.005 for summer, and 0.021 and 0.057 for winter, for BA × 2 and BA × 4, respectively. Randerson et al. (2006)  Recent studies of woody encroachment in the arctic have highlighted the importance of evapotranspiration (ET)driven atmospheric water vapor feedbacks (Swann et al., 2010;Bonfils et al., 2012). Swann et al. (2010) calculated this forcing to be 50 % larger than albedo feedbacks when bare ground in the arctic, which has zero transpiration, was replaced with deciduous broadleaf trees. Conversely, our study suggests that ET feedbacks are relatively unimportant for changing boreal forests: albedo-driven cooling from younger stands switched the sign of uncoupled ET anomalies and reduced atmospheric water content. This argues that the dominant biophysical driver of high-latitude vegetation feedbacks to climate is context-dependent. Researchers must consider the dynamics of all relevant energy fluxes when assessing terrestrial-atmosphere coupling in boreal and arctic biomes.

Summary of results and implications
The coupled model results presented here provide evidence for a near-linear response of boreal North American winter and spring surface temperatures to surface energy budget changes associated with a decrease or increase in burning up to about 100 %. With even greater burning, surface forcings were somewhat dampened, yet temperature responses were amplified. Driven by higher albedos, cooler temperatures in the increased burn scenarios reduced evaporative demand and led to a slightly drier atmosphere with a shallower boundary layer for most of the year. Consequently, although uncoupled simulations suggested greater summer evaporation with younger deciduous forests, specific humidity and precipitation displayed slight negative trends.
Although our modeling experiments were conducted at equilibrium, they may offer insight into past and future feedbacks to climate change in boreal North America because of the similar operating timescales for forest re-growth and climate changes (several decades). During the last 50 yr, this region warmed by about 0.3-0.4 • C per decade in summer , and nearly twice that in winter (Arctic Climate Impact Assessment, 2004). Based on the NFDBs, burn area has increased by approximately 0.39 Mha yr −1 per decade during this same time period. Our modeling results suggest that land surface changes due to fires may have cooled the surface by roughly 0.04 • C per decade in the winter, providing a negative feedback of approximately 5-7 %. Looking towards the future, recent syntheses project a 3-5 • C increase in annual surface air temperature, and 5-8 • C increase during the winter, over boreal North America by 2100 (Arctic Climate Impact Assessment, 2004;Christensen et al., 2007). Coincident with this warming, numerous research groups project increased burning in North American boreal forests. Based on our simulations, a doubling of burn area could cool the land surface by 0.23 • C during the winter and provide a negative feedback of 3-5 % across the domain, while a quadrupling could cool by 1.1 • C and counteract 14-23 % of the warming. For areas in central Alaska and Canada that experience frequent burning, we found cooling as high as 1.2 • C (25-34 % feedback) for a doubling of burn area, and 3 • C (38-60 % feedback) for a quadrupling, between February and April.
Our simulated climate impacts will not occur independently, but instead develop amidst an evolving 21st century landscape and climate. Ocean and sea ice feedbacks will respond non-linearly to an overall warming trend, and may not interact with this smaller cooling footprint to the same degree. A longer snow-free season will decrease the magnitude of spring albedo effects. Boreal and arctic vegetation will also likely respond to 21st century climate via mechanisms other than fire. With permafrost thawing and surface air warming, we have witnessed increases in the northern boundaries and productivities of boreal shrubs during the last several decades (Sturm et al., 2001;, leading to an estimated +0.2 W m −2 of atmospheric heating per decade since 1950 . It has been estimated that shrubland coverage under future climates has the potential to add over 6 W m −2 of atmospheric heating , with additional warming coming from atmospheric water vapor feedbacks via shrub transpiration . Boreal forests, in contrast, have become less productive due to drought stress from higher temperatures and longer growing seasons , and may take much longer than shrubs to migrate northward (Chapin III and Starfield, 1997;Rupp et al., 2001).
Boreal fires also affect climate in ways other than land surface characteristics. Pyrogenic aerosols released during combustion interact with the atmosphere and planetary radiation balance in complex ways. Typically, direct effects from fire aerosols produce warming in the local atmosphere but cooling at the surface (Stone et al., 2008;Tosca et al., 2012). In the Arctic, however, top of atmosphere energy balance depends on the albedo of the underlying surface (Ramanathan et al., 2001;Stone et al., 2008), and deposited black carbon can produce strong local warming and melting of snow and sea ice (Flanner et al., 2007;Ramanathan and Carmichael, 2008). When indirect effects on cloud formation are considered, the atmosphere is further cooled and radiative forcing becomes strongly negative (Ward et al., 2012).
Fires also release greenhouse gases: mostly CO 2 , but also CH 4 and N 2 O, which are well-mixed and produce positive global radiative forcings. A 50 % increase in circumpolar boreal forest fires has the potential to release 27-51 Pg C, or 0.33-0.8 Pg C yr −1 during the next 50-100 yr (Kasischke et al., 1995). Model estimates suggest that Canada's fire emissions may increase considerably more than this, driven mostly by increasing burn area (Amiro et al., 2009). The overall biogeochemical forcing, however, includes ecosystem recovery rates, which are different for different carbon pools and depend on climate, fire severity, and species composition.
Surface biophysical, aerosol, and biogeochemical effects from fires must all be considered when assessing the overall impacts of fires on climate. For one set of fires in Alaska, Randerson et al. (2006) found that surface albedo dominates the radiative budget, approximately doubling the other positive forcings in magnitude and leading to an overall negative radiative forcing. O'Halloran et al. (2012) documented a similar result using six sites in Manitoba, Canada, although the albedo and biogeochemical forcings were closer in magnitude. However, aerosol impacts were neglected in O' Halloran et al. (2012), and aerosol indirect effects were not included in Randerson et al. (2006). Neither study considered downwind land and ocean ecosystem responses to nutrient deposition or fire-induced changes in regional climate, including changes in the fraction of diffuse photosynthetically active radiation (Mahowald, 2011).

Caveats
While we attempted to quantify several influential yet accessible uncertainties, a variety of potential biases and errors remain untreated. One outstanding issue involves post-fire succession. Simulated vegetation distributions were direct functions of succession, which is highly uncertain. Our approach supported the hypothesis that forests in boreal North America primarily re-grow through relay floristics, wherein deciduous trees prevail during mid-succession. Although this has a large potential for error, comparisons of energy budgets to observations were favorable, and halving the prescribed deciduous fraction had a relatively minor impact on energy budgets (Table 5). More pertinent for uncertainty, however, is the assumption that these trajectories are constant across space, time, and levels of fire severity. Observations suggest post-fire succession is influenced by geography and topography, and depends on pre-fire conditions. Specifically, higher proportions of deciduous trees establish and survive during succession on south-facing slopes with improved drainage, and with greater pre-fire deciduous fractions, higher fire severity, and increased fire frequency. Deciduous forests are less flammable. Using a statistical model of fire occurrence and spread for a region in Alberta, Krawchuk and Cumming (2011) estimated a relative negative biotic feedback on burn area of 23 % under a scenario that otherwise predicted an increase of 61 %. Our simulations suggest a 22 % relative reduction with a 100 % increase in burn area. However, it has been suggested that an intensifying fire regime may not just increase deciduous fractions, but lead to an alternate stability point of deciduous forest dominance (Johnstone et al., 2010a). This would tend to amplify our modeled temperature effects, although a landscape dominated by deciduous trees may display entirely different fire dynamics.
Fire severity, or carbon combusted per unit burn area, may also change in the future and alter mortality and succession. While some studies suggest an increase in severity due to climate (e.g., Harden et al., 2000;Rogers et al., 2011) and drier forest floor organic matter (OM) (Bonan et al., 1990), less flammable deciduous forests with increasingly smaller pools of OM may decrease severity. Moreover, changes to boreal vegetation may alter the emissions of organic vapors that impact cloud condensation nuclei and regional cloud formation (Spracklen et al., 2008).
We did not directly characterize biases from our main driving datasets, including MODIS land cover and burn area from the NFDBs. There are known omissions in the NFDBs (Kasischke et al., 2002) and biases in high-latitude vegetation from MCD12Q1 (Friedl et al., 2010). However, our varying burn area prescriptions largely bypass minor errors in historical burn area, and vegetation was indirectly validated through comparisons of post-fire energy exchange. Our domain classification is also imperfect, undoubtedly including omission and commission errors in boreal forest 500 m pixels. Our scheme contains some forests not traditionally identified as boreal, such as those in British Columbia, and neglects forests in the continental US However, the forests of British Columbia tend to burn infrequently and contributed little to continental energy budget anomalies under our altered burn area scenarios. Forests of the continental US that may be climatically classified as boreal occur on the northern Rocky and Cascade Mountain Ranges, and would not be expected to contribute substantially to continental climate anomalies within our experimental design. We also did not attempt to address overarching model biases in either the CLM land model or the CESM climate model. Other Earth system models with different formulations for energy and water fluxes, atmospheric circulation, sea ice dynamics, ocean fluxes, and climate sensitivity may display different quantitative and qualitative responses. We also include no representation of timescale in our experiments, but instead simulate climate responses to a land surface at equilibrium.

Future directions
Further work that expands the spatial domain, fire forcing agents, and modeling techniques considered here is needed to better clarify the impacts of boreal forest fires on climate and ecosystems over the next century. Eurasia must be considered, as it contains approximately twice the coverage of boreal forests as does North America, with distinct and heterogeneous ecology and fire regimes. Although surface fluxes may dominate in high latitudes, ecosystem carbon budgets and fire emissions of aerosols and greenhouse gasses must be included in order to adequately represent the full climate impacts. With the continuing development of global land models, dynamic fire and vegetation modules should be used to capture the non-linear and interacting responses of succession, burn area, and fire severity. Additionally, studies that better characterize spatial and temporal relationships between succession, fire intensity, burn severity, and surface fluxes have the potential to reduce important sources of uncertainty.

Conclusions
Modeling results suggest that increases in fire activity across North American boreal forests, similar to those projected by various research groups for the 21st century, would result in a younger landscape that cools the high-latitude atmosphere. Amplified by ocean and sea ice feedbacks, impacts from fireinduced vegetation change were dominant in February-April when solar insolation was moderate and albedo changes were large. Annually, this would provide only a minor reduction to high-latitude climate warming across the domain with a doubling of historical burn area, but could decrease warming by nearly 25 % in the winter with a quadrupling of burn area. Although this result implies a negative feedback to climate change, future work is necessary to adequately quantify the spatial and seasonally varying effects from all fire forcing agents: landscape biophysical, aerosol emissions, and greenhouse gases. While landscape effects cool locally, the impacts of aerosol emissions are uncertain and pyrogenic greenhouse gases warm the global atmosphere.
This work implies that representing fire and vegetation succession in Earth system models is crucial in order to project high-latitude climate and ecosystem responses to anthropogenic forcings. Dynamic vegetation models that simulate distinct age, functional, and structural classes can help characterize contemporary succession in locations with limited data, and project how succession and migration will ulti-mately shape future distributions. These should be coupled to mechanistic fire models that respond to climate and represent spatially varying fuel dynamics to simulate fires of different types and intensities. Until such models are validated against historical and paleo-archives, and fully integrated into climate models, we can only provide boundary conditions for possible future impacts.
Finally, while this study suggests that boreal forest fires can provide regional negative feedbacks to warming, it must be noted that fires are associated with a number of harmful non-climate impacts. Fires decrease air quality, leading to deleterious health effects (Bowman and Johnston, 2005;Johnston et al., 2012), and damage homes and tourism industries throughout the world. While regular low-intensity fires may increase the health of forests in certain regions, large conflagrations negatively impact many ecosystem services, flora, and fauna.