Future threat to boreal ecosystem health from wildfire air pollution 1 2

Biomass burning is an important source of tropospheric ozone (O 3 ) and aerosols. These air pollutants can affect vegetation photosynthesis through stomatal uptake (for O 3 ) and light scattering (for aerosols). Climate change will significantly increase wildfire activity in boreal North America by the midcentury, while little is known about the impacts of enhanced emissions on the terrestrial carbon budget. Here, combining site-level and satellite observations and a carbon-chemistry-climate model, we estimate the impacts of fire emitted O 3 and aerosols on net primary productivity (NPP) over boreal North America. Fire emissions are calculated based on an ensemble projection from 13 climate models. In the present day, wildfire enhances surface O 3 by 2 ppbv (7 %) and aerosol optical depth (AOD) at 550 nm by 0.03 (26 %) in the summer. By midcentury, boreal area burned is predicted to increase by 66 %, contributing more O 3 (13 %) and aerosols (37 %). Fire O 3 causes negligible impacts on NPP because ambient O 3 concentration is far below the damage threshold of 40 ppbv. Fire aerosols reduce surface solar radiation but enhance atmospheric absorption, resulting in enhanced air stability and intensified regional drought. The domain of this drying is confined to the North in the present day, but extends southward by 2050 due to increased fire emissions. Consequently, wildfire aerosols enhance NPP by 72 Tg C yr -1 in the present day but decrease NPP by 118 Tg C yr -1 in the future, mainly because of the soil moisture perturbations. Our results suggest that future wildfire may accelerate boreal carbon loss, not only through direct emissions, but also through the biophysical impacts of fire aerosols.


Introduction
Wildfire is becoming more active in recent decades over North America boreal regions (Stocks et al., 2002;Kasischke and Turetsky, 2006).Fire activity is closely related to weather conditions and large-scale atmospheric oscillations (Gillett et al., 2004;Duffy et al., 2005), and is projected to increase significantly in the future due to climatic changes (Flannigan et al., 2005;Balshi et al., 2009).More frequent wildfires are accelerating carbon loss in boreal North America (Bond-Lamberty et al., 2007;Turetsky et al., 2011).Meanwhile, fire-induced air pollution, including ozone (O 3 ) and aerosols, is predicted to increase in boreal and downwind regions by midcentury (Yue et al., 2013;Yue et al., 2015).Wildfire emissions have large impacts on air quality (Wotawa and Trainer, 2000;Morris et al., 2006), weather/climate conditions (Randerson et al., 2006;Zhao et al., 2014), and public health (Zu et al., 2016;Liu et al., 2017).However, little is known about how these pollutants affect ecosystem carbon assimilation, and how this impact will change with the increased wildfire activity in the future.
Surface O 3 is detrimental to plant health because it damages photosynthesis through stomatal uptake (Sitch et al., 2007).In the present climate state, fire-induced O 3 enhancements are predicted to reduce net primary productivity (NPP) in the Amazon forest by 230 Tg C yr -1 , a magnitude comparable to the direct release of CO 2 from fires in South America (Pacifico et al., 2015).The aerosol effects are more uncertain because both positive and negative feedbacks occur.Appearance of aerosols increases diffuse light, which is beneficial for shaded leaves in the lower canopy.Consequently, photosynthesis of the whole ecosystem will increase as long as the total light availability is not compromised (Kanniah et al., 2012).Rap et al. (2015) estimated that biomass burning aerosols increase Amazon NPP by  Tg C yr -1 , which offsets about half of the damage caused by fire O 3 (Pacifico et al., 2015).In contrast, strong light attenuation associated with high aerosol loading may decrease canopy photosynthesis (Cohan et al., 2002;Oliveira et al., 2007;Cirino et al., 2014).Furthermore, the aerosol radiation changes indirectly influence land carbon uptake through concomitant meteorological perturbations that are only beginning to be examined (Yue et al., 2017).
Future wildfire activity is projected to increase over boreal North America but with large uncertainties (Flannigan et al., 2005;Tymstra et al., 2007;Girardin and Mudelsee, 2008;Nitschke and Innes, 2008;Amiro et al., 2009;Balshi et al., 2009;Bergeron et al., 2010; Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.Wotton et al., 2010;de Groot et al., 2013;Wang et al., 2016).For example, Amiro et al. (2009) predicted an increase of 34% in Canadian area burned for a 2×CO 2 scenario (2040-2060) relative to a 1×CO 2 condition , using the Canadian Fire Weather Index (CFWI) and output from Canadian global climate model (CGCM) version 1.On the other hand, Balshi et al. (2009) projected that area burned in boreal North America would double by the year 2045-2050 relative to 1991-2000, using the Multivariate Adaptive Regression Splines (MARS) approach and meteorological output from CGCM version 2. The increasing rate in Balshi et al. (2009) is higher than that in Amiro et al. (2009), indicating substantial uncertainties in fire projections originating from both fire models and simulated future climate.Recently, Yue et al. (2015) developed stepwise regressions over boreal ecoregions between area burned and multiple meteorological variables as well as CFWI indexes.They used output from 13 climate models to drive these regression models and predicted an average increase of 66% in boreal area burned at 2046-2065 relative to 1981-2000 under the IPCC A1B scenario (Solomon et al., 2007).The multi-model ensemble projection reduces uncertainties associated with single model climate projections.Yue et al. (2015) calculated that the wildfire emission increase by the 2050s would increase mean summertime surface O 3 by 5 ppbv in Alaska and 3 ppbv in Canada.The study found regional maximum O 3 enhancements as high as 15 ppbv, suggesting the potential for possible vegetation damage and land carbon loss due to the enhanced boreal fire-related air pollution.Wildfire aerosols are also expected to increase significantly but not predicted in Yue et al. (2015).
In this study, we quantify the impacts of O 3 and aerosols emitted from boreal wildfires on the land carbon uptake in North America in the present climate state and in the future world at 2050, taking advantage of the ensemble projection of future wildfire emissions by Yue et al. (2015).We first analyze relationships between gross primary production (GPP) and aerosol optical depth (AOD) at 550 nm over the boreal regions based on observations.We then perform a suite of Earth system model simulations using NASA GISS ModelE2 that embeds the Yale Interactive Terrestrial Biosphere model (YIBs), a framework known as ModelE2-YIBs (Yue and Unger, 2015).Future projections of wildfire emissions from Yue et al. (2015) are applied as input to ModelE2-YIBs model to project fire-induced O 3 and aerosol concentrations in the 2010s and 2050s.The impacts of the boreal fire O 3 on forest photosynthesis are predicted using the flux-based damage algorithm proposed by Sitch et al. (2007).Fire aerosols induce perturbations to radiation, meteorology, and hydrology, leading to multiple influences on the land carbon uptake.Sensitivity experiments are performed using

Observed GPP-AOD relationships
Following the approach by Strada et al. (2015), we investigate the GPP sensitivity to diffuse radiation and AOD variability in boreal regions.First, we identify study sites in Canada and Alaska from the AmeriFlux (AMF) network (http://ameriflux.lbl.gov/).The number of available boreal sites is much fewer than those in temperate regions.We select AMF sites providing hourly (or half-hourly) simultaneous measurements of GPP (non gap-filled) and photosynthetically active radiation (PAR, total and diffuse) for at least 3 consecutive years.
Only two Canadian sites meet the criteria: Groundhog River (CA-Gro, 82.2°W, 48.2°N), a mixed forest (MF), and Quebec Mature Boreal Forest Site (CA-Qfo, 73.4°W, 49.7°N), an evergreen needleleaf forest (ENF).At the two selected sites, we calculate the Pearson's correlation coefficients between half-hourly GPP and different components of PAR.We then apply instantaneous Level 2 Collection 6 of AOD pixels at 3-km resolution retrieved by the Moderate Resolution Imaging Spectroradiometer (MODIS, https://ladsweb.nascom.nasa.gov/)onboard the Aqua and Terra satellites (Levy et al., 2013).The MODIS 3-km AOD product has been fully validated against ground-based sun photometers at both global (Remer et al., 2013) and urban/suburban (Munchak et al., 2013)  productivity (Niyogi et al., 2004).For each summer day, we select instantaneous MODIS 3km AOD pixels that are (a) located within a distance of 0.03° (about 3 km) from the targeted AMF site and (b) "quasi-coincident" with AMF data, which are available each half-hour.
Because of the unavoidable temporal differences between MODIS overpass and AMF data availability, we name this selection "quasi-coincident".Cloud mask, applied in the MODIS retrieval procedure, conveniently filters out cloudy instants and should reduce the effect of clouds in the scattering process.We calculate both the correlation and regression coefficients between "quasi-coincident" GPP and AOD at the selected sites.Negative GPP is considered as a missing value.To further reduce the influence of cloud cover, we discard instants (both AMF and MODIS data) when precipitation is non-zero.The GPP-AOD sampling pairs are much fewer than GPP-PAR, because we select instants when both instantaneous AOD and GPP data are available.In addition, AOD is screened for clear instants to exclude the impacts of clouds.
Fuel consumption, the dry mass burned per fire area, is the product of fuel load and burning severity.For fuel load in Alaska, we use 1-km inventory from the US Forest Service (USFS)  Fuel Characteristic Classification System (FCCS, McKenzie et al., 2007).For fuel load in Canada, we use a 1-km fuel type map from the Canadian Fire Behavior Prediction (FBP) system (Nadeau et al., 2005), combined with fuel-bed definition from the FCCS.Burning severity, the fraction of fuel load burned by fires, is calculated with the USFS CONSUME model 3.0 following the approach described in Val Martin et al. (2012).With both fuel load and burning severity, we derive fuel consumption and further calculate biomass burned in boreal North America with the predicted area burned.Fire emissions for a specific species are then estimated as the product between biomass burned and the corresponding emission factor, which is adopted from measurements by Andreae and Merlet (2001) except for NO x .We use the average values from six studies as NO x emission factor, because the number reported in Andreae and Merlet (2001) is much higher than field observations (Alvarado et al., 2010).
Based on projected area burned and observation-based fuel consumption and emission factors, we derive fire emissions of NO x , carbon monoxide (CO), non-methane volatile organic compounds (NMVOCs, Alkenes and Alkanes), NH 3 , SO 2 , black (BC) and organic carbon (OC) in the present day and midcentury.

NASA ModelE2-YIBs model
The NASA ModelE2-YIBs is an interactive climate-carbon-chemistry model, which couples the chemistry-climate model NASA ModelE2 (Schmidt et al., 2014) and the YIBs vegetation model (Yue and Unger, 2015).NASA ModelE2 is a general circulation model with horizontal resolution of 2°×2.5°latitude by longitude and 40 vertical layers up to 0.1 hPa.It dynamically simulates both the physical (emissions, transport, and deposition) and chemical (production, conversion, and loss) processes of gas-phase chemistry (NO x , HO x , O x , CO, CH 4 , and NMVOCs), aerosols (sulfate, nitrate, ammonium, BC, OC, dust, and sea salt), and their interactions.In the model, oxidants influence the photochemical formation of secondary aerosol species (e.g., sulfate, nitrate, and biogenic secondary organic aerosol), in turn, aerosols alter photolysis rates and influence the online gas-phase chemistry.The model also considers interactions between climate and atmospheric components.(Wild et al., 2013), atmospheric composition (Shindell et al., 2013b), and radiative effects (Shindell et al., 2013a).
YIBs is a process-based vegetation model that dynamically simulates changes in leaf area index (LAI) through carbon assimilation, respiration, and allocation.Coupled photosynthesisstomatal conductance is simulated with the Farquhar-Ball-Berry scheme (Farquhar et al., 1980;Ball et al., 1987).Leaf-level photosynthesis is upscaled to canopy level by separating diffuse and direct light for sunlit and shaded leaves (Spitters, 1986).Plant respiration considers thermal dependence as well as acclimation to temperature (Atkin and Tjoelker, 2003).Soil respiration is calculated based on the carbon flows among 12 biogeochemical pools (Schaefer et al., 2008).Net carbon uptake is allocated among leaves, stems, and roots to support leaf development and plant growth (Cox, 2001).An interactive flux-based O 3 damage scheme proposed by Sitch et al. ( 2007) is applied to quantify the photosynthetic responses to ambient O 3 (Yue and Unger, 2014).The YIBs model has been benchmarked against in situ GPP from 145 eddy covariance flux tower sites and satellite retrievals of LAI and phenology (Yue and Unger, 2015).

Simulations
We perform 6 time-slice simulations, three for present-day (2010s) and three for midcentury (2050s), with atmosphere-only configuration to explore the impacts of fire emissions on NPP in boreal North America (Table 1).Simulations F10CTRL and F50CTRL turn off all fire emissions as well as O 3 vegetation damage for the 2010s and 2050s, respectively.However, climatic feedbacks of aerosols from other sources (both natural and anthropogenic) and related photosynthetic responses are included.Simulations F10AERO and F50AERO consider the responses of plant productivity to perturbations in radiation and meteorology caused by aerosols, including emissions from wildfires and other sources, but do not include any O 3 vegetation damage.In contrast, simulations F10O3 and F50O3 calculate offline O 3 damage based on the simulated O 3 from all sources including fire emissions.Vegetation NPP, the net carbon uptake by biosphere, is the impact metric.The difference between AERO and CTRL runs isolates the impacts of fire aerosols on NPP, and the difference between O3 and CTRL runs isolates O 3 vegetation damage caused by fire and non-fire emission sources.average well-mixed greenhouse gas concentrations and anthropogenic emissions of shortlived species, both at present day and midcentury, are adopted from the RCP8.5 scenario (Table 2).Natural emissions of soil and lightning NO x , biogenic volatile organic compounds (BVOC), dust, and sea salt are climate-sensitive and simulated interactively.The RCP8.5 land cover change dataset shows limited changes in land cover fractions between 2010s and 2050s (Oleson et al., 2010).For example, relative to the 2010s, a maximum gain of 5% is predicted for grassland in the 2050s, resulting from a 1% loss in deciduous forest and another 1% loss in needleleaf forest over boreal North America.As a result, a land cover dataset derived from satellite retrievals (Hansen et al., 2003) is applied for both the 2010s and 2050s.
To isolate the impact of individual aerosol-induced climatic perturbations on NPP, we perform 10 sensitivity experiments using the YIBs vegetation model driven with offline meteorology simulated by ModelE2-YIBs model (Table 3).For example, the offline run Y10_CTRL is driven with variables from the online simulation of F10CTRL (Table 1).The run Y10_TAS adopts the same forcing as Y10_CTRL except for temperature, which is simulated by the climate simulation of F10AERO.In this case, we quantify the NPP responses to individual and/or combined climate feedback (mainly in temperature, radiation, and soil moisture) by fire aerosols.Each offline run is conducted for 12 years and the last 10 years are used for analyses.

Observation datasets
We use observations to evaluated GPP, AOD, and O 3 in boreal North America simulated by ModelE2-YIBs.For GPP, we use a benchmark data product upscaled from FLUXNET eddy covariance data using an ensemble of regression trees (Jung et al., 2009).For AOD observations, we use satellite retrieval at 550 nm from Terra MODIS Level 3 data product.

Observed GPP-AOD relationships
Positive correlations between GPP and diffuse PAR are found at the two boreal sites (Figs 1b-1c).The magnitude of diffuse PAR is similar for these sites, possibly because they are located at similar latitudes (Fig. 1a).GPP values at CA-Gro are generally higher than that at CA-Qfo, likely because deciduous broadleaf forest (DBF) has higher photosynthetic rates.
Consequently, the slope of regression between GPP and PAR dif is higher at CA-Gro than that at CA-Qfo, suggesting that GPP of DBF (or MF) is more sensitive to changes in diffuse PAR than that of ENF.We find almost zero correlation between GPP and PAR dir at the two sites (Table 4), indicating that photosynthesis is in general light-saturated for sunlit leaves at these sites during boreal summer noontime.As a result, modest reductions in direct light by aerosols will not decrease GPP of the whole canopy.
With satellite-based AOD, we find positive correlations between GPP and AOD at both sites (Figs 1d-1e).However, the slope of regression between GPP and AOD is lower at CA-Gro compared with that at CA-Qfo, opposite to the GPP-PAR dif regressions.The cause of such discrepancy might be related to the limitation of data availability.For the same reason, the GPP-AOD correlation is insignificant at CA-Gro site.On average, we estimate GPP sensitivity of 3.5 ± 1.1 µmol m -2 s -1 per unit AOD at lower latitudes of boreal regions in the summer.

Evaluation of model photosynthesis and atmospheric composition
Simulated summer GPP shows high values in mid-western Canada (Alberta and Saskatchewan) and the Southeast (Ontario) (Fig. 2a).Forest GPP at high latitudes is low because of the cool weather and light limitation there.Compared against the data product, America and high values in the western and eastern U.S. (Fig. 3a).This pattern is consistent with surface observations (Fig. 3b), but the model overestimates the measured surface O 3 by 40%.The Canadian measurement sites are located near the southern boundary, and as a result do not represent the average state over the vast boreal region at higher latitudes.

Simulation of wildfire O 3 and aerosols
Compared to the present day, wildfire area burned is projected to increase by 66% in boreal North America at midcentury, mainly because of the higher temperature in future fire seasons (Yue et al., 2015).Consequently, enhanced fire emissions increase concentrations of surface O 3 and column AOD, especially over Alaska and central Canada (Fig. 4).The maximum centers of air pollutants are collocated for O 3 and AOD but with unproportional magnitudes, suggesting non-linear conversion among fire emission species as well as the interactions with natural emission sources (e.g., lightning/soil NO x and BVOC).On average, wildfire emissions contribute 7.1 ± 3.1% to surface O 3 and 25.7 ± 2.4% to AOD in the summer over boreal North America in the present day.By midcentury, these ratios increase significantly to 12.8 ± 2.8% for O 3 and 36.7 ± 2.0% for AOD.

Simulation of fire pollution impacts on NPP
Surface O 3 , including both fire and non-fire emissions (Table 2), causes limited damages to summer GPP in boreal North America (Fig. 5).The most significant damage is predicted over eastern U.S., where observed [O 3 ] is high over vast forest ecosystems (Fig. 3).In the western U.S., [O 3 ] is also high but the O 3 -induced GPP reduction is trivial because low stomatal conductance in the semi-arid ecosystems limits O 3 uptake there (Yue and Unger, 2014).Over boreal regions, the mean [O 3 ] of 28 ppbv is much lower than the damaging threshold of 40 ppbv for most tree species (Yue et al., 2016).Statistics in Yue et al. (2015) show that maximum daily 8-hour average (MDA8) [O 3 ] can be higher than 40 ppbv in Alaska Fire aerosols cause significant perturbations in shortwave radiation at surface (Fig. 6).The direct light is largely attenuated especially over Alaska and central Canada, where fire aerosols are most abundant (Fig. 4).In contrast, diffuse light widely increases due to particle scattering.In the present day, the average reduction of 5.6 W m -2 in the direct light component is in part offset by the enhancement of 2.6 W m -2 in the diffuse light component, leading to a net reduction of 3.0 W m -2 in solar radiation over boreal North America.By the midcentury, a stronger reduction of 9.5 W m -2 in direct light is accompanied by an increase of 4.0 W m -2 in diffuse light, resulting in a net reduction of 5.5 W m -2 in solar radiation.Fireinduced BC aerosols strongly absorb solar radiation in the atmospheric column (Figs 7a-7b).
On average, fire aerosols absorb 1.5 W m -2 in the present day and 2.6 W m -2 by the midcentury.7c-7d).Surface radiative cooling and atmospheric heating together increase air stability and induce anomalous subsidence.In the present day, such descending motion is confined to 55-68°N, accompanied by a rising motion at 52-55°N (Fig. 7c).As a result, fire aerosols induce surface warming at higher latitudes but cooling at lower latitudes in boreal regions (Fig. 8a).

Atmospheric circulation patterns respond to the aerosol-induced radiative perturbations (Figs
Meanwhile, precipitation is inhibited by the subsidence in northwestern Canada but is promoted by the rising motion in the Southwest (Fig. 8c).By the midcentury, the range of subsidence expands southward to 42°N (Fig. 7d) due to strengthened atmospheric heating (Fig. 7b).The downward convection of warm air offsets surface radiative cooling (Fig. 6b), leading to a significant warming in the Southwest (Fig. 8b).The expanded subsidence further inhibits precipitation in vast domain of Canada (Fig. 8d).Soil moisture is closely related to rainfall and as a result exhibits dipole changes (drier north and wetter south) in the present day (Fig. 8e) but widespread reductions (Fig. 8f) by the midcentury.
In response to the climatic effects of fire aerosols, boreal NPP shows distinct changes between the present day and midcentury (Fig. 9).In the 2010s, forest NPP increases by 5-15% in Alaska and southern Canada, but decreases by 5-10% in northern and eastern Canada.This pattern of NPP changes (ΔNPP) is connected to the climatic effects of aerosols, especially changes in soil moisture (Fig. 8).The correlation between ΔNPP (Fig. 9a) and changes in soil moisture (Fig. 8e) reaches R = 0.56 (n = 356), much higher than the values of R = -0.11for temperature change (Fig. 8a) and R = 0.22 for precipitation change (Fig. 8c).At the continental scale, the patchy responses of NPP offset each other.Since the dominant fraction of carbon uptake occurs in southern Canada (Fig. 2a), where positive NPP change is predicted (Fig. 9a), wildfire aerosols enhance the total NPP by 72 Tg C yr -1 in the present day (Table 5).In contrast, increased wildfire emissions in the 2050s inhibit precipitation and decrease soil moisture in boreal North America (Fig. 8), leading to widespread NPP reductions and a total NPP loss of 118 Tg C yr -1 (Fig. 9b, Table 5).

Roles of aerosol climatic feedback
The contrasting sign of NPP responses in the present day and midcentury are closely related to the aerosol-induced surface climatic feedback.Sensitivity experiments using offline YIBs model (Table 3) allowed assessment of the impacts of individual changes in the major meteorological drivers, including temperature, radiation (diffuse and direct), and soil moisture (Table 5).The offline simulations driven with changes in all three variables yield ΔNPP of 126 Tg C yr -1 for the 2010s and -97 Tg C yr -1 for the 2050s.These values are different from the online simulations, which predict ΔNPP of 72 Tg C yr -1 for the 2010s and -118 Tg C yr -1 for the 2050s.Missing of other aerosol climatic feedbacks in the offline model, for example, changes in relative humidity, surface pressure, soil temperature, and turbulence momentum, may cause such discrepancy between the online and offline simulations.
Seasonal analyses show that summertime ΔNPP is 99 Tg C at present day and -95 Tg C at midcentury, dominating the NPP changes all through the year, because both wildfire emissions and ecosystem photosynthesis maximize in boreal summer.
Observations show that aerosols can promote plant photosynthesis through increasing diffuse radiation (Niyogi et al., 2004;Cirino et al., 2014;Strada et al., 2015).Our analyses with ground and satellite data also show positive correlations between GPP and AOD (Fig. 1).
The estimated NPP changes of 8 Tg C yr -1 by the radiative effects of boreal fire aerosols are much weaker than the enhancement of 78-156 Tg C yr -1 by fires in Amazon basin (Rap et al., 2015).There are at least two reasons for such discrepancies in the DFE between boreal and Amazon fire aerosols.First, wildfire emissions and associated impacts on radiation are much smaller in boreal regions.Wildfires in Alaska and Canada directly emit 68 Tg C yr -1 at the 2010s, resulting in enhancement of summer AOD by 35% and diffuse radiation by 1.7%.These boreal emissions are much smaller than the ~240 Tg C yr -1 in Amazon basin (van der Werf et al., 2010), where fires enhances regional PM2.5 concentrations by 85% and diffuse radiation by 6.2% in dry seasons (Rap et al., 2015).Second, larger solar insolation in lower latitudes allows stronger DFE for the same unit change of diffuse radiation.In our prediction, most of NPP changes occur at high latitudes of boreal regions (Fig. 9), where total insolation is not so abundant as that at the tropical areas.Consequently, decline of direct radiation in boreal regions more likely converts the light availability of sunlit leaves from light-saturation to light-limitation, offsetting the benefit from enhanced diffuse radiation for shaded leaves.
For this study, we do not find GPP reduction by the decline of direct light at the two Ameriflux sites (Table 4), possibly because these sites are located at middle latitudes (<50°N).
In the future, more observations at higher latitudes (> 55°N) are required to explore the sensitivity of GPP to AOD at the light-limited conditions.
Simulations have shown that absorbing aerosols can cause regional drought by increasing air stability (Liu, 2005;Cook et al., 2009;Tosca et al., 2010).Our results confirm such tendency but with varied range of hydrological responses depending on the magnitude of wildfire  8f).Observations suggest that precipitation (and the associated soil moisture) is the dominant driver of the changes in GPP over North America, especially for the domain of cropland (Beer et al., 2010).Sensitivity experiments with offline YIBs model show that changes in soil moisture account for 82.5% of ΔNPP at present day and 70.5% of ΔNPP at midcentury (Table 5).These results suggest that aerosol-induced changes in soil water availability, instead of temperature and radiation, dominantly contribute to the changes of boreal NPP, consistent with observational and experimental results (Ma et al., 2012;Girardin et al., 2016;Chen et al., 2017).

Limitations and uncertainties
For this study, we use the ensemble projected fire emissions from Yue et al. (2015).Area burned is predicted based on the simulated meteorology from multiple climate models.Such an approach may help reduce model uncertainties in climatic responses to CO 2 changes, but cannot remove the possible biases in the selection of climate scenarios and fire models.All predictions in Yue et al. (2015) are performed under the IPCC A1B scenario.With two different scenarios, A2 of high emissions and B2 of low emissions, Balshi et al. (2009) showed that future area burned in boreal North America increases at a similar rate until the 2050s, after which area burned in A2 scenario increases much faster than that in B2 scenario.
On average, boreal area burned in Balshi et al. (2009) increases by ~160% at 2051-2060 compared with 2001-2010, much higher than the change of 66% in Yue et al. (2015).In contrast, Amiro et al. (2009) predicted that boreal area burned at the 2×CO 2 scenario increases only by 34% relative to the 1×CO 2 scenario.This ratio is only half of the estimate in Yue et al. (2015), which compared results between periods with 1.44×CO 2 and 1×CO 2 .
The discrepancies among these studies are more likely attributed to the differences in fire models.Although both Amiro et al. ( 2009) and Yue et al. (2015) developed fire-weather regressions in boreal ecoregions, the former study did not include geopotential height at 500 hPa and surface relative humidity as predictors, which make dominant contributions to area burned changes in the latter study.On the other hand, Balshi et al. (2009) developed nonlinear regressions between area burned and climate at grid scale, which helps retain extreme values at both the temporal and spatial domain.Compared to previous estimates, Yue et al. (2015) predicted median increases in future fire emissions over boreal North America.Predicted surface [O 3 ] is much higher than observations over boreal North America (Fig. 3).This bias does not affect main conclusions of this study, because predicted O 3 causes limited damages to boreal GPP even with the overestimated [O 3 ] (Fig. 5).The result confirms that fire-induced O 3 vegetation damage is negligible in boreal North America.For aerosols, the model captures reasonable spatial pattern of AOD but with a background value of ~0.1 outside fire-prone regions, where the actual AOD is usually 0.1-0.2 (Fig. 2).This discrepancy may be related to the insufficient representations of physical and chemical processes in the model, but may also result from the retrieval biases in MODIS data due to the poor surface conditions (Liu et al., 2005) and small AOD variations (Vachon et al., 2004) at high latitudes.
Simulated aerosol climatic effects depend on radiative and physical processes implemented in the climate model.We find that present-day boreal fire aerosols on average absorb 1.5 W m -2 in the atmosphere (Fig. 6), which is much smaller than the value of 20.5 ± 9.3 W m -2 for fires in equatorial Asia (Tosca et al., 2010).This is because boreal fires enhance AOD only by 0.03 while tropical fires increase AOD by ~0.4.Previous modeling studies showed that fire plumes induce regional and downwind drought through enhanced atmospheric stability (Feingold et al., 2005;Tosca et al., 2010;Liu et al., 2014).Most of these results were based on the direct and/or semi-direct radiative effects of fire aerosols.Inclusion of the indirect aerosol effect may further inhibit precipitation and amplify drought, but may also introduce additional uncertainties for the simulations.The fire-drought interaction may promote fire activity, especially in a warmer climate.Ignoring this interaction may underestimate future area burned and the consequent emissions.

Implications
Inverse modeling studies have shown that the land ecosystems of boreal North America are carbon neutral in the present day, with the estimated land-to-air carbon flux from -270 ± 130 Tg C yr -1 to 300 ± 500 Tg C yr -1 (Gurney et al., 2002;Rodenbeck et al., 2003;Baker et al., 2006;Jacobson et al., 2007;Deng et al., 2014).Here, we reveal a missing land carbon source due to future wildfire pollution, taking into account full coupling among fire activity, climate change, air pollution, and the carbon cycle.Fire pollution aerosol increases boreal NPP by 72 Tg C yr -1 in the present day, comparable to the direct carbon loss of 68 Tg C yr -1 from wildfire CO 2 emissions.By midcentury, increasing fire emissions instead cause a NPP reduction of 118 Tg C yr -1 due to the amplified drought.Although NPP is not a direct   forest (MF at Gro) and evergreen needleleaf forest (ENF at Qfo) (Table 1).Data cover summer days (June-August).AmeriFlux diffuse PAR and GPP (in µmol m −2 s −1 ) are halfhourly observations (10:00-14:00 LT).Instantaneous MODIS Aqua and Terra 3-km AOD are selected in a time span centered on AmeriFlux record time.For each plot: the red line indicates the regression line, black lines depict the 1-σ interval; the regression slope and correlation coefficient are both included for each site (in bold if statistically significant at 95% confidence level).Blue dots in (b, c) show instants when MODIS Aqua and Terra 3-km AODs overlap AmeriFlux data.(d, e) GPP and MODIS AOD at (a) two boreal sites: Groundhog River (Gro) and Quebec Mature Boreal Forest Site (Qfo).The two sites are from the AmeriFlux network in Canada and are dominated by mixed forest (MF at Gro) and evergreen needleleaf forest (ENF at Qfo) (Table 1).Data cover summer days (June-August).AmeriFlux di use PAR and GPP (in µmol m 2 s 1 ) are half-hourly observations (10:00-14:00 LT).Instantaneous MODIS Aqua and Terra 3-km AOD are selected in a time span (30-minutes long) centered on AmeriFlux record time.For each plot: the red line indicates the regression line, black lines depict the the 1-interval; the regression slope and correlation coe⇥cient are both included for each site (in bold if statistically significant at 95% confidence interval).Blue dots in (b, c) show instants when MODIS Aqua and Terra 3-km AODs overlap AmeriFlux data.observations.Simulation results are from F10AERO (Table 1).Each point on the (e, f) scatter 915 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License. the YIBs model in offline mode to isolate the contributions of changes in the individual meteorological drivers.
scales.Strada et al. (2015) used groundbased AOD observations from the Aerosol Robotic Network (AERONET) near AMF sites to validate the sampling technique of MODIS 3-km AOD product.For this study, the validation against ground-based AOD observations was not possible because no AERONET stations exist near to the selected AMF sites.Every day, MODIS satellite sensors pass a specific region between 10:00 and 14:00 Local Time (LT), leaving patchy signals around the AmeriFlux sites.Most of MODIS AOD data at high latitudes are available only in boreal summer; as a result, we narrow our explorations of the GPP-AOD relationships to the noontime (10:00-14:00 LT) from June to August.The chosen noontime window limits the contributions that confounding factors such as low solar angles and high diffuse fraction may have on the amount of diffuse PAR and plant Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.
Simulated climate affects formation, transport, and deposition of atmospheric components, in turn, both O 3 and aerosols influence climate by altering radiation, temperature, precipitation, and other climatic variables.Both observation-based evaluations and multi-model inter-comparisons indicate that ModelE2 demonstrates skill in simulating climatology (Schmidt et al., 2014), radiation Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.All simulations are conducted for 20 years and outputs for the last 15 years are used for analyses.The simulations apply sea surface temperatures (SSTs) and sea ice distributions from previous NASA GISS experiments under the IPCC RCP8.5 scenario (van Vuuren et al., 2011) archived for the Coupled Model Intercomparison Project Phase 5 (CMIP5).Decadal average monthly-varying SST and sea ice of 2006-2015 are used as boundary conditions for present-day (2010s) runs while that of 2046-2055 are used for future (2050s) runs.Decadal For O 3 , gridded datasets are not available.We use site-level observations from 81 U.S. sites at the Clean Air Status and Trends Network (CASTNET, https://www.epa.gov/castnet) and Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.202 Canadian sites at the National Air Pollution Surveillance (NAPS, http://www.ec.gc.ca/rnspa-naps/) program.All datasets are averaged over the 2008-2012 period to represent present-day climatological conditions.Gridded datasets are interpolated to the same 2°×2.5°resolution as ModelE2-YIBs model.
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.simulated GPP reasonably captures the spatial distribution with a high correlation coefficient of 0.77 (p << 0.01) and relatively small biases within 20%.Simulated AOD reproduces the observed spatial pattern including the high values in boreal forests (Fig. 2b).In contrast to the MODIS observations, predicted AOD is relatively uniform over the West with a background value of ~0.1.This discrepancy explains the low correlation coefficient (R = 0.25, p<0.01) between the model and MODIS data.Simulated [O 3 ] shows low values in boreal North Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.and Canada.However, such episodes appear at 95 percentile for present day and 90 percentile for midcentury, suggesting that O 3 vegetation damage is rare in boreal North America and fire-induced O 3 enhancement does not exacerbate such damages.Therefore, we do not consider O 3 damage effects further.
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.indicator of the land carbon sink, reduction of NPP is always accompanied with the decline of net ecosystem exchange (NEE) and the enhanced carbon loss.In combination with the enhanced carbon emission of 130 Tg C yr -1 , future boreal wildfire presents an increasing threat to the regional carbon balance and global warming mitigation.Furthermore, the NPP reductions are mostly located in southern Canada, where cropland is the dominant ecosystem, newly exposing the future wildfire-related air pollution risk to food production.SST, [CO 2 ], and emissions are adopted from RCP8.5 scenario, with the average of 2006-2015 for the 2010s and that of 2046-2055 for the 2050s.For fire emissions, values at the 2010s are predicted based on meteorology for 1981-2000 and those at the 2050s are for 2046-2065.Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2017-319Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 23 June 2017 c Author(s) 2017.CC BY 3.0 License.846 847 Figure 1.Relationships between (b, c) GPP and diffuse PAR and (d, e) GPP and MODIS AOD at (a) two boreal sites: Groundhog River (Gro) and Quebec Mature Boreal Forest Site (Qfo).The two sites are from the AmeriFlux network in Canada and are dominated by mixed

Figure 1 :
Figure1: Relationships between (b, c) GPP and di use PAR and (d, e) GPP and MODIS AOD at (a) two boreal sites: Groundhog River (Gro) and Quebec Mature Boreal Forest Site (Qfo).The two sites are from the AmeriFlux network in Canada and are dominated by mixed forest (MF at Gro) and evergreen needleleaf forest (ENF at Qfo) (Table1).Data cover summer days (June-August).AmeriFlux di use PAR and GPP (in µmol m 2 s 1 ) are half-hourly observations (10:00-14:00 LT).Instantaneous MODIS Aqua and Terra 3-km AOD are selected in a time span (30-minutes long) centered on AmeriFlux record time.For each plot: the red line indicates the regression line, black lines depict the the 1-interval; the regression slope and correlation coe⇥cient are both included for each site (in bold if statistically significant at 95% confidence interval).Blue dots in (b, c) show instants when MODIS Aqua and Terra 3-km AODs overlap AmeriFlux data.

1Figure 2 .
Figure 2. Evaluation of simulated summer (a) GPP and (b) AOD at 550 nm with (c, d) 914 Figure 3. Evaluation of simulated summer surface [O 3 ] with observations for 2008-2012.Observations are collected from 81 U.S. sites at the Clean Air Status and Trends Network (CASTNET) and 202 Canadian sites at the National Air Pollution Surveillance (NAPS) program.The number of points (n), correlation coefficient (r), and mean bias (b) for the evaluation are presented on the plot.Values over Canada and Alaska are denoted with blue points.

Figure 4 .Figure 5 .(
Figure 4. Changes in summer (a, b) [O 3 ] and (c, d) AOD at 550 nm induced by wildfire 967 emissions in (a, c) the 2010s and (b, d) the 2050s over boreal North America.Only 968 significant changes (p<0.05) are shown.969 970 971 Figure 8. Predicted changes in summertime (a, b) surface air temperature, (c, d) precipitation, and (e, f) soil water content at surface caused by aerosols from wildfire emissions at (a, c, e) present day and (b, d, f) midcentury.Results for temperature and precipitation are shown as absolute changes.Results for soil water are shown as relative changes.Results for the 2010s are calculated as (F10AERO -F10CTRL).Results for the 2050s are calculated as (F50AERO -F50CTRL).Significant changes (p<0.05) are marked with black dots.