Estimating PM_2.5 in Southern California using satellite data: factors that affect model performance

The study area encompasses several counties in southern California and includes the metropolitan areas of Los Angeles, Long Beach, Riverside, and San Diego. Our modeling domain measures approximately 460 km × 245 km adjacent to the U.S.-Mexico border, with a total population of more than 15 million. It contains a mixture of terrain including coastal lowlands, highlands, mountain valleys, and both lowland and highland deserts. Additionally, several areas in the region tend to have higher concentrations of pollution, including a few of the country's most polluted cities [46]. Some of these areas can be found in inland urban centers (such as Riverside) and constrained valleys (i.e. Imperial Valley). A map of the study area with corresponding cities, EPA ground monitoring stations, and underlying topography is shown in figure 1.

Zoom In Zoom Out Reset image size

All 24-hour average PM_2.5 and PM₁₀ concentrations in 2012 were acquired from the EPA Air Quality System Data Mart (http://www.epa.gov/ttn/airs/aqsdatamart) [47]. A variable representing the ratio between PM_2.5 and PM₁₀ concentrations was calculated for use in the first stage of the modeling process to account for presence of dust and other coarse particles. To date, PM10 has yet to be regularly included in satellite-driven PM2.5 exposure models. Although MODIS AOD is generally most sensitive to smaller particles generally best characterized as PM2.5, coarse particles may also scatter or absorb light. As shown in this analysis, the southern California domain appears to have consistently high PM10 levels that warrant its inclusion as a predictor in our model to enhance PM2.5 predictions [48–51] Primary PM_2.5 emissions (tons per year) and the number of major point sources for the study area were acquired from the EPA National Emissions Inventory and the total emissions were summed by grid cell. The aerosol optical depth (AOD) data at 1 km spatial resolution were extracted from the NASA MODIS Multi-Angle Implementation of Atmospheric Correction (MAIAC, MCD19A2) product (https://search.earthdata.nasa.gov/) [31, 52, 53]. The meteorological data such as air temperature, wind speed, and relative humidity were extracted from the North American Land Data Assimilation System Phase 2 (NLDAS-2) at 1/8-degree grid (∼ 12 km) resolution [54]. Planetary boundary layer height data were derived from the analysis fields of the North American Regional Reanalysis (NARR) at ∼ 32 km spatial resolution and 3-hour temporal frequency. NARR fields were spatially interpolated to the 1 km grid and then temporally disaggregated to the NLDAS hourly frequency [55].

Elevation was derived from the 3-arc-second (90-meter) Shuttle Radar Topography Mission (SRTM) dataset distributed by USGS Earth Resources Observation and Science (EROS) Data Center (https://www.usgs.gov/centers/eros). Additional land cover variables, including forest cover and impervious surfaces, were retrieved from the 2011 National Land Cover Database (NLCD, https://catalog.data.gov/dataset/usgs-2011-national-landcover). The spatial resolution of the NLCD coverage is 30 × 30 m². Coverage maps were generated for forest pixels (pixels assigned 1 for forest and 0 for non-forest) and for percent imperviousness across the study area. Additional distance variables were included to account for potential effect modifiers unique to the region, including distance to the coastline and distance to Mexicali, Mexico. Road length data were obtained from ESRI StreetMap USA (Environmental Systems Research Institute, Inc. Redland, CA). The sum of the road segment lengths was determined in ArcGIS for each modeling grid cell. All model input data were mapped to the 1 km modeling grid using a spatial averaging procedure. Each PM_2.5 monitoring site was matched to the nearest cell AOD, temperature, relative humidity, and wind speed. Land use variables were either averaged (forest cover and elevation) or summed (emissions, roads, point sources).

We calibrated the relationship between PM_2.5 and AOD using a two-stage modeling framework, which allowed the relationship to vary in both space and time [33, 56, 57]. In the first stage, a linear mixed effects model (LME) was utilized with daily random slopes and intercepts for AOD, relative humidity, and wind speed. Since each of these variables are time varying, their inclusion as random effects aids in accounting for any temporal variation in the overall relationship between AOD and PM_2.5. In addition to the random effects terms, the model included several fixed effects. Fixed effects in the model help to estimate the mean values and the included random effects help to account for daily variability in the relationship between dependent and independent variables. We considered multiple land use and meteorological predictors during the model selection process; however, we chose to eliminate some of the predictors from the final model due to lack of significance. The first stage of the model can be expressed by the following:

$$\begin{align} P{M_{2.5,sd}} & = \left( {{b_0} + {b_{0,d}}} \right) + \left( {{b_1} + {b_{1,d}}} \right)AO{D_{sd}} \nonumber\\ & \quad + \left( {{b_2} + {b_{2,d}}} \right)RelHumidit{y_{sd}} \nonumber\\ & \quad + \left( {{b_3} + {b_{3,d}}} \right)WindSpee{d_{sd}} + {b_4}PMRati{o_{sd}} \nonumber\\ & \quad + {b_5}Tem{p_{sd}} + {b_7}\% Cultivate{d_s} + {b_8}Fores{t_s} \nonumber\\ & \quad + {\varepsilon _{st}}\left( {{b_{0,d}}{b_{1,d}}{b_{2,d}}{b_{3,d}}} \right)\sim N\left[ {\left( {0,0,0} \right),{\mathbf{\Psi }}} \right] \end{align} \tag{ 1 }$$

where PM_2.5,sd represents ground-level PM_2.5 concentrations in µg m⁻³ at each site (s) for each day (d); b₀ and b_0,d are respectively the fixed and random intercepts for the model; AOD_sd denotes the retrieved MAIAC AOD values at site s and day d with fixed and random day-specific slopes (b₁ and b_1,d); RelHumidity_sd represents the measured relative humidity at site s and day d with fixed and random day-specific slopes (b₂ and b_2,d); WindSpeed_sd is the average wind speed at site s and day d with fixed and random day-specific slopes (b₃ and b_3,d); PMRatio_sd represents the ratio of PM_2.5 to PM₁₀ values; Temp_sd is the average daily temperature at each site; %Cultivated_s is the percentage cultivated to non-cultivated land at each site; Forest_s is an indicator variable denoting pixel forest cover; b_0,d, b_1,d, b_2,d, b_3,d are multi-variate normally distributed and Ψ represents the variance-covariance matrix for all random effects. The specific fixed effects for AOD, relative humidity, and wind speed aid in accounting for the average effects of these variables on the PM_2.5 concentrations and the random effects are included to account for the daily variability between both the dependent and independent variables. Other potential confounders were assessed (boundary layer height, emission point sources, heat index, etc.). However, these did not influence the results and were omitted in the final model.

The purpose of the first stage of the model is to estimate the temporal PM_2.5/AOD relationships with included effects of added covariates. However, we expect that the relationship will vary in space as well. To account for this potential additional variation, we added a second stage to the model utilizing geographically weighted regression (GWR) methods to create a continuous surface of estimates for parameters at each location. This incorporated using adaptive bandwidth selection methods in order to minimize the Akaike Information Criterion (AIC) value and aid in model selection. The GWR model structure can be expressed as:

$$\begin{align}P{M_{2.5}}.residual{s_{sd}} & = {\beta _{0,s}} + {\beta _{1,s}}Coast.distance \nonumber\\ & + {\beta _{2,s}}Mexicali.distance \nonumber\\ & \quad + {\beta _{3,s}}Elev + {\beta _{4,s}}Roads + {\varepsilon _s}\end{align} \tag{ 2 }$$

where PM_2.5.residuals represents the residual values from stage one of the model at site s for each day d; Coast.distance is the Euclidean distance calculated from each site s to the coast;, Mexicali.distance represents the distance from each site to Mexicali, Mexico; Elev is the elevation at each site s in meters;; and Roads is the sum of primary roads and highways within each grid cell; β_0,s, refers to the location-specific intercept with β_1,s, β_2,s, β_3,s, β_4,s representing location-specific slopes for each of the parameters. The results from the second stage were then used as a calibration measure for the measurements obtained from the first stage of the model.

In order to assess the goodness of fit for the final model, we compared the outcome of the fitted model with the observed values using cross-validated coefficients of determination (CV R²) and root mean squared error (RMSE). A 10-fold cross-validation of the model was done by first randomly splitting the dataset into 10 equal subsets. The model was then run 10 times—each time one subset was kept in reserve as a test sample while the other 9 subsets were used to train the model. Since the calculation of the PMRatio variable included interpolation of PM₁₀ observations, we chose to recalculate the PMRatio parameter for each of the 10 runs to avoid including information from the left-out subset. We then estimated predicted values for the remaining subset. Finally, the agreement between the predicted and observed values was then tested using the R² and RMSE values and a comparison was made between the cross-validated model and the observations in order to assess agreement and/or potential over-fitting. Additionally, we conducted sensitivity analyses by running the full 2-stage model leaving out key parameters (AOD and PMRatio). We conducted all modeling and analyses in R 3.6.0 (2019) and ESRI ArcGIS^® 10.6 (2018).

The annual mean PM_2.5 concentration for all monitoring sites was 10.7 µg m⁻³, with maximum values as high as 78.8 µg m⁻³ present during the study period. The coverage for AOD values was 62% and overall mean AOD value was 0.08 with a maximum value of 2.96 during the study period. Maximum wind speeds reached 18.30 m s⁻¹ with mean annual wind speeds measured at 4.78 m s⁻¹. Additional parameters and corresponding mean, standard deviation, and range of the statistics are included in table 1.

The model was fitted and the fixed effects estimated from the stage 1 linear mixed effects model are provided in table 2. The contributions of the intercept and parameters in the model are all significant at α = 0.05 level. The positive β values are indicative of a positive relationship between AOD, PM ratios, temperature, relative humidity, and cultivated land cover. Negative β values shown for wind speed and forest coverage suggest a negative relationship with PM_2.5 concentrations. β values represent the change in the variable (keeping all others constant) that would increase the PM_2.5 concentration by 1 µg m⁻³.

The first stage of the model has a CV R² of 0.77 (model fitting R² = 0.78) and CV RMSE of 2.41 µg m⁻³ (original model RMSE = 2.38 µg m⁻³). Compared with the original model, the results from the CV suggest slight model over-fitting as shown by the changes in statistical measures of R² and RMSE. The second stage of the model (or full model) incorporates the first stage and GWR methods and resulted in a CV R² increase to 0.80 (original model R² = 0.85) and a change in CV RMSE to 2.25 µg m⁻³ (original model RMSE = 1.96). As shown in these validation results, the second stage of the model improved the overall prediction performance with a CV R² increase of 0.03 from the first stage to the second stage of the model and a 0.16 µg m⁻³ decreased CV RMSE between the two stages. This improvement in accuracy could indicate the ability of the GWR methods for capturing more of the spatial variation in the data than is possible with the LME stage alone. These results suggest that adding a second stage to account for spatial as well as temporal variability can substantially improve model accuracy.

Annual mean PM_2.5 at ground stations alone and across the entire domain based on inverse distance weighting can be seen in figure 2. This comparison shows an increased coverage of exposure with this simple interpolation technique; however, this method lacks the spatial detail necessary to be confident in the estimated exposures located at greater distances from monitoring stations. For example, with the simple interpolation method, major differences in PM_2.5 concentrations can be found, such as those near the Mexican border. The annual mean PM_2.5 estimated from model fitting for the 1 km × 1 km grid is shown in figure 3, with figure 3(A) showing results of the analysis including the PMRatio parameter and figure 3(B) showing results without including PMRatio. As seen in figure 3(A), annual averaged PM_2.5 tends to result in high concentrations seen in population centers, along some major highways, and in valleys and canyons. Estimated concentrations align closely with the topography of the surrounding areas (see figure 1). In the Los Angeles area, lower concentrations of PM_2.5 are seen on the coast, with increasing concentrations to inland population centers (i.e. Riverside). In San Diego, the same trends occur, but with higher coastal PM_2.5 than Los Angeles. The Imperial Valley area is subject to higher concentrations due to transport from Mexico, major highways, high airborne dust, and surrounding topography.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During the process of model fitting, it was found that the presence of dust (or larger sizes of PM) in the troposphere is an integral part of the AOD-PM_2.5 relationship in southern California (β = 14.79, p ≤ 0.00). We accounted for this by adding a model parameter representing the observed PM_2.5 divided by the monthly mean PM₁₀ for each site. The pattern of PM_2.5 to PM₁₀ ratio is shown by the location of EPA monitoring stations in Supplemental figure 1. There was a strong positive correlation between PMRatio and observed PM_2.5 (r = 0.7) and a weak negative association with observed PM₁₀ (r = −0.4, see Supplemental table 1 for complete correlation results). The importance of this parameter was also tested by running the model with and without the PM_2.5 to PM₁₀ ratio. In figure 3(A), the ratio of PM_2.5 to PM₁₀ was not included in the model, resulting in full model CV R² was 0.70 with a CV RMSE of 2.81 µg m⁻³. However, with the inclusion of this ratio parameter in figure 3(B), accounting for larger particles significantly improved the model performance, resulting in the CV R² = 0.80 and a RMSE = 2.25. Detailed results comparing the full model with models without the PM_2.5 to PM₁₀ ratio and without AOD can be seen in Supplemental table 1, and temporal patterns of monthly AOD and the PM_2.5 to PM₁₀ ratio are shown in Supplemental figure 2. Notably, the spatial patterns of the PMRatio variable can change greatly from monitor to monitor—even between monitors relatively close in proximity to one another. This is likely due to the shorter airborne residence time of larger particles, which vary greatly depending on current conditions.

Seasonality of PM_2.5

Figures 4(A)–(D) illustrates the seasonal patterns of PM_2.5 concentrations in the study domain. The mean predicted concentrations varied by quarter with a mean of 4.24 µg m⁻³ in the first quarter, 6.77 µg m⁻³ in the second, 7.37 µg m⁻³ in the third, and 7.75 µg m⁻³ in the final quarter. Precipitation during California's rainy season (∼ October-April) can contribute to low PM_2.5 concentrations, as seen in the first and second quarter mean results. These means also follow the general patterns of PM_2.5 concentrations and temperature. Maximum concentrations were found in the final quarter (max = 22.77 µg m⁻³) and we see a decrease in first and second quarterly overall PM_2.5 compared to later levels. However, concentrations surrounding populated areas may still be pronounced.

Zoom In Zoom Out Reset image size