Satellite-based estimation of the impacts of summertime wildfires on PM_2.5 concentration in the United States

Frequent and widespread wildfires in the northwestern United States and Canada have become the “new normal” during the Northern Hemisphere summer months, which significantly degrades particulate matter air quality in the United States. Using the mid-visible Multi Angle Implementation of Atmospheric Correction (MAIAC) satellite-derived aerosol optical depth (AOD) with meteorological information from the European Centre for Medium-Range Weather Forecasts (ECMWF) and other ancillary data, we quantify the impact of these fires on fine particulate matter concentration (PM_2.5) air quality in the United States. We use a geographically weighted regression (GWR) method to estimate surface PM_2.5 in the United States between low (2011) and high (2018) fire activity years. Our results indicate an overall leave-one-out cross-validation (LOOCV) R² value of 0.797 with root mean square error (RMSE) between 3 and 5 µg m⁻³. Our results indicate that smoke aerosols caused significant pollution changes over half of the United States. We estimate that nearly 29 states have increased PM_2.5 during the fire-active year and that 15 of these states have PM_2.5 concentrations more than 2 times that of the inactive year. Furthermore, these fires increased the daily mean surface PM_2.5 concentrations in Washington and Oregon by 38 to 259 µg m⁻³, posing significant health risks especially to vulnerable populations. Our results also show that the GWR model can be successfully applied to PM_2.5 estimations from wildfires, thereby providing useful information for various applications such as public health assessment.

1 Introduction

The United States (US) Clean Air Act (CAA) was passed in 1970 to reduce pollution levels and protect public health and has led to significant improvements in air quality (Hubbell et al., 2010; Samet, 2011). However, the northern part of the US continues to experience an increase in surface PM_2.5 due to fires in the northwestern United States and Canada (hereafter NWUSC), especially during the summer months, and these aerosols are a new source of “pollution” (Coogan et al., 2019; Dreessen et al., 2016). The smoke aerosols from these fires increase fine particulate matter (PM_2.5) concentrations and degrade air quality in the United States (Miller et al., 2011). Moreover, several studies have shown that from 2013 to 2016, over 76 % of Canadians and 69 % of Americans were at least minimally affected by wildfire smoke (Munoz-Alpizar et al., 2017). Although wildfire pre-suppression and suppression costs have increased, the number of large fires and the burnt areas in many parts of western Canada and the United States have also increased (Hanes et al., 2019; Tymstra et al., 2019). Furthermore, in a changing climate, as surface temperature increases and humidity decreases, the flammability of land cover also increases and thus accelerates the spread of wildfires (Melillo et al., 2014). The accumulation of flammable materials like leaf litter can potentially trigger severe wildfire events, even in those forests that hardly experience wildfires (Calkin et al., 2015; Hessburg et al., 2015; Stephens, 2005).

Wildfire smoke exposure can cause small particles to be lodged in lungs, which may lead to exacerbations of asthma chronic obstructive pulmonary disease (COPD), bronchitis, heart disease and pneumonia (Apte et al., 2018; Cascio, 2018). According to a recent study, a 10 µg m⁻³ increase in PM_2.5 is associated with a 12.4 % increase in cardiovascular mortality (Kollanus et al., 2016). In addition, exposure to wildfire smoke is also related to massive economic costs due to premature mortality, loss of workforce productivity, impacts on the quality of life and compromised water quality (Meixner and Wohlgemuth, 2004).

Surface PM_2.5 is one of the most commonly used parameters to assess the health effects of ambient air pollution. Given the sparsity of measurements in many parts of the world, it is not possible to use interpolation techniques between monitors to provide PM_2.5 estimates on a square kilometer basis. Since surface monitors are limited, satellite data have been used with numerous ancillary datasets to estimate surface PM_2.5 at various spatial scales. Several techniques have been developed to estimate surface PM_2.5 using satellite observations from regional to global scales, including simple linear regression, multiple linear regression, mixed-effect models, chemical transport models (scaling methods), geographically weighted regression (GWR), and machine-learning methods (see Hoff and Christopher, 2009, for a review). The commonly used global satellite data product is the 550 nm (mid-visible) aerosol optical depth (AOD), which is a unitless columnar measure of aerosol extinction. Simple linear regression methods use satellite AOD as the only independent variable, which shows limited predictability compared to other methods, and correlation coefficients vary from 0.2 to 0.6 from the western to eastern United States (Zhang et al., 2009). Multiple linear regression methods use numerous variables along with AOD data, and the prediction accuracy varies with different conditions, including the height of the boundary layer and other meteorological conditions (Goldberg et al., 2019; Gupta and Christopher, 2009b; Liu et al., 2005). For both univariate and multi-variate models, the AOD shows stronger correlation with PM_2.5 during fire episodes compared to pre-fire and post-fire periods (Mirzaei et al., 2018). Chemistry transport models (CTMs) that scale the satellite AOD by the ratio of PM_2.5 to AOD simulated by models can provide PM_2.5 estimations without ground measurements, which are different from other statistical methods (Van Donkelaar et al., 2006, 2019). However, the CTM models depend on reliable emission data, show limited predictability at shorter timescales and are largely useful for studies that require annual averages (Hystad et al., 2012). Different machine-learning methods, including neuron networks, random forests, and deep belief networks, show improvements in prediction accuracy (with CV R² values larger than 0.8), which is hard to accomplish for other parametric regression models (Gupta and Christopher, 2009a; Hu et al., 2017; Li et al., 2017; Wei et al., 2019, 2020, 2021). However, these methods also require a large number of samples to train the model, which means it is more suitable for daily PM_2.5 estimation rather than short-term wildfire events with relatively low occurrence frequency.

The relationship among PM_2.5, AOD and other meteorological variables is not spatially consistent (Hoff and Christopher, 2009; Hu, 2009). Therefore, methods that consider spatial variability can replicate surface PM_2.5 with higher accuracy. One such method is the GWR, which is a non-stationary technique that models spatially varying relationships by assuming that the coefficients in the model are functions of locations (Brunsdon et al., 1996; Fotheringham et al., 1998, 2003). In 2009, satellite-retrieved AOD was introduced in the GWR method to predict surface PM_2.5 (Hu, 2009), followed by the use of meteorological parameters and land use information (Hu et al., 2013). Meteorological variables are crucial for simulating surface PM_2.5 since they interact with PM_2.5 through different processes (Chen et al., 2020), which will be discussed in detail in the data section. Several studies (Guo et al., 2021; Ma et al., 2014; You et al., 2016a) successfully applied the GWR model in estimating PM_2.5 in China by using AOD and meteorological features as predictors. Similarly to all the statistical methods, however, the GWR relies on an adequate number and density of surface measurements (Chu et al., 2016; Gu, 2019; Guo et al., 2021), underscoring the importance of adequate ground monitoring of surface PM_2.5.

In this paper, we use satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) and surface PM_2.5 data combined with meteorological and other ancillary information to develop and use the GWR method to estimate PM_2.5. The use of the GWR method is not novel, and we merely use a proven method to estimate surface PM_2.5 from forest fires. We calculate the change in PM_2.5 between high fire activity (2018) and low fire activity (2011) periods during summer to assess the role of NWUSC wildfires in surface PM_2.5 in the United States. The paper is organized as follows. We describe the datasets used in this study, followed by the GWR method. We then describe the results and discussion, followed by a summary with conclusions.

2 Data

A 17 d period (9 to 25 August) in 2018 (high fire activity) and 2011 (low fire activity) was selected based on analysis of total fires (details in the methodology section) to assess surface PM_2.5 (Table 1).

Table 1 Datasets used in the study with sources.

https://www.epa.gov/outdoor-air-quality-data, https://earthdata.nasa.gov/, https://earthdata.nasa.gov/, https://www.ecmwf.int/en/forecasts (last access: 26 July 2021).

2.1 Ground-level PM_2.5 observations

Daily surface PM_2.5 from the Environmental Protection Agency (EPA) is used in this study. These data are from Federal Reference Methods (FRM), Federal Equivalent Methods (FEM), or other methods that are to be used in the National Ambient Air Quality Standards (NAAQS) decisions. A total of 1003 monitoring sites in the US are included in our study, with 949 having valid observations in the study period in 2018 and a total of 873 sites with 820 having valid observations in the study period in 2011. PM_2.5 values less than 2 µg m⁻³ are discarded since they are lower than the established lower detection limit (LDL) (EPA, 2011, 2018).

2.2 Satellite data

AOD, which represents the total column aerosol mass loading, is related to surface PM_2.5 as a function of aerosol properties (Koelemeijer et al., 2006):

$$\begin{array}{}\text{(1)}& \mathrm{AOD}={\mathrm{PM}}_{\mathrm{2.5}}Hf\left(\mathrm{RH}\right){\displaystyle \frac{\mathrm{3}{Q}_{\mathrm{ext},\mathrm{dry}}}{\mathrm{4}\mathit{\rho }{r}_{\mathrm{eff}}}}={\mathrm{PM}}_{\mathrm{2.5}}HS,\end{array}$$

where H is the aerosol layer height, f(RH) is the ratio of ambient and dry extinction coefficients, Q_ext,dry is the extinction efficiency under dry conditions, r_eff is the particle effective radius, ρ is the aerosol mass density and S is the specific extinction efficiency (m² g⁻¹) of the aerosol under ambient conditions. Although AOD usually has a positive correlation with PM_2.5, this relationship depends on various meteorological parameters which will be discussed in Sect. 2.3.

The MODIS mid-visible AOD from the Multi-Angle Implementation of Atmospheric Correction (MAIAC) product (MCD19A2 version 6) is used in this study. We used the MAIAC-retrieved Terra and Aqua MODIS AOD product at 1 km spatial resolution (Lyapustin et al., 2018), where different orbits are averaged to obtain mean daily values. Since thick smoke plumes generated by wildfires can be misclassified as clouds, only AOD less than 0 will be discarded. Validation with AERONET studies shows that 66 % of the MAIAC AOD data agree within ±0.5 to ±0.1 AOD (Lyapustin et al., 2018). Largely due to cloud cover, grid cells may have a limited number of AOD observations within a certain period. On average, cloud-free AOD data are available about 40 % of the time during 9 to 25 August 2018 when fires were active in the region bounded by 25–50^∘ N, 65–125^∘ W. A smoke flag from the same product is used as a predictor in estimating surface PM_2.5. The smoke detection is performed using MODIS red, blue and deep blue bands, and smoke pixels are separated from dust and clouds based on absorption parameter, size parameter and thermal thresholds (see Lyapustin et al., 2012, 2018, for further discussion). The smoke flag data can provide the percentage of a smoke pixel in each grid, which is related to smoke coverage.

We also use the MODIS level-3 daily FRP (MCD14ML, fire radiative power) product which combines the Terra and Aqua fire products to assess wildfire activity. The fire radiative energy indicates the rate of combustion, and thus FRP can be used for characterizing active fires (Freeborn et al., 2014). For the purposes of the study we sum the FRP within every 2.3^∘ × 3.5^∘ box to represent the total fire activity in different locations.

2.3 Meteorological data

Meteorological information, including boundary layer height (BLH), 2 m temperature (T2M), 10 m wind speed (WS), surface relative humidity (RH) and surface pressure (SP), is obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA5) product, with a spatial resolution of 0.25^∘ and temporal resolution of 1 h, and is matched temporally with the satellite overpass time. The meteorological parameters provide important information on different processes affecting surface PM_2.5 concentration, which can also be seen as supplements of the AOD–PM_2.5 relationship as previously discussed.

The BLH can provide information on aerosol layer height (H in Eq. 1) as aerosols are often found to be well mixed within the boundary layer (Gupta and Christopher, 2009b). With the same amount of pollution within the boundary layer, the higher the BLH, the more PM_2.5 is distributed within that layer and vice versa (Miao et al., 2018; Zheng et al., 2017). Therefore, PM_2.5 usually has an anticorrelation with BLH. However, for wildfire events, the aerosol layer height is sometimes higher than the BLH (Haarig et al., 2018), which leads to lower correlation between AOD and PM_2.5 since we use only BLH to present the aerosol layer height. Thus BLH can provide aerosol vertical information in most cases except for suspended high-layer aerosol caused by fires, which leads to higher bias of the model for high-layer aerosols near the fire sources. Surface temperature (T2M) can affect PM_2.5 through convection, evaporation, temperature inversion and secondary pollutant generation processes (Chen et al., 2020). The first two processes are negatively related to PM_2.5 concentration: (1) higher temperature increases turbulence, which accelerates the dispersion of pollution and a decrease in PM_2.5; (2) higher temperature increases evaporation loss of PM_2.5, including ammonium nitrate and other volatile or semi-volatile components (Wang et al., 2017). The latter two processes are positively related to PM_2.5 by limiting vertical motion and promoting photochemical reactions under high temperature (Xu et al., 2019; Zhang et al., 2015). High wind speeds are often negatively related to PM_2.5 since it increases the dispersion of pollutants. However, unique geographical conditions such as mountains with certain wind directions can cause accumulation of pollutants (Chen et al., 2017). RH may promote hygroscopic growth of particles to increase PM_2.5 (Trueblood et al., 2018; Zheng et al., 2017). Changes in surface pressure (SP) may influence the diffusion or accumulation of pollutants through formation of low-level wind convergence (You et al., 2017). Precipitation is another factor that largely influences surface PM_2.5 since it can accelerate the wash-out of suspended particles, but AOD values are not available when clouds are present.

3 Methodology

To assess the impact of NWUSC fires on PM_2.5 in the United States, we first estimate PM_2.5 over the study region during a time period with high fire activity (2018). We then use the same method during a year with low fire activity (2011) to compare the differences between the two years. The two years are selected based on the total FRP in August calculated within Canada (49–60^∘ N, 55–135^∘ W) and the northwestern (NW) United States (35–49^∘ N, 105–125^∘ W). Table 2 shows the total FRP in Canada and the northwestern US in August from 2010 to 2018. The total FRP in the two regions is lowest in 2011 and highest in 2018 during the nine years, which provides the basis for the study. In order to create a 0.1^∘ surface PM_2.5, the GWR model is used to estimate the relationships of PM_2.5 and AOD using various parameters. Detailed processing steps for the GWR model are shown in Fig. 1.

Figure 1 Flow chart for the geographically weighted regression model used. All satellite, ground, and meteorological data are gridded to 0.1 by 0.1^∘.

Table 2 Total FRP in Canada and the northwestern US in August of different years (unit: 10⁴ MW).

3.1 Data preprocessing

The first step is to resample all datasets to a uniform spatial resolution by creating a 0.1^∘ resolution grid covering the continental United States. During this process, we collocate the PM_2.5 data and average the values if there is more than one value in one grid. Then the MAIAC AOD and smoke flag are averaged into 0.1^∘ grid cells. Meteorological datasets are also resampled to the 0.1^∘ grid cells by applying the inverse distance method.

3.2 Time selecting and averaging

Next we select data where AOD and ground PM_2.5 are both available (AOD > 0 and PM_2.5 > 2.0 µg m⁻³) and average them for the study period (since the LDL for the FRM method is 2 µg m⁻³ in 2011 and 3 µg m⁻³ in 2018, we decided to use the LDL for 2011) (EPA, 2011, 2018). This is to ensure that the AOD, PM_2.5 and other variables match with each other, because PM_2.5 is not a continuous measurement for some sites and AOD has missing values due to cloud cover and other reasons. Therefore, it is important to use data from days where both measurements are available to avoid sampling biases.

3.3 GWR model development and validation

The adaptive bandwidth selected by the Akaike information criterion (AIC) is used for the GWR model (Loader, 1999). For locations that already have PM_2.5 monitors, we calculate the mean AOD of a 0.5 × 0.5^∘ box centered at the ground location and estimate the GWR coefficients (β) for AOD and meteorological variables to estimate PM_2.5. The model structure can be expressed as

$$\begin{array}{rl}{\mathrm{PM}}_{\mathrm{2.5}i}& ={\mathit{\beta }}_{\mathrm{0},i}+{\mathit{\beta }}_{\mathrm{1},i}{\mathrm{AOD}}_i+{\mathit{\beta }}_{\mathrm{2},i}{\mathrm{BLH}}_i+{\mathit{\beta }}_{\mathrm{3},i}\mathrm{T}\mathrm{2}{\mathrm{M}}_i\\ & +{\mathit{\beta }}_{\mathrm{4},i}\mathrm{U}\mathrm{10}{\mathrm{M}}_i+{\mathit{\beta }}_{\mathrm{5},i}{\mathrm{RH}}_{sfci}+{\mathit{\beta }}_{\mathrm{6},i}{\mathrm{SP}}_i+{\mathit{\beta }}_{\mathrm{7},i}{\mathrm{SF}}_i+{\mathit{\epsilon }}_i,\end{array}$$

where PM_2.5i (µg m⁻³) is the selected ground-level PM_2.5 concentration at location i; β_0,i is the intercept at location i; ${\mathit{\beta }}_{\mathrm{1},i}\sim {\mathit{\beta }}_{\mathrm{8},i}$ are the location-specific coefficients; AOD_i is the resampled AOD selected from MAIAC daily AOD data at location i; BLH_i, T2M_i, U10M_i, RH_sfci, and SP_i are selected meteorological parameters (BLH, T2M, WS, RH, and PS) at location i; SF_i (%) is the resampled smoke flag data at location i, and ε_i is the error term at location i.

We perform leave-one-out cross-validation (LOOCV) to test the model predictive performance (Kearns and Ron, 1999). Since the GWR model relies on an adequate number of observations, the prediction accuracy will be lower if we preserve too many data for validation. Therefore, we choose the LOOCV method, which preserves only one datum for validation at a time and repeats the process until all the data are used. In addition, R² and root mean square error (RMSE) are calculated for both model fitting and model validation processes to detect overfitting, which leads to low predictability.

3.4 Model prediction

While predicting the ground-level PM_2.5 for unsampled locations, we make use of the estimated parameters for sites within a 5^∘ radius to generate new slopes for independent variables based on the spatial weighting matrix (Brunsdon et al., 1996). The closer to the predicted location, the closer to 1 the weighting factor will be, while the weighting factor for sites farther than 5^∘ in distance is zero. It is important to note that AOD and other independent variables used for prediction in this step are averaged values for days that have valid AOD, which is different from the data used in the fitting process since PM_2.5 is not measured every day in all locations.

4 Results and discussion

We first discuss the surface PM_2.5 for a few select locations that are impacted by fires followed by the spatial distribution of MODIS AOD and the FRP for August 2018. We then assess the spatial distribution of surface PM_2.5 from the GWR method. The validation of the GWR method is then discussed. To further demonstrate the impact of the NWUSC fires on PM_2.5 air quality in the United States, we show the spatial distribution of the difference between August 2018 and August 2011 and quantify these results for 10 US EPA regions.

4.1 Descriptive statistics of satellite data and ground measurements

The 2018 summertime Canadian wildfires started around the end of July in British Columbia and continued until mid-September. The fires spread rapidly to the south of Canada during August, causing high concentrations of smoke aerosols to drift down to the US and affecting particulate matter air quality significantly. From late July to mid-September, wildfires in the northwestern US that burnt forest and grassland also affected air quality in the United States. Starting with the Cougar Creek Fire and then the Crescent Mountain and Gilbert fires, different wildfires in NWUSC caused severe air pollution in numerous US cities. Figure 2a shows the rapid increase in PM_2.5 in selected US cities from 1 July to 31 August, due to the transport of smoke from these wildfires. For all sites, July had low PM_2.5 concentrations (<10 µg m⁻³) and rapidly increases as fire activity increases. Calculating only from the EPA ground observations, the mean PM_2.5 of the 17 d for the entire US is 13.7 µg m⁻³ and the mean PM_2.5 for Washington (WA) is 40.6 µg m⁻³, which indicates that the PM pollution is concentrated in the northwestern US for these days. This trend is obvious when comparing the mean PM_2.5 of all US stations (black line with no markers) and the mean PM_2.5 of all WA stations (grey line with no markers). Ground-level PM_2.5 reaches peak values between 17 and 21 August, and daily PM_2.5 values during this time period far exceed the 17 d mean PM_2.5. For example, mean PM_2.5 in Washington on 20 August is 86.75 µg m⁻³, which is more than 2 times the 17 d average of this region. On 19 August, Omak, which is located in the foothills of the Okanogan Highlands in WA, had PM_2.5 values that exceed 250 µg m⁻³. According to a review of US wildfire-caused PM_2.5 exposures, 24 h mean PM_2.5 concentrations from wildfires ranged from 8.7 to 121 µg m⁻³, with a 24 h maximum concentration of 1659 µg m⁻³ (Navarro et al., 2018).

Figure 2 (a) Variations of EPA ground-observed PM_2.5 in different cities from July to August 2018 (Omak – Washington, Seattle – Washington, Chicago – Illinois, Portland – Oregon, Billings – Montana). Black line without markers shows the mean variation of all the US stations, and the grey line without markers shows the mean variation of stations in Washington State. (b) Mean MAIAC satellite AOD distribution from 9 to 25 August 2018. AOD values equal to or larger than 0.5 are shown as the same color (yellow). Also shown are circles with FRP. Black arrow shows the wind direction, and the length of it represents the wind speed. The round spots of different colors on the map show the locations of the five selected cities (green – Omak, black – Seattle, yellow – Chicago, blue – Portland, red – Billings).

Table 3 shows relevant statistics for 15 states that have at least one daily record of non-attainment of EPA standards (>35 µg m⁻³). From the frequency records of non-attainment in the 17 d period (last column), four states (Montana, Washington, California and Idaho) were consistently affected by the wildfires, and large portions of ground stations in these states were influenced by smoke aerosols. Most of the neighboring states also experienced significant air pollution (third column). Noticeable from these records is that the total number of ground stations in some of the highly affected states (such as Idaho) is not sufficient for capturing the smoke. Although there are a total of eight EPA stations in Idaho, only two of them have consistent observations during the fire event; the other two stations have no valid observations, and the remaining four stations have only two to six observations during the 17 d period. Limited availability of valid data along with unevenly distributed stations makes it hard to quantify smoke pollution in the northwestern US during the fire event period. Therefore, we utilize satellite data to enlarge the spatial coverage and estimate pollution at a finer spatial resolution.

Table 3 Statistics of 15 states that violate EPA standards (35 µg m⁻³) during the 17 d wildfire period.

The spatial distribution of AOD shown in Fig. 2b indicates that the smoke from Canada is concentrated mostly in northern US states such as Washington, Oregon, Idaho, Montana, North Dakota and Minnesota. The black arrow shows the mean 800 hPa-level mean wind for 17 d, and the length of the arrow represents the wind speed in m s⁻¹. Also shown in Fig. 2b are wind speeds close to the fire sources which are about 4–5 m s⁻¹, and according to the distances and wind directions, it can take approximately 28–36 h for the smoke to be transported southeastward to Washington. Then the smoke continues to move east to other northern states such as Montana and North Dakota. In addition, the grey circle represents the total FRP of every 2.3×3.5^∘ box. The reason for not choosing a smaller grid for the FRP is to not clutter Fig. 2b with information from small fires. The bigger the circle is, the stronger the fire is in that grid, and different sizes and its corresponding FRP values are shown in the lower-right corner. It is clear that the strongest fires in 2018 are located in Tweedsmuir Provincial Park of British Columbia in Canada (53.333^∘ N, 126.417^∘ W). The four separate lightning-caused wildfires burnt nearly 301 549 ha of the boreal forest. The total FRP of August 2018 in Canada is about 5362 (⋅ 1000 MW), while the total FRP of August 2011 in Canada is 48 (⋅ 1000 MW). The 2011 fire was relatively weak compared to the 2018 Tweedsmuir Complex fire, and we therefore use the 2011 air quality data as a baseline to quantify the 2018 fire influence on PM_2.5 in the United States.

4.2 Model fitting and validation

The main goal for using the GWR model is to help predict the spatial distribution of PM_2.5 for places with no ground monitors while leveraging the increased spatial resolution of satellite AOD, and therefore it is important to ensure that the model is robust. Figure 3a and b show the results for 2018 for GWR model fitting for the entire US and the LOOCV models, respectively. The color of the scatter plots represents the probability density function (PDF) which calculates the relative likelihood that the observed ground-level PM_2.5 would equal the predicted value. The lighter the color is, the more points are present, with a higher correlation. The model fitting process estimates the slope for each variable, and therefore the model can be fitted close to the observed PM_2.5, and using this estimated relationship, we are able to assess surface PM_2.5 using other parameters at locations where PM_2.5 monitors are not available. The LOOCV process tests the model performance for predicting PM_2.5. If the results of LOOCV have a large bias from the model fitting, then the predictability of the model is low. Higher R² and RMSE differences indicate that the model is overfitting and therefore not suitable. The R² for the model fitting is 0.834, and the R² for the LOOCV is 0.797, while the RMSE for the GWR model fitting is 3.46 µg m⁻³, and for LOOCV the RMSE is 3.84 µg m⁻³. There are minor differences between fitting R² and validation R² (0.037) and between fitting RMSE and validation RMSE (0.376 µg m⁻³), suggesting that the model is not overfitting and has stable predictability, further indicating that the model can predict surface PM_2.5 reliably. In addition, we also performed a 20-fold cross-validation by splitting the dataset into 20 consecutive folds, and each fold is used for validation, while the 19 remaining folds form the training set. The 20-fold cross-validation has an R² of 0.745 and a RMSE of 4.3 µg m⁻³. The increase/decrease in the cross-validated R² and RMSE indicates that sufficient data are used for fitting since a small decrease in the number of fitting data can reduce the model prediction accuracy. Overall, the prediction error of the model is between 3 and 5 µg m⁻³, which is a reasonable error range for 17 d average prediction of PM_2.5. For data greater than the daily mean EPA standard (35 µg m⁻³), the model has a RMSE of 12.07 µg m⁻³, which is a lot larger than the RMSE when using the entire model. Therefore, the model has a tendency to underestimate PM_2.5 exceedances by around 12.07 µg m⁻³. The larger the PM_2.5 is, the more the model underestimates. To examine the model performance for high- and low-pollution areas, the results are divided into two parts (larger than 35 µg m⁻³ and less than 35 µg m⁻³). Areas with high pollution have an R² of 0.64 and areas with low pollution have an R² of 0.67; therefore, the model performance is relatively stable for both large and small PM_2.5 values. Also, the inclusion of low aerosol concentration areas does not influence the model performance for high values (seen in the Supplement in Figs. S1 and S2), which means that the high R² is not a reason for the large number of low values. The GWR model fitting and validation results for the 17 d in 2011 are shown in Fig. S3.

Figure 3 Results of model fitting and cross-validation for the GWR model for the entire US region averaged from 9 to 25 August 2018. (a) GWR model fitting results; (b) GWR model LOOCV results. The dashed line is the 1 : 1 line as a reference and the black line shows the regression line. The color of the scatter plots represents the probability density function which provides a relative likelihood that the value of the random variable would equal a certain sample.

4.3 Predictors' influence during wildfires

Table 4 shows the GWR model mean coefficients for the whole US region and for different selected regions. The selected boxes are shown in Fig. 4c in different colors: box1 (red) located in the NW US includes major fire sources in the US; box2 (gold) located in Montana is influenced by both neighboring states and smoke from Canada; box3 (green) in Minnesota is located further from the fires and has a minor increase in PM_2.5 due to remote smoke; box4 (black) in the NE (northeastern) US is furthest from the fires and has no obvious pollution increase. The second column of the tables shows the conditions for sample selection, and the third column shows the number of pixels selected for each box. By comparing the coefficients of samples selected in these boxes, predictors have different influences in different locations. AOD has a stronger influence on predicting PM_2.5 closer to fire sources, but local emissions become more dominant if the distances are large enough. The smoke flag is overall positively related to surface PM_2.5, while it could slightly negatively relate to PM_2.5 around fire sources and the northeastern coasts. PBL is negatively related to PM_2.5 when the pollution is concentrated near the surface (fires or human-made emissions), while it appears to be positively related to PM_2.5 in locations where the main pollution source comes from remote wildfire smoke. Surface temperature has a relatively stable positive correlation with surface PM_2.5; however, surface pressure and wind speeds are negatively correlated with PM_2.5. Relative humidity, on the other hand, shows large variations in PM_2.5 influence across the nation. Around the wildfires where the RH is relatively low, RH has a positive correlation with PM_2.5 since hygroscopicity would increase and leads to accumulation of PM_2.5, but increasing RH can also decrease PM_2.5 concentration by overgrowing the PM_2.5 particles to deposition in a high-RH environment (Chen et al., 2018).

Figure 4 (a) EPA ground-observed PM_2.5 distribution over the US averaged from 9 to 25 August 2018. (b) GWR-predicted 17 d mean PM_2.5 distribution. (c) Difference map of predicted ground PM_2.5 of the 17 d mean values between 2018 and 2011. PM_2.5 values equal to or larger than 60 µg m⁻³ are shown as the same color (red). Note that the D-PM_2.5 has a different color scale to make the negative values more apparent (blue).

Table 4 Coefficients of different predictors.

From Table 4, we know that the weighting for AOD is much larger than other predictors, but predictors other than AOD are important for the prediction. We tested our model with AOD as the only predictor to conduct a comparison to the original model, and the R² decreases from 0.83 to 0.79 and RMSE increases from 3.46 to 3.8. This is consistent with a previous study (Jiang et al., 2017) which shows improvements of R² from 0.69 to 0.78 and RMSE from 7.25 to 6.18 by adding four meteorological parameters in summer in eastern China. Other predictors have higher weighting in the fire source region (box1), where BLH cannot provide the aerosol vertical distribution information since smoke tends to be injected to higher levels. For high AOD regions where aerosol tends to be suspended at high levels, adding predictors other than AOD tends to have lower improvement of the model compared to low AOD values, because adding BLH can significantly improve the prediction for low-level aerosols. For regions with AOD less than 35, R² increases by 0.09 from the AOD-only model (0.6 to 0.69), while R² increases by 0.05 for areas with AOD larger than 35. RMSE decreases by 12 % and 7 % for AOD less and larger than 35 conditions, respectively. Overall, the meteorological factors have larger improvements for low-pollution areas (low-level aerosol in this case).

4.4 Predicted PM_2.5 distribution

The mean PM_2.5 distributions over the United States shown in Fig. 4a are calculated by averaging the surface PM_2.5 data from ground monitors for the 17 d, which matches well with the GWR model-predicted PM_2.5 distributions shown in Fig. 4b. The model estimation extends the ground measurements and provides pollution assessments across the entire nation. Comparing the AOD map (Fig. 2b) with the PM_2.5 estimations (Fig. 4b) demonstrates the differences between columnar and surface-level pollution. Differences between the AOD and PM_2.5 distributions are for various reasons, including (1) areas with high PM_2.5 concentrations in Fig. 4b corresponding to low AOD values in Fig. 2b (southern California, Utah, and the southern US) and (2) and high AOD regions in Fig. 2b corresponding to low PM_2.5 concentrations in Fig. 4b (Minnesota). The first situation usually occurs at the edge of polluted areas that are relatively far from the fire source, which is consistent with previous studies that reported smaller particles (< 10 µg) being able to travel longer distances compared to large particles (> 10 µg) (Gillies et al., 1996) and larger particles tending to settle closer to their source (Sapkota et al., 2005; Zhu et al., 2002).

We use the same method for 9 to 25 August in 2011 that had low fire activity, ensuring consistency in estimating coefficients for different variables for 2011. Figure 4c shows the difference in spatial distribution of mean ground PM_2.5 of the 17 d between 2018 and 2011. Larger differences in PM_2.5 are in the northwestern and central parts of the United States, with the southern states having very little impact due to the fires. Of all 48 states within the study region, there are 29 states that have a higher PM_2.5 value in 2018 than 2011, and 15 states have a 2018 PM_2.5 value of more than 2 times their 2011 value (shown in Fig. 5). The mean PM_2.5 for WA increases from 5.87 in 2011 to 46.47 µg m⁻³ in 2018, which is about 8 times more than 2011 values. The PM_2.5 values in Oregon increase from 4.97 in 2011 to 33.3 µg m⁻³ in 2018, which is nearly a 7-fold increase. For states from Montana to Minnesota, the mean PM_2.5 decreases from east to west, which reveals the path of smoke transport. As shown in Fig. 4c, there is a clear transport path of smoke from North Dakota all the way to Texas. Along the path, smoke increases PM_2.5 concentrations by 168 % in North Dakota and 27 % in Texas. Smoke aerosols transported over long distances typically contain fine-fraction PM, which significantly affects the health of children, adults, and vulnerable groups.

Figure 5 Mean PM_2.5 from 9 to 25 August 2018 and 2011 of most affected states.

Figure 6 Mean PM_2.5 of EPA regions from 9 to 25 August 2011 and 2018. Inset shows the map of 10 EPA regions in different colors. Yellow column represents the 2018 mean PM_2.5, and green column represents for 2011 mean PM_2.5.

Figure 6 shows the mean PM_2.5 predicted from the GWR model of different EPA regions for the 17 d in 2011 and 2018 (Hawaii and Alaska are not included). The most influenced region is region 10, which has a 2018 mean PM_2.5 value of 34.2 µg m⁻³ that is 6 times larger than the values in 2011 (5.8 µg m⁻³). The PM_2.5 of regions 8 and 9 has 2.4 and 2.6 times increases in 2018 compared to 2011. Regions 1–4 have lower PM_2.5 in 2018 than 2011, possibly due to Clean Air Act initiatives, absence of any major fire activities and being further away for transported aerosols. The emission reduction improves the US air quality and lowers the PM_2.5 every year, but 6 out of the 10 EPA regions show significant increases in PM_2.5 during the study period, which indicates that the long-range transported wildfire smoke has become the new major pollutant in the US.

4.5 Estimation of Canadian fire pollution

To evaluate the pollution caused only by Canadian fires, we did a rough assessment according to the total FRP and PM_2.5 values. There are three states in the US that have wildfires during the study period, California, Washington and Oregon, and they have total FRPs of 1186, 518 and 439 (⋅ 1000 MW), respectively. Assuming that California was only influenced by the local fires, then fires of 1186 (⋅ 1000 MW) cause a 13 µg m⁻³ increase in PM_2.5. Accordingly, wildfires in Washington State and Oregon State will cause 6 and 5 µg m⁻³ increases in state mean PM_2.5. Therefore, Canadian fires caused PM_2.5 increases in Washington and Oregon of about 35 and 23 µg m⁻³, respectively. Since the FRP of Canadian wildfires is approximately 5 times larger than that of the California fires, which are the strongest fires in the US, we assume that the pollution affecting the states located in the downwind directions other than the three states is mainly coming from Canadian wildfires. States with no local fires such as Montana, North Dakota, South Dakota and Minnesota have PM_2.5 increases of 18.31, 12.8, 10.4 and 10.13 µg m⁻³. The decrease in these numbers reveals that the smoke is transported in a southeasterly direction. This influence of Canadian wildfires on US air quality is only a rough quantity estimation, and thus additional work is needed to understand long-range transport smoke pollution and its impact on public health.

4.6 Comparison to previous studies

Compared to the Bayesian ensemble model developed by Geng et al. (2018) using MAIAC AOD and CMAQ (Community Multiscale Air Quality) model and ground PM_2.5 measurements, our GWR model has larger R², but with the CTM, their method can provide more vertical distribution information, which is important for wildfire smoke. GWR usually has better accuracy than the CTM since there are large uncertainties related to different CTM inputs such as emission and meteorological and land cover data, but for regions with fewer or no ground measurements, the CTM provides a good approach for estimating surface PM_2.5. Other studies which used machine-learning methods to predict surface PM_2.5 have better performance for long-term prediction rather than monthly estimations (Liang et al., 2020; Xiao et al., 2018) but can better resolve complex relationships between different predictors than statistical models (Geng et al., 2020). For wildfire events, the available data are much fewer than the long-term aerosol analysis, so the performance of a machine-learning method could be less accurate compared to long-term prediction. Our study also shows slightly larger R² compared to other GWR studies (Hu et al., 2013; Ma et al., 2014; You et al., 2016b) due to the inclusion of more meteorological and other related predictors.

4.7 Model uncertainties and limitations

There are various sources of uncertainties and limitations for studies that use satellite data to estimate surface PM_2.5 concentrations. Since wildfires develop quickly, it is important to have continuous observations to capture the rapid changes. This study uses polar-orbiting high-quality satellite aerosol products, but the temporal evolution can only be estimated by geostationary datasets. Although satellite observations have excellent spatial coverage, missing data due to cloud cover are a limitation. As discussed in the paper, the prediction error (RMSE) of the model is between 3 and 5 µg m⁻³, while the RMSE increased for locations with high aerosol concentration. This is partly due to lack of accurate vertical distribution information, which is very important for wildfire smoke. The GWR model is largely influenced by the distribution of ground stations, and the prediction error will be different in different locations due to unevenly distributed PM_2.5 stations. For locations that have a dense ground-monitoring distribution, the prediction error will be low, while the prediction error will be relatively larger in other places with sparse surface stations. Although there are obvious limitations, complementing surface data with satellite products and meteorological and other ancillary information in a statistical model like the GWR has provided robust results for estimating surface PM_2.5 from wildfires. We also note that we did not consider some variables used in other studies, such as NDVI, forest cover, vegetation type, industrial density, visibility and chemical constituents of smoke particles (Van Donkelaar et al., 2015; Hu et al., 2013; You et al., 2015; Zou et al., 2016). Visibility mentioned in some studies may improve the model performance, but unlike AOD, it has limited measurement across the nation, which will restrict the applicability of training data. Another uncertainty comes from the 2011 wildfires, which we assumed to be zero fire events, but there are actually few fire events in EPA regions 6, 8, 9 and 10, and this will lead to underestimation of PM_2.5 increase due to the 2018 fires in these regions.

One limitation of this study is that analysis based on 17 d mean values cannot capture daily pollution variations, which is also very important for pollution estimation during rapidly changing wildfire events. To extend this analysis to daily estimation, the cloud contaminations of satellite observations become a major problem. Therefore, future work is needed using chemistry transport models and other data to fill in the gaps on missing AOD data due to cloud coverage.

5 Summary and conclusions

We estimate the surface mean PM_2.5 for 17 d in August for a high fire activity year (2018) and compare that to a low fire activity year using the geographically weighted regression (GWR) method to assess the increase in PM_2.5 in the United States due to smoke transported from fires. The difference in PM_2.5 between the two years indicates that more than half of the United States (29 states) is influenced by the NWUSC wildfires, and half of the affected states have 17 d mean PM_2.5 increases larger than 100 % of the baseline value. The peak PM_2.5 during the wildfires can be much larger than the 17 d average and can affect vulnerable populations susceptible to air pollution. Some of the most affected states are Washington, California, Wisconsin, Colorado and Oregon, all of which have populations greater than 4 million. According to the Centers for Disease Control and Prevention (CDC), 8 % of the population has asthma (CDC, 2011). Therefore, for asthma alone, there are about 3 million people facing significant health issues due to the long-range transport of smoke in these states.

For states that show a decrease in PM_2.5 due to the Clean Air Act, the mean decrease is about 16 % of the baseline after 7 years. This is consistent with the EPA's report that there is a 23 % decrease in PM_2.5 on national average from 2010 to 2019 (U.S. Environmental Protection Agency, 2019). Compared to the dramatic increase (132 %) caused by wildfires, pollution from the fires counteracts our effort on emission controls. Although wildfires are often episodic and short term, high frequency of fire occurrence and increasingly longer durations of summertime wildfires in recent years have made them now a long-term influence on public lives. Our results show a significant increase in pollution in a short time period in most of the US states due to the NWUSC wildfires, which affects millions of people. With wildfires becoming more frequent during recent years, more effort is needed to predict and warn the public about the long-range transported smoke from wildfires.