Wildfire plumes in the Western US are reaching greater heights and injecting more aerosols aloft as wildfire activity intensifies

By producing a first-of-its-kind, decadal-scale wildfire plume rise climatology in the Western U.S. and Canada, we identify trends toward enhanced plume top heights, aerosol loading aloft, and near-surface smoke injection throughout the American West. Positive and significant plume trends suggest a growing impact of Western US wildfires on air quality at the local to continental scales and support the notion that wildfires may have an increasing impact on regional climate. Overlap of identified trends with regions of increasing wildfire emissions and burn severity suggests a link to climate driven trends toward enhanced wildfire activity. Further, time series of plume activity point to a possible acceleration of trends over recent years, such that the future impacts to air quality and regional climate may exceed those suggested by a linear fit to the multi-decadal data. These findings have significant implications for human health and exacerbate concern for the climate–wildfire connection.


Results
We use the Freitas et al. 47,53 plume rise model (hereafter abbreviated as "F2010", see "Freitas et al. plume rise model" section in "Methods" section), and an entrainment-mass conservation approach (see "Entrainment-mass conservation" section in "Methods" section) to estimate the vertical mass distribution of detrained smoke. Thus, we construct a plume rise climatology of ~ 4.6 million plumes occurring within the Western US and Western Canada during the months of August and September for the years 2003-2020. Model outputs are combined with PM 2.5 emissions from the Quick Fire Emissions Dataset (QFED) 61 and are considered in the context of PBL heights adopted from the Stochastic Time-Inverted Lagrangian Transport model (STILT; see "PBL estimates from STILT" section in "Methods" section) 62,63 . Linear regression and hinge fitting analyses (two-piece linear regressions, see "Trend analyses" section in "Methods" section) are used for insight into the trend of wildfire plume top heights and the vertical distribution of wildfire injected aerosol mass at the ecoregion scale. Ecoregions featuring prominently within the results are labeled in Fig. 1, with their corresponding number of simulated plumes per year presented in Table S1. Model evaluation is performed via comparison to the MISR plume heights dataset 64 . Emphasis is placed on August and September plumes given previously-identified wildfire emissions trends 17 and the computational/storage expense of generating regional-decadal scale model inputs.
Results of trend analyses should be considered within the context of potentially large uncertainties. Model inputs characterizing the wildfire heat flux and active burn area are influenced by satellite spatial resolution, the ability of clouds to obscure satellite observations, the assumption of a constant heat release per unit biomass burned, and the imperfect synthesis of burned area and heat release datasets ("Wildfire inputs" section). Further, though the model framework captures the core of the MISR observed plume top height distribution, there is evidence of a systemic underestimation of the frequency of the tallest plume tops, likely impacting the presented trend analyses.

Trends in plume top height and aerosol mass injection. Ecoregion scale analyses within the Western
US and Canada indicate that wildfire plume top heights are increasing throughout much of the mountainous Western US between the years of 2003 and 2020. In terms of the maximum plume top height during the months of August and September, we find trends as large as ~ 230 m per year in the Sierra Nevada ecoregion and trends in excess of 100 m per year in the Southern Rockies, Eastern Cascades, Arizona-New Mexico Mountains, and the Wasatch-Uinta ecoregions (Fig. 1a). If performed for the 5, 10, 30, or 100 largest plume top heights per year in each ecoregion, robust trends toward enhanced plume top heights persist in the Sierra Nevada and Eastern Cascades. Trend analyses considering the PM 2.5 flux (kg PM 2.5 emitted within the plume per second)-weighted average plume top height (Fig. 1b) further support the growth of plume top heights for the strongest emitting fires across much of the Western US, including notable trends in the mountain adjacent Colorado Plateau and Mojave Basin and Range ecoregions (> 30 m/year). While statistically insignificant, hinge fits (see "Trend analyses" section in "Methods" section) to the data suggest an acceleration of flux-weighted average plume top trends from ~ 2018 to 2020, highlighting the unprecedented nature of recent Western US wildfire activity relative to 2003-mid-2010s. In the Klamath Mountains, Sierra Nevada, and Southern Rockies, we find hinge fit trends (2017/2018-2020) of 99.0, 201.2, and 256.1 m/year, respectively, representing enhancements by a factor of 9-15 relative to results from linear regression (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). The presence of positive, significant (p < 0.1), and potentially accelerating trends throughout much of the mountainous Pacific Northwest and California, as well as the Southern Rockies, overlaps with previously identified potential wildfire emissions "hotspots" 17 .
We further undertake an examination of trends in the most extreme fire activity, as represented by pyrocumulonimbus (pyroCb) formation, defined by especially tall plumes with a plume top temperature < − 38 °C (the homogenous freezing temperature for liquid water) being leveraged as a pyroCb proxy 65 . We find the possibility of increasingly frequent pyroCb activity within the Colorado Plateau (~ 0.12 pyroCb hours/year). Although trends for other ecoregions and the Western US/Canada as a whole are statistically insignificant, this climatology finds the first occurrence of pyroCb development in 6 ecoregions between 2017 and 2020, demonstrating an acceleration in recent years. These ecoregions include the Colorado Plateau, Southern Rockies, Sierra Nevada, Wasatch-Uinta, Mojave Basin and Range, and Northwestern Great Plains. Hinge fits further support the possibility of a recent uptick in pyroCb activity, with hinged trends (p > 0.1; 2018-2020) of 2.72, 2.15, and 0.43 pyroCb hours/year for the Colorado Plateau, Southern Rockies, and Sierra Nevada, respectively. A comparison of simulated pyroCb activity for 2013 to the observational inventory of pyroCb activity constructed by Peterson et al. 65 indicates that this model framework reasonably captured the pyroCb formation above the Pony/Elk and Beaver Creek fires in Idaho, while failing to reflect the observed pyroCb activity above the Rim fire in California. Should  www.nature.com/scientificreports/ pyroCb trends present greater clarity over the coming years, this result would support the growing concerns surrounding pyroCb induced warming of the stratosphere 38,39 . An examination of aerosol mass injection indicates trends toward enhanced injection aloft throughout much of the Western US. Trends in aerosol mass injection above or within the PBL (see "PBL estimates from STILT" section in "Methods" section) largely follow total emissions trends (Fig. 2), highlighting wildfire activity in the Klamath Mountains and the Sierra Nevada. This result is noteworthy given the evidence for climate driven enhancements, via atmospheric aridity, to wildfire activity in these regions [18][19][20][21] . When limiting timeseries to 2003-2019 (a rather quiet fire year in the Western US), trend analyses further support the presence of positive and significant trends across numerous ecoregions (p < 0.1; p < 0.05-above the PBL in the Sierra Nevada and for all trends in the Klamath Mountains), bolstering confidence in this result given the effect of extreme 2020 values (Fig. 3) on trends fitted to 2003-2020. Altitude-based thresholds highlight significant injection trends up to 7 km altitude for the Sierra Nevada, Southern Rockies, Colorado Plateau, and Arizona-New Mexico mountains (Fig. 4i). A diurnal breakdown of the aforementioned PBL injection trends, based on a 3-h rolling window, indicates that trends are most robust near 18:00 local time (Fig. S1), a result that corresponds to the late afternoon maximum in wildfire activity 54,55 and the expectation of midday PBL growth. The relevance of enhanced injection above the PBL and at altitudes from 1 to 8 km (Fig. 4) to climate impacts is drawn from recognition of the importance of altitude to extending the aerosol atmospheric lifetime 29,31,32 , enhancing the direct radiative effects of an absorbing aerosol layer 29,68 , and expanding the potential for impact on cloud processes. Recent aircraft-based observations indicate the ability of smoke aerosols to suppress downwind precipitation by acting as a source of cloud condensation nuclei, lowering the average cloud droplet size and staving off initiation of collision-coalescence 12 . Further, enhanced aerosol injection aloft increases the potential for regional-continental scale smoke transport. Thus, Western US and Canadian wildfires will have increasing potential to deteriorate air quality and visibility across much of North America 69 . Concurrent trends toward enhanced aerosol injection within the PBL (Fig. 2b) contribute to elevated risks of wildfire deteriorated air quality at the local scale, a feature of these results that is particularly concerning given a growing Western US population 70 . Hinge fit analyses of aerosol mass injections amplify concern for the potential climate and air quality impacts of Western US wildfires, as results depict the possibility of an emergent acceleration of trends relative to linear trends spanning 2003-2020 ( Fig. 3, Fig. S2). While currently statistically insignificant, hinge fit results point toward the possibility of future Western US wildfire impacts to air quality and climate that are in excess of those suggested by a linear fit to 2003-2020. This possible acceleration of trends is in accordance with the extraordinarily large 2020 fire season 71 , and we strongly suspect analysis of the 2021 fire season would lend additional support.
Investigation into the evolution of above PBL emissions as a fraction of total wildfire injected aerosol reveals significant trends toward a growing percentage of emitted mass being injected above the PBL across the Sierra Nevada, Southern Rockies, and Northern Rockies ecoregions. We hypothesize that climate driven trends toward enhanced wildfire activity are responsible for increasing high altitude smoke injections, particularly across the Sierra Nevada and Southern Rockies. Warming enabled trends toward increased high altitude burning throughout the Sierra Nevada and portions of the Rockies 18 suggest a possible explanation for trends toward preferential smoke injection above the PBL.

Evaluation against MISR observed plume top heights.
Despite the uncertainties associated with generating a regional-decadal scale plume rise climatology (~ 4.6 million plumes), and those inherent to the current state of plume rise modeling 44 , we find that the modeled plumes evaluate well against MISR observed plume top heights when considering the context of previous evaluation efforts 45,56,72 . In an aggregate sense, we see that the distribution of model-derived plume top heights occurring during the late morning MISR overpass window captures much of the distribution of MISR observed plumes (Fig. 5a). While there is a systematic under-estimation of the proportion of high-altitude plumes, this result is in line with prior analyses employing the F2010 model 45 and is a current reality of the challenges facing physics-based plume top estimates. We find that this model-derived climatology captures the range of MISR observed plume tops, though the numerous high-altitude plumes are proportionally less frequent than what is seen in the observations. As a caveat to this underestimation, it should be noted that the MISR datasets preferentially samples and digitizes larger plumes 51 .
On a plume-by-plume basis, we find a standard deviation of modeled plume top heights from MISR observations of 46.9 m for 425 sets of reasonably collocated plumes (within 30 min and 4 km). Of these 425 sets of collocated plumes, we see that 47.8% of modeled plume tops fall within the ± 500 m observational uncertainty of MISR and that there is little geographic clustering of model uncertainty (Fig. 5b), therefore suggesting consistency of model results across ecoregions. Prior analyses considering smaller spatial and/or temporal domains have reported up to 64% of modeled plume top heights falling within the MISR observational uncertainty 56,72 .
With relevance to above/within PBL injections, 41.4% (176 plumes) of the 425 sets of collocated plumes are found to straddle the PBL, with the modeled plume top above the PBL and the observed plume top within the PBL, or vice versa. In 71% of such cases, the modeled plume top exceeded the PBL, with PBL heights averaging 77.9% of the modeled plume top heights. In terms of the climatological average of the modeled vertical distribution of smoke below plume top, ~ 78% of the PM 2.5 associated with these modeled plumes is being injected within the PBL. Given that the MISR observed plume top for these cases was contained within the PBL, this equates to an average model underestimate of the PM 2.5 injected within the PBL of ~ 22% for such cases. Similarly, if it is assumed that the smoke is distributed evenly throughout a MISR observed plume and that the retrieved plume top heights should follow a normal distribution for any individual plume, an average of ~ 77% of the injected PM 2.5 is contained within the PBL for the 29% of cases in which the MISR observation exceeded the PBL and the modeled plume top did not. Based on the findings of Val Martin et al. 51 , this assumption of a normal distribution for MISR retrieved heights is reasonable in many cases, though a lognormal distribution is likely more www.nature.com/scientificreports/ www.nature.com/scientificreports/ common. As such, it is reasonable to expect that > 77% of the PM 2.5 associated with these MISR observed plumes is retained within the PBL, enhancing the agreement with modeled injections. Overall, this analysis supports the idea of largely cancelling over/underestimation of PM 2.5 injection relative to the PBL, though given the greater frequency of modeled plume tops exceeding the PBL when MISR observations do not, it appears likely that the net effect is a modest overestimation of above PBL PM 2.5 injection.

Discussion
We find that trends toward enhanced wildfire activity correspond to elevated wildfire plume top heights and aerosol injection aloft for the majority of mountainous ecoregions across the Western US. August-September plume top maxima (Fig. 1a) and PM 2.5 flux-weighted average plume tops (Fig. 1b) exhibit statistically significant increasing trends across many Western U.S. ecoregions. These trends were collocated with increases in wildfire emissions (Fig. 2c), most notably across the Sierra Nevada ecoregion. Recast in terms of the aerosol mass injected, we find that the spatial variability of injection trends above and within the PBL (Fig. 2a,b) is largely explained by variability in wildfire emissions trends, albeit the Sierra Nevada and portions of the Rockies demonstrate a shift towards preferential smoke injections above the PBL. We hypothesize that the Sierra Nevada and Southern Rockies are exhibiting trends toward a greater fraction of smoke being emitted above the PBL due to climate enabled enhancements in high altitude burning 18,73 . Furthermore, there is a potentially emergent acceleration of flux-weighted average plume top, aerosol mass injection, and pyroCb trends over recent years (2017/2018-2020). Enhanced high-altitude injections of smoke aerosol mass (Fig. 4) suggests evolution of the vertical distribution of wildfire emitted aerosols in response to increasingly intense wildfire activity. In light of the climate changewildfire connection across the Western US 14-24 , we suspect a causal chain that, in part, implicates climate change in elevated wildfire plume top heights and enhanced aerosol injection aloft. Given the link between aerosol injection altitude and the corresponding air quality impacts, our findings suggest that wildfire activity in the Western US presents a growing risk in terms of long-range smoke transport and air quality degradation. Concurrent upward trends in the mass of biomass burning aerosols injected above and within the PBL present a threat to human health across local to continental scales. Recent wildfire activity stands testament to this, with the 2020 fire season likely producing thousands of smoke-related deaths 71 , and satellite imagery suggesting continental scale smoke transport akin to the documented presence of Western US/ Canadian wildfire smoke in New York City in August of 2018 74 .
From a climate perspective, these results point to Western US wildfire activity as becoming increasingly relevant to regional climate. Injection of a greater aerosol mass at higher altitudes, potentially into the upper troposphere-lower stratosphere, corresponds to an enhanced direct radiative effect given longer atmospheric lifetimes and increased absorption efficacy 9,29,31,32,68 . While less understood, there are also potential indirect effects related to aerosol-cloud microphysical interactions 12,13,34 . Though speculative, the potential for wildfire-derived aerosols to act as cloud condensation nuclei 13 and stave off initiation of rainfall 12,13 is further cause for concern given the risk of a positive feedback on wildfire activity. www.nature.com/scientificreports/ As caveats to the results presented within this study, we note that due to the multi-decadal timeframe considered and the associated computational cost, the model framework neglects the possibility of multi-core plume updrafts within an active fire 75 , as well as fire-atmosphere and aerosol-PBL feedbacks that can alter the PBL height and aerosol mass fraction contained within the PBL [76][77][78] . Additionally, we do acknowledge that vegetation type produces variability in the heat release per unit dry matter burned by wildfires. However, this heat release remains relatively constant between 16 and 22 MJ/kg as shown in pervious literature 79 . Furthermore, the F2010 model exhibits relatively low sensitivity to heat flux inputs 53 . Therefore, we suspect that the assumed heat release constant of 18.6 MJ/kg has not skewed the results presented here.
While the questions addressed herein provide a critical initial view of how the vertical distribution of wildfire emitted aerosols has evolved over the Western US and Canada over the past two decades, the presented hinge fit results suggest that this topic will need to be revisited in future studies. Further development of wildfire plume rise modeling and the inclusion of additional years of data should provide for an improved understanding of trends in wildfire aerosol injection. Furthermore, given global variability in wildfire activity 80 , it is apparent that investigation into plume height trends at the global scale would better elucidate the potential climate impacts of increasingly intense wildfire activity [14][15][16][17][18] . Such an investigation would be particularly valuable should plume height trends be considered within a radiative modelling framework. An expanded understanding of the evolving vertical distribution of wildfire emitted aerosols would hold direct relevance to air-quality and climate change projections, expression of wildfire emissions in near-field modeling efforts, and the accuracy of aerosol layer height assumptions required for satellite retrievals of absorbing aerosols.

Methods
Data. The spatial information, temporal information, and usage of datasets included in this study are presented in Table 1. Individual dataset descriptions can be found in Supplemental Information (SI) Sect. S2.

Models. Freitas et al. plume rise model. The F2010 model is a 1-dimensional, physics-based model that
relies on the vertical momentum equation, the first law of thermodynamics, and mass continuity for water phases. Processes related to lateral and dynamic entrainment are also parameterized within this model. The F2010 model estimates the wildfire plume top height as a function of wildfire heat flux density, active burn area, and environmental conditions ingested from an atmospheric transport model. By iteratively solving for a steady state plume rise solution, the F2010 model is able to approximate the plume top height and a thermodynamic profile of the final plume at 100 m vertical resolution. The plume top height is defined as the altitude at which the plume is neutrally buoyant, and is approximated as a vertical velocity < 1 m/s. Entrainment coefficients, output as part of the thermodynamic profile, can be used to estimate the vertical profile of aerosol detrainment below plume top when it is assumed that mass is conserved within the column for a steady state solution. While uncertainties remain in the model, F2010 was selected for its adherence to a physical basis and its ability to account for latent heat release by cloud microphysical processes 47,53 . Empirical plume rise schemes neglecting latent heat release through condensation and freezing are thus at a disadvantage relative to F2010. Evaluating the performance of the F2010 model is difficult given limited observations and the relatively poor constraint on wildfire heat flux and active burn area inputs. A breakdown of errors due to input uncertainty relative to the model structure is unknown at this time. In an evaluation of the F2010 model against MISR plume top observations, Val Martin et al. 45 take a range of approaches to model input, ultimately finding a maximum   86 , and sub-daily temporal downscaling profiles based on satellite active fire detections 54,80 .
WRF-atmospheric inputs. The environmental condition inputs to the F2010 model were supplied by a regional-decadal scale WRF simulation. For this study, we used an updated version of the Freitas model described and evaluated by Mallia et al. 57 , which can ingest WRF input files directly. The WRF simulations generated for this study used a nested domain setup where a 4 km resolution domain covered most of the Western US and was nested within a 12 km resolution domain covering western North America (Fig. S3). WRF simulations were reinitialized weekly using initial and boundary conditions from the Climate Forecast System Reanalysis 83 . The physical parameterizations used in these WRF simulations can be found in SI Sect. S3 and were based on configurations recommended by Mallia et al. 87 . Simulated plume top heights were not sensitive to the resolution of WRF domains as seen in the sensitivity analysis carried out in SI Sect. S3.
To evaluate the performance of these WRF simulations for surface winds, temperature, and humidity, we used METSTAT, a software package designed to evaluate meteorological inputs to air quality models 88 . Using METSTAT, WRF gridded variables were interpolated to the locations of surface meteorological observations from 1258 sites in the TDL U.S. and Canada Surface Hourly Observations dataset 81 to assess the performance of our WRF model simulations. We calculate the bias and root mean square error (RMSE) for surface wind speed, wind direction, temperature, and humidity (Table 2). When compared to benchmark values 89 , the WRF meteorology generated here outperforms the suggested benchmarks for surface wind speed bias, wind speed RMSE, wind direction bias, temperature bias, and humidity bias. Emery et al. 89 do not present RMSE benchmarks for surface wind direction, temperature, or humidity. The bias and RMSE of surface WRF variables compare favorably to those presented in related work by Mallia et al. 87 , though given overlap with the WRF setup used by Mallia et al. 87 , this likely owes to greater averaging of model errors over space and time.
The upper-air evaluation of WRF wind speed, wind direction, and potential temperature was carried out using data from rawinsonde launches at 24 sites throughout the Western US and Canada 82 . For each balloon launch, the vertical profiles of U, V, and potential temperature were extracted from the corresponding 12 km WRF grid cell at the time of the balloon launch. Horizontal winds simulated by WRF were then converted to wind speeds and wind directions, while vertical profiles of rawinsonde data were interpolated by pressure to match WRF vertical levels. To focus the evaluation on simulated upper air values and to avoid the effects of greater turbulence within the PBL, evaluation was limited to vertical levels situated above the simulated PBL height. Furthermore, our evaluation was only performed for rawinsonde observations at which the balloon had been advected less than 4 km from the original launch site. This 4 km limit was selected to balance between the need for ample observations to evaluate WRF output and the fact that the quality of the evaluation would be degraded when using observations of a balloon advected into adjacent WRF grid cells. The bias and RMSE of wind speed, wind direction, and potential temperature were calculated (Table 2).
Wildfire inputs. Gridded QFED PM 2.5 emissions (0.1° resolution) 61 are provided at daily resolution for four biome types (extratropical forests, grasslands, savannahs, and tropical forests) and were converted to a heat flux dataset using PM 2.5 emissions factors (PM 2.5 emitted/kg dry matter burned) 90 and assuming that the heat release per unit dry matter burned is 18.6 MJ/kg 57,91 . This conversion is depicted in Eq. (1), where F heat is the gridded heat flux, and PM 2.5,i and EF i are the gridded PM 2.5 flux and emission factor associated with biome type i , respectively. www.nature.com/scientificreports/ Heat flux datasets were then overlaid with daily burned area polygons from Moderate Resolution Imaging Spectroradiometer (MODIS) collection 6 to identify probable linkages between wildfire heat flux and satellite observed burned area. Since QFED emissions are based on MODIS observed thermal hotspots (MOD14A1/ MYD14A1 product), the literature suggests an expectation of collocated burned area and wildfire heat fluxes within a ± 8 day window (algorithm nominal uncertainty 92 ), with the majority of collocations occurring within ± 1 day 86,93 . Given this temporal window, we linked wildfire heat fluxes to burned area polygons using three complementary approaches that account for the temporal uncertainty in wildfire burned area observations and the idea of a directional fire spread. Specifics regarding this process are provided in SI Sect. S4, with the relevant methods depicted within a detailed visual (Fig. S5).
After developing relationships between gridded heat fluxes and burned area, we assumed that each cluster of continuous burned area polygons would produce its own plume. Within a given QFED grid cell, daily heat fluxes were then attributed to each intersecting cluster of burned area polygons based on the fraction of the total burned area within the QFED grid cell comprised by the burned area cluster. Burned area clusters intersecting multiple grid cells with non-zero heat fluxes were linked to heat fluxes from each grid cell based on the area of overlap within that grid cell. Ultimately, the heat flux value attributed to each burned area cluster was divided by the area of that cluster to estimate the heat flux density.
In total, 78% of the MODIS burned area was paired with a gridded heat flux, with this value increasing to 89% when considering only burned area polygons that intersect a grid cell with non-zero heat flux within the ± 8 day window. Similarly, 91% of QFED derived heat fluxes were linked to burned area polygons. Possible explanations for the imperfect collocation of heat fluxes and MODIS burned area polygons could stem from clouds obscuring satellite observation and differences between the spatial resolution of the MCD64A1 burned area product and the MOD14A1 thermal hotspot product used by QFED. Generally, years characterized by greater wildfire activity present greater collocation of heat fluxes and burned area (Table S2). We hypothesize that this relationship is the result of resolution limitations to observation, such that more intense wildfire activity is more readily observed by both the MCD64A1 and MOD14A1 products, increasing the odds of collocated observations. Temporal downscaling from daily to 3-h emissions and burned area was performed using diurnal cycles of wildfire activity based on gridded (0.25° resolution) and monthly aggregated satellite active fire detections provided by the Global Fire Emissions Database (GFED) 54,80 . For each set of plume inputs, the centroid of the associated burned area cluster was referenced to the GFED grid, and the relevant diurnal cycle was applied to both the heat flux and active burn area (such that the heat flux density remained unchanged). The resultant values were assumed to apply to each hour within the 3-h window, meaning additional plume top/thermodynamic variability within any 3-h downscaling window would strictly be a result of varied meteorological conditions. Entrainment-mass conservation. To estimate the vertical distribution of smoke detrainment below plume top, thermodynamic F2010 outputs for each plume were subject to vertical level-by-level mass conservation calculations following Eq. (2), where D mass,Z is the mass of in plume air detrained from vertical level Z , E mass,Z is the mass of environmental air entrained into the plume at vertical level Z , and the remaining terms characterize the vertical flux of mass into the layer from the layer below and the vertical flux of mass out of the layer to the layer above, respectively (left-to-right). The vertical mass flux between any two layers was taken to be a product of their mean density ( ρ ), mean vertical velocity ( W ), and mean plume horizontal area ( A ), ideally characterizing the interface between layers. At the plume top, it was assumed that the vertical flux to the above layer is zero, such that any remaining vertical mass transport is redirected as detrainment within the uppermost layer.
To express level-by-level detrainment values in terms of a vertical distribution of wildfire emitted aerosol, vertical profiles of detrainment were linearly interpolated from 100 m vertical resolution to 10 m vertical resolution, and then normalized. For interpolation, it is assumed that the detrainment at the surface and plume top boundaries is equal to 0. Resultant profiles effectively describe the fraction of total plume aerosol emissions that are injected at each 10 m vertical level.
PBL estimates from STILT. STILT calculates the PBL height using a modified Richardson number method that is applicable under stable, neutral, and unstable conditions 62,94 . Here, STILT is used to estimate the PBL height relevant to each modeled plume rise by initializing a model receptor on the centroid of the active burn area at the time of the plume rise. PBL heights are then extracted from the STILT output and linked with the plume climatology, allowing identification of penetrative plumes and quantification of the aerosol mass injected above the PBL. Meteorological inputs to STILT are sourced from the nested WRF simulations described above.
Trend analyses. Trend analyses aimed at understanding the evolution of plume top heights and the vertical distribution of wildfire generated aerosol mass, were performed for US level 3 ecoregions 66 and Canadian ecoprovinces 67 . For each ecoregion/ecoprovince West of − 100° E, the annual maximum plume top height, the PM 2.5 flux-weighted average plume top height, the PM 2.5 mass injected relative to threshold altitudes (> 18 km, > PBL, < PBL), and the total emissions were subjected to linear regression. Statistical significance was addressed via 100,000 bootstrapped repetitions. Each linear regression was reapplied to 100,000 resampled ver- (2) D mass,Z = E mass,Z + ρ Z,Z−1 W Z,Z−1 A Z,Z−1 − ρ Z,Z+1 W Z,Z+1 (A Z,Z+1 ) www.nature.com/scientificreports/ sions of the annually aggregated values, allowing replacement of selected samples, to determine if the regression to the actual data was distinguishable from the scatter produced by resampling. Hinge fitting, or change point analysis, was also used to elucidate potential emergent trend accelerations over recent years. Using the statistical analysis software R 95 , the "chngpt" package 96 was used to fit 2-segmented linear trends to PM 2.5 flux-weighted average plume top height, PM 2.5 mass injection relative to the PBL, and PM 2.5 emissions timeseries. Significant hinge-points were selected using a bootstrapping procedure that identifies the optimal hinge, with an associated p-value based on the quality of model fit. Hinge fitting has previously been used to compute climate normals amidst a changing climate 97 . Comparison to MISR plume top heights. MISR derived plume top heights were used to evaluate the performance of the F2010 model by comparing the distribution of modeled and observed plume top heights. Additionally, a plume-by-plume analysis was performed, which considers the deviation between modelled and observed plume tops in the context of instrument uncertainty. From the archived plumes maintained as part of the MISR Plume Heights Project 2, we found 480 plumes with a quality flag of "good" or "fair" within western North America (Fig. 5b). These observed plumes characterize the years 2008-2011 and 2017-2018. Using bin widths of 250 m, we developed a histogram comparing these 480 observed plumes to the model-derived late morning plumes (corresponding to MISR overpass times) spanning the same years (251,262 plumes; Fig. 3a). Distribution level comparisons between F2010 model output and MISR observations have previously identified a strong sensitivity to model inputs and a tendency for the model to underestimate the range of observed plume top heights 45 .
On a plume-by-plume scale, we found 425 MISR observed plumes within 30-min and 4 km of a modelderived plume. A cutoff distance of 4 km was used to compromise between the importance of local meteorology to plume top heights and the uncertainty in plume location as a result of using burned area centroids to denote plume locations. In cases when hourly resolution model outputs result in a plume top estimate both before and after the MISR overpass time, modeled plume top heights are temporally interpolated to the timing of the MISR overpass. Analysis of matched sets of plumes focusses on the standard deviation of modeled plume top heights from their observed counterparts and the fraction of modeled plumes falling within the MISR observational uncertainty.