Satellite-Based Long-Term Spatiotemporal Patterns of Surface Ozone Concentrations in China: 2005-2019.

BACKGROUND: Although short-term ozone (O3) exposure has been associated with a series of adverse health outcomes, research on the health effects of chronic O3 exposure is still limited, especially in developing countries because of the lack of long-term exposure estimates.
OBJECTIVES: The present study aimed to estimate the spatiotemporal distribution of monthly mean daily maximum 8-h average O3 concentrations in China from 2005 to 2019 at a 0.05° spatial resolution.
METHODS: We developed a machine learning model with a satellite-derived boundary-layer O3 column, O3 precursors, meteorological conditions, land-use information, and proxies of anthropogenic emissions as predictors.
RESULTS: The random, spatial, and temporal cross-validation R2 of our model were 0.87, 0.86, and 0.76, respectively. Model-predicted spatial distribution of ground-level O3 concentrations showed significant differences across seasons. The highest summer peak of O3 occurred in the North China Plain, whereas southern regions were the most polluted in winter. Most large urban centers showed elevated O3 levels, but their surrounding suburban areas may have even higher O3 concentrations owing to nitrogen oxides titration. The annual trend of O3 concentrations fluctuated over 2005-2013, but a significant nationwide increase was observed afterward. DISCUSSION: The present model had enhanced performance in predicting ground-level O3 concentrations in China. This national data set of O3 concentrations would facilitate epidemiological studies to investigate the long-term health effect of O3 in China. Our results also highlight the importance of controlling O3 in China's next round of the Air Pollution Prevention and Control Action Plan. https://doi.org/10.1289/EHP9406.

Introduction

Numerous epidemiological studies have shown that short-term ozone () exposure is associated with a series of adverse health outcomes, including all-cause nonaccidental mortality (Yin et al. 2017) and respiratory morbidity (Barry et al. 2019). In 2015, pollution contributed to an estimated asthma emergency department visits globally (Anenberg et al. 2018). A handful of studies have also investigated the chronic health effects of , but their conclusions differed, as reviewed by Nuvolone et al.(2018). For example, Turner et al. (2016) reported a positive association between long-term exposure and all-cause mortality in a large prospective study in the United States. However, a meta-analysis reported that this association existed only in the warm season rather than the whole year (Atkinson et al. 2016). One potential reason for this inconsistency was the bias in the exposure matrices, especially in large long-term studies. Seltzer et al. (2018) highlighted the value of a dense monitoring network in measuring long-term exposure, but such a network is unavailable in most developing countries. As the world’s most populous nation, China began to establish its national air quality monitoring network in 2013. To date, this network covers most Chinese cities, but the rural and suburban areas remain largely unmonitored (Xiao et al. 2020). The insufficient spatiotemporal coverage of measurements presents a major hurdle to retrospective epidemiological studies in China, especially those established before the 2010s (Wang et al. 2017).

Satellite remote sensing has become a promising tool to extend the records of measurements in space and time. For example, the Ozone Monitoring Instrument (OMI) is a nadir-viewing ultraviolet (UV)-visible () solar backscatter spectrometer aboard the National Aeronautics and Space Administration (NASA) Aura satellite designed to measure the total column and other trace gases (Levelt et al. 2006). Since its launch in July 2004, researchers have been exploring ways to estimate boundary-layer levels using the OMI and other satellite data. Based on the OMI’s measurements, Liu et al. (2010a) developed an optimal estimation technique to retrieve the profile from the surface to to produce the Smithsonian Astrophysical Observatory (SAO) OMI Ozone Profile (OMPROFOZ) product. The OMPROFOZ product has been shown to capture enhancements of lower tropospheric over China (Shen et al. 2019) and East Asia (Hayashida et al. 2015), but its performance varies in space and by season owing to various factors affecting the horizontal and vertical distribution of tropospheric (Huang et al. 2017). In China, the daily correlation between the OMPROFOZ tropospheric column and ground-level measurements varies from in high-latitude regions to in low-latitude regions (Shen et al. 2019).

Surface-level is formed by complex photochemical reactions between volatile organic compounds (VOCs) and nitrogen oxides () in the presence of heat and solar radiation (Li et al. 2020; Pu et al. 2017). In addition to the abundance of precursors, the production of surface can be strongly influenced by meteorology. For example, high temperature boosts chemistry by increasing the decomposition rate of peroxyacytyl nitrate (PAN), thus preventing the sinkage of and peroxyl radicals (Fischer et al. 2014; Lu et al. 2019a). Water vapor may affect surface production by modulating the hydrogen oxide () radicals essential to production from oxygen (Lu et al. 2016). Besides these persistent effects, lightning would result in a surge in emission, significantly elevating short-term ground-level concentrations (DeCaria et al. 2005; Kang et al. 2020).

Although photochemical reactions predominantly determine tropospheric concentrations (Lelieveld and Dentener 2000; Monks 2005), in some circumstances the stratospheric -rich air may penetrate rapidly into the lower troposphere and cause a sharp increase in ground-level pollution (Knowland et al. 2017). The peak of this cross-tropopause mass transport, known as a stratospheric intrusion (SI), in the northern hemisphere usually occurs in springtime, especially in high-altitude regions, such as the Qinghai-Tibet plateau (Appenzeller et al. 1996; Itahashi et al. 2020; Lin et al. 2012; Lu et al. 2019b; Wang et al. 2020b). Changes in large-scale climate patterns, such as El Niño, are associated with increased SIs and surface concentrations (Shen and Mickley 2017; Xie et al. 2014). In addition, vegetation has complicated effects on surface . On one hand, plants can remove surface through dry deposition (Clifton et al. 2020); however, on the other hand, plants may emit VOCs to the atmosphere and affect pollution (Kigathi et al. 2019). Therefore, land cover types are also important in determining ambient concentrations.

In this study, we developed a national-scale machine learning model to estimate historical ambient ground-level concentrations in China at a monthly level from 2005 to 2019 at a 0.05° spatial resolution. In addition to the OMPROFOZ profile, we included meteorological factors, land-use information, precursors, and indicators of anthropogenic emissions to account for the complicated formation and removal processes of surface . We first present our model development strategy, then evaluate model performance using statistical techniques, as well as ground measurements not included in model training. Finally, we investigate the spatiotemporal trend of and discuss the drivers of these patterns.

Data and Methods

Model Development

We used a random forest framework to estimate monthly mean daily maximum 8-h average (MDA8) concentrations in China from 2005 to 2019. The overall study workflow is illustrated in Figure S1. Briefly, we first extracted the fraction of the boundary-layer column from the OMPROFOZ profile. The depth of this layer for each grid cell was determined dynamically by the tropopause pressure (Liu et al. 2010a). Missing fraction values were imputed with daily random forest models incorporating the Modern-Era Retrospective Analysis for Research and Applications (MERRA-2) meteorological fields and surface flux. The gap-filled boundary-layer fractions were then used to calculate the full-coverage boundary-layer column (in Dobson units). The surface-level monthly average MDA8 concentration was generated by a separate random forest model trained with the gap-filled boundary-layer column, meteorological fields, land-use terms, elevation, and population density. The details of these predictors are described below.

Ground Measurements

Ground-level monitoring data from 2013 to 2019 were obtained from the China National Environmental Monitoring Center (CNEMC; http://www.cnemc.cn/). China’s national air quality monitoring network measures hourly concentrations with either UV absorption (for point analyzers) or differential optical absorption spectroscopy (for open-path analyzers). MDA8 concentration was defined as the maximum 8-h moving average concentration (containing at least six valid hourly values) within a calendar day. We chose to use MDA8 rather than MDA1 or 24-h average because MDA8 has been widely used in health effects research (Lu et al. 2020; Lyu et al. 2019).

The quality assurance of the CNEMC data was conducted primarily based on its official standard (CNEMC 2012). Briefly, a valid 8-h moving average concentration must contain at least six hourly measurements. No regulatory standard has yet been set up for MDA8 , but 20 hourly measurements were generally required for daily air pollutant concentrations. Therefore, we also removed all the observations during a given day from the stations with hourly measurements to balance data abundance and data quality. A total of 2,443 (3%) monthly MDA8 concentrations were removed from the data set.

We further separated rural stations from urban stations to test our model performance in different settings. In compliance with China’s official regulations (NBS 2008), we defined areas with a population density of ( per grid cell) as rural. Consequently, 406 monitoring stations of 1,532 (27%) were identified as rural stations, contributing to a total of 24,405 (31%) data points. Note that population density is not the only determinant of rural/urban status listed in China’s official standards. We used population to identify rural/urban stations primarily because other economic and political determinants were less quantifiable in this modeling study.

monitoring data from before 2014 was obtained from the Tropospheric Ozone Assessment Report (TOAR) (Xu et al. 2020). This data set contains MDA8 concentrations from eight different stations, including Mt. Waliguan (WLG, located at 36.30°N, 100.9°E), Shangdianzi (SDZ, at 40.39°N, 117.00°E), Lin’an (LAN, at 30.30°N, 119.73°E), Longfengshan (LFS, at 44.73°N, 127.60°E), Xianggelila (XGLL, at 28.01°N, 99.68°E), Akedala (AKDL, at 47.10°N, 87.93°E), Gucheng (GCH, at 39.13°N, 115.67°E), and China Meteorological Administration (CMA, at 39.95°N, 116.32°E). The locations of the TOAR sites can be found in Figure S2. The TOAR historical monitoring data was not used to train the final model but, rather, served as an independent external validation data set.

Satellite Remote Sensing Data

In this study, we used the OMI OMPROFOZ product developed at the Harvard SAO, publicly available at the NASA Aura Validation Data Center (AVDC). The profile is retrieved at 24 layers ( per layer) from the surface to from the spectral range , using the optimal estimation approach. It is based on the initial retrieval algorithm of Liu et al. (2010a) with modifications described by Kim et al. (2013). The layers between the surface and the tropopause were defined on a daily basis. That is to say, the layers’ pressure boundaries were initially set at atm for , and . For each individual day, the National Centers for Environmental Prediction (NCEP) tropopause pressure were used to replace the pressure level closest to it. The layers between the surface and the tropopause were then reassigned based on equal logarithmic pressure intervals (Liu et al. 2010a). OMPROFOZ’s retrieval errors due to precision (instrument random noise) and smoothing errors (insufficient vertical resolution) ranged from 1–6% in the stratosphere to 6–35% in the troposphere. The retrieval was performed at a spatial resolution of at nadir and gridded to a 0.5° resolution for ease of use.

For more intuitive interpretations, we defined OMPROFOZ’s boundary layer (from surface pressure to hPa) as Layer 24 (L24). Correspondingly, L23–L1 represents the second-lowest layer to the top layer. Our model initially used the boundary-layer (L24) partial column from the OMPROFOZ profile. To better understand OMPROFOZ’s role in modeling long-term pollution, we also tested whether using additional tropospheric columns (i.e., the summation of L22–L24) or replacing the retrievals with an a priori profile (both L24 and the summation of L22–L24) would influence the model performance.

Tropospheric column nitrogen dioxide () concentration from 2005 to 2019 was extracted from the OMI global product named OMNO2d (Krotkov et al. 2017). It provides global daily tropospheric column concentration at a 0.25° resolution. We extracted the tropospheric column concentration with cloud screening (i.e., pixels with a cloud fraction were removed for quality assurance) as a measurement of precursor in the present study.

MERRA-2 Assimilated Data

We used the MERRA-2 meteorological data in the present study. The MERRA-2 provides the latest NASA atmospheric reanalysis data starting from 1980. It has a native resolution of and 72 vertical layers (Gelaro et al. 2017). We extracted surface-level meteorological fields, as well as those between surface and 150 hPa, to account for the effects of stratospheric intrusion (Knowland et al. 2017). The detailed list of MERRA-2 meteorological and chemistry fields we used is provided in Table S1.

Lightning Flash Density

Previous studies have shown that lightning flash is an important enhancer of tropospheric because of its considerable contribution to (DeCaria et al. 2005; Finney et al. 2016), especially in springtime (Lu et al. 2019b). We obtained global monthly lightning flash density data from the Harvard-NASA Emissions Component at a resolution, adopting an optimal regional scaling algorithm to reduce the bias of satellite-driven tropical lightning data (Murray et al. 2012).

South and Southeast Asia Wildfire

The massive human-initiated biomass burning in South and Southeast Asia greatly enhances springtime pollution over China, especially in the southwest (Ni et al. 2018; Wang et al. 2011). To test whether incorporating foreign wildfire data would increase model performance, we obtained the Moderate Resolution Imaging Spectrometer (MODIS) daily active fire data from 2005 to 2019 via the Fire Information for Resource Management System (FIRMS) data achieve (NASA Earthdata 2019). Fire points from 10 South and Southeast Asian countries—namely, Bangladesh, Bhutan, Cambodia, India, Laos, Myanmar, Nepal, Pakistan, Thailand, and Vietnam—were selected to capture most fire points in China’s southern neighbors. Fire radiative power (FRP) was used as a quantitative proxy of fire-related emissions (Li et al. 2019a; Wooster et al. 2003).

Land Use, Population, Road Length, and Digital Elevation

Annual land cover maps were obtained from the European Space Agency (ESA) Climate Change Initiative (CCI) (ESA 2017) for 2005–2015 and the Copernicus Climate Change Service (C3S) Climate Data Store (CDS) for 2016–2019 (CDS 2021). The C3S land-use product used the same methodology as the ESA CCI land cover maps to guarantee long-term continuity, according to the product manual (C3S 2021). Both products provide 23 types of land cover at a spatial resolution of .

LandScan global population data were obtained from the Oak Ridge National Laboratory (https://landscan.ornl.gov/). This data set provides annual population density from 2005 to 2019 at a resolution. Road networks were obtained from the Global Roads Open Access Data Set (CIESIN and ITOS 2013). This data set was compiled from sources before 2010 (specific date unavailable). Total road length was included in our model as a proxy of traffic emissions. We also used elevation data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM), version 3 (NASA and METI 2019). This latest version of ASTER GDEM has increased accuracy compared with previous versions.

Missing Value Imputation

Equation 1 Illustrates the relationship between the total column and the boundary-layer (L24) column . where denotes the boundary-layer column (in Dobson units); is the boundary-layer fraction of total column ; and is the total column . The OMPROFOZ product has a nonnegligible portion of missing values that will reduce the spatial coverage of and affect predicted ground-level concentrations. We first filled the data gap of OMI , then multiplied it with MERRA-2 to get the final . We chose this approach over directly imputing the OMPROFOZ product because OMI total column and Microwave Limb Sounder (MLS) stratospheric profiles have been assimilated into the MERRA-2 data after 2004 (Wargan et al. 2017). The correlation coefficient between MERRA-2 column and OMI column concentration was 0.97 in the present study. The missing values of OMI were imputed with random forest models that incorporated MERRA-2 meteorological fields and surface flux measurements. These models were trained on a daily basis at the native resolution of the OMPROFOZ product.

Other OMI-derived partial column concentrations were processed similarly. To specify, missing values in other retrieved partial column concentrations were imputed independently with the same process to the boundary layer but with different meteorological fields corresponding to their pressure levels (i.e., 500–700 hPa for L23, 350–500 hPa for L22). The a priori profile (L22, L23, and L24) was processed independently but with the same methodology as the retrieved profile.

Data Integration

We created a 0.05° resolution modeling grid across China for data integration and model construction (Figure 1). A buffer region was added to China’s national boundaries to ensure data sensitivity at the border area. The total number of grid cells was 399,513. Three megacity clusters were selected to study the regional patterns of pollution—namely, the North China Plain (NCP), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD). We used an inverse distance weighting method to resample data at a spatial resolution coarser than 0.05°, including all OMI-derived column concentrations, tropospheric column concentrations, MERRA-2 meteorological fields, and lightning flash density. In addition, we calculated the percentage coverage of different land-use categories, average elevation, total road length, and total population for each pixel. Daily data were then aggregated to the monthly level. Population and land use were processed at the annual level, whereas elevation and road length were fixed during the entire study period. After data integration, we selected grid cells that contained air quality monitoring stations to generate the training data set. All valid ground-level MDA8 observations within each grid cell were averaged on a monthly basis to match other model parameters.

Figure 1.

The study domain and three major city clusters. The study domain covered China plus a buffer region that extends outside the national boundary. The points are the location of the China National Environmental Monitoring Center monitoring sites. The color scale represents elevation (km). Note: NCP, North China Plain; PRD, Pearl River Delta region; YRD, Yangtze River Delta region.

The foreign wildfire data was processed in a different way because our study domain did not extend to China’s neighboring countries other than the buffer area. We assumed that the foreign fire points have an additive and distance-dependent influence on China’s concentrations. Therefore, the contributions of the South and Southeast Asian countries were quantified with the equation below (Equation 2). where denotes the influence of foreign fire emissions on grid cell at day ; represents the FRP for fire point at day ; and is the distance between fire point and grid cell .

Model Training, Validation, and Parameter Comparison

We divided the study period into the training period (2014–2019) and the hindcast period (2005–2013). The year 2013 was excluded from model training owing to fewer numbers of monitors and unstable data quality at the onset of the Chinese national monitoring network. We trained two separate random forest models with the same set of predictors for springtime (March–April–May) and the rest of the year because of the significantly different pattern of springtime in the northern hemisphere (Lin et al. 2012; Ni et al. 2018). To specify, the aforementioned rest of the year includes summer (June–July–August), fall (September–October–November), and winter (December–January–February). The detailed list of predictors included in our random forest models can be found in Table S1.

Note that some parameters were not included in the initial model. To specify, our original model was established with the boundary-layer (L24) column concentrations from the OMPROFOZ retrieved profile. The OMI L22–L24 retrieval and the a priori partial columns were used separately to compare model performance with different OMI-derived fields. In addition, the MODIS FRP was not initially used because it was generated completely outside of the study domain. Therefore, the original model hereafter represents the model with the retrieved profile L24 and without foreign fire emissions.

All models were validated with 10-fold random cross-validation (CV), that is, we randomly divided the original data set into 10 equal-sized subsets, used 9 of them to train a model, and made predictions on the left-out subset. This process was repeated 10 times so that each monthly mean MDA8 measurement would have a corresponding predicted value. Model performance metrics including and root-mean-square error (RMSE) were calculated using the measurement–prediction pairs. Similarly, we conducted a 10-fold spatial CV to test whether our model could make reliable predictions at locations without ground monitors. In the spatial CV, the original data set was randomly divided on the basis of each data point’s location. Model predictions at a given location were generated by a model trained with data elsewhere. Finally, a temporal CV was conducted in which models trained for a given year would be validated with data from other years in the training period to test the reliability of model predictions in the hindcast period.

Given that each variable provides a different contribution to the overall model performance, the predictors’ importance was evaluated with a permutation method (Altmann et al. 2010). Briefly, a variable’s importance represents the percentage increase in the model’s total mean squared error if this variable is replaced by its random permutation.

We used R (version 3.6.3; R Development Core Team) to process data and perform statistical analyses. Package ranger was used for training the random forest models.

Coupled Trend Evaluation between and fine particulate matter

To reduce severe air pollution, the State Council of China enacted the Air Pollution Prevention and Control Action Plan (APPCAP) in 2013. This policy has resulted in a substantial reduction in ambient fine particulate matter [PM in aerodynamic diameter ()], but the contemporary concentrations increased unexpectedly (Huang et al. 2018; Wang et al. 2020a). We obtained model-estimated nationwide concentrations at a resolution (Liang et al. 2020; Xiao et al. 2021) to evaluate the nationwide long-term coupled change between and under the APPCAP. The correlation between and was examined with a partial correlation analysis, controlling for temperature, relative humidity, and total precipitation from MERRA-2.

Results

Model Performance and Parameter Comparison

The performance of the original model is illustrated in Figure 2. The predicted monthly average MDA8 concentrations from the combined springtime and non-spring model were in good agreement with the ground-based observations, with a random CV of 0.87 and an RMSE of (Figure 2A). The spatial CV had an almost identical performance with the random CV , (Figure 2B). The temporal CV had a slightly lower of 0.76 and a higher RMSE of (Figure 2C). The regression lines between the predicted and observed concentrations were close to the 1:1 line for all three types of CV.

Figure 2.

The integrated model performance of the original model with OMI L24 retrieval. (A) Random CV; (B) spatial CV; and (C) temporal CV. The functions on the bottom-right corners are the regression functions between the predicted and observed monthly mean MDA8 ozone concentrations. Red dashed lines represent the regression line between the predictions and observations; black solid lines represent the . The color scale represents the density of the points. Note: CV, cross-validation; L24, boundary layer; MDA8, daily maximum 8-h average; OMI, Ozone Monitoring Instrument.

As shown in Figure S3, replacing OMI L24 retrieval with the summation of L22–L24 retrievals would not substantially impact the model performance. The and RMSE were almost identical with the original model for all three CV types. Similarly, using the a priori profile (L24 or L22–L24) also had a minimum influence on the overall performance. Although using a priori L22–L24 would increase the temporal CV by 0.6%, such an attempt also had a negative impact on the predicted spatial distribution of . Figure S4A shows an example of the spatial artifact (horizontal gap in the predicted concentrations) from the model with a priori profile L22–L24. This kind of artifact was not seen in the original model (Figure S4B). Given that using the alternative partial column amount would not improve the model performance, all results presented hereafter used the model with OMI L24 retrieval unless otherwise specified.

The season-specific model performance showed that our model had a lower performance in spring than in other seasons (Figure S5). To specify, the springtime were 0.72, 0.71, and 0.53 for random, spatial, and temporal CV, respectively. In addition, the RMSE for all types of CV in spring was around higher than in fall, but the overall concentrations were comparable in these two seasons.

In addition, our model performed better in urban regions than in rural regions (Figure S6). The random, spatial, and temporal CV were 0.88, 0.87, 0.76 in urban regions and 0.83, 0.81, 0.72 in rural regions. The prediction errors (i.e., RMSE) were also higher in rural regions than urban areas ( for all CV types).

Our predicted MDA8 concentrations also agreed well with the TOAR historical data monitoring data before 2014 (overall , ), except for the XGLL station. Site-specific time-series comparison (Figure S7) showed that the predicted trends were mostly identical with the observations at stations CMA, GCH, and LFS. Although our model may have underestimated concentrations at the LAN, SDZ, and WLG stations, especially in springtime, it still captured most of the ’s temporal variation over these locations. The worst agreement was observed at the XGLL station, where our model had almost no sensitivity to the springtime peak concentrations. The monitoring data was mostly incomplete at the AKDL station, but it may also indicate some springtime underestimations in that region. Incorporating the MODIS FRP data would neither significantly improve the springtime nor the overall agreement with the TOAR historical data (Figure S8). Therefore, we excluded MODIS FRP from the model given that introducing data fields external to the study domain may incur unexpected uncertainty.

Although the OMI L24 retrieval had a low raw correlation with the CNEMC data (, as presented in Figure S9), removing this parameter would result in a slight decrease in the model performance compared with the original model ( decrease in for all three types of CV, as shown in Figure 2; Figure S10). Besides, the agreement between model predictions and the TOAR historical monitoring data would also be worse if OMI L24 retrieval is removed from the model (overall dropped from 0.73 to 0.70, RMSE increased from 20.68 to ) although the predicted temporal trends were generally similar (Figures S7 and S11).

Predictor Importance Ranking

Most of the top predictors in both the spring and non-spring models were meteorological factors, but their relative orders differed (Figure S12). Non-meteorological variables were more important in the non-spring model. For example, the five most important non-meteorological variables in the spring-excluded model were OMI , population, the proportion of rainfed cropland, irrigated or post-flooding cropland, and elevation. However, their importance ranking all dropped in the spring model. The gap-filled OMI boundary-layer was one of the most important variables in the spring model, and those that gained the most importance were associated with stratospheric intrusion (e.g., tropopause pressure and vertical wind columns) and lightning flash activity.

Seasonality and Spatial Heterogeneity of Levels in China

As shown in Figure 3, the spatial distribution and severity of surface pollution in China varied by season. The mean springtime MDA8 concentrations over 2005–2019 were mostly . Moderate pollution () was observed in Central-East China, especially around the boundary region between the NCP and the YRD. Southwest China, including the Sichuan Basin and Yunnan Province, was also moderately polluted in this season (Figure 3A). In summer, heavy pollution was widespread in China, except for the southwest. The NCP had the worst pollution, with the highest 15-y mean MDA8 concentrations around Beijing approaching the national level-2 air quality standard () (Figure 3B). In fall, the PRD region became the most polluted area after a sharp decrease in concentrations in North China (Figure 3C). Winter has the lowest levels, especially around the NCP region, and only the low-latitude regions may have an concentration approaching (Figure 3D).

Figure 3.

Spatial distribution of seasonal average MDA8 concentrations during 2005–2019. The concentrations were predicted by the original model with OMI L24 retrieval and averaged seasonally over 2005–2019. (A) Spring: March–April–May; (B) summer: June–July–August; (C) fall: September–October–November; and (D) winter: December–January–February. Note: MDA8, daily maximum 8-h average; , ozone; OMI, Ozone Monitoring Instrument.

Different regions appeared to have different seasonal patterns (Figure 4D). For example, concentrations in the NCP (Figure 4B) and nearby northern regions started to rise in spring and peaked in summer. After a sharp decrease in fall, these regions would have the lowest levels in the winter. In contrast, high concentrations could persist from spring to fall in southern China, without a clear peak in summer. Sporadic hot spots may even be found in winter (Figure 3D). The concentration in the YRD (Figure 4C) was a mixture of the two patterns above, that is, it had a flatter summer high and a winter low, but high concentrations also persisted in spring and fall. Most regions in China followed the aforementioned patterns, with a few exceptions. For example, Yunnan province saw levels peak in spring and then drop substantially in summer. concentrations on the Qinghai-Tibet plateau were relatively stable throughout the year, with only a mild peak occurring in summer.

Figure 4.

Region-specific seasonal mean MDA8 concentrations during 2005–2019 for (A) China as a whole, (B) the North China Plain (NCP), (C) the Yangtze River Delta region (YRD), and (D) the Pearl River Delta (PRD). The concentrations were predicted by the original model with OMI L24 retrievals. The bars and the corresponding numbers represent the season-specific regional mean MDA8 concentrations () over 2005–2019. Note: Fall, September–October–November; MDA8, daily maximum 8-h average; , ozone; OMI, Ozone Monitoring Instrument; spring, March–April–May; summer, June–July–August; winter, December–January–February.

’s spatial heterogeneity was observed not only nationwide but also at the city level. Figure 5A shows the model-estimated concentrations in the YRD region in August 2019, whereas Figure 5B shows population densities. Some population centers, such as the subregions labeled A (Bengbu and Huainan City) and B (Nanjing Metropolitan Area), had lower levels than their surrounding areas. On the contrary, the subregions labeled C (Anqing City) and D (Quzhou and Jinhua City) were more polluted than their surrounding areas. As shown in Figure S13, the tropospheric column concentrations were higher in subregions A and B than in C and D.

Figure 5.

Summer peak and population in the Yangtze River Delta (August 2019). (A) The model-predicted monthly mean MDA8 concentrations in August 2019; and (B) the 1-km population data for 2019 from LandScan. The boxes represent some YRD cities and their surrounding areas: A, Bengbu and Huainan City; B, Nanjing Metropolitan Area; C, Anqing City; and D, Quzhou and Jinhua City. Note: MDA8, daily maximum 8-h average; , ozone; YRD, the Yangtze River Delta.

The Long-Term Trend of in China

Figure 6 shows the long-term trend of concentrations in the season (defined as March–November) in China. Before 2014, mean concentrations during the season fluctuated from year to year but generally stayed at the same level nationally, at . levels in the NCP, YRD, and PRD regions were higher than the national average but were also stable around their respective long-term averages. A sharp decrease was observed in the YRD and PRD from 2014 to 2016. After 2016, levels nationwide started to rise at various paces. For example, seasonal mean concentrations were almost identical for the YRD and NCP before 2014, but the former experienced a sharper increase from 2015 to 2019 owing to the significant drop from 2014 to 2015. After 2018, the mean MDA8 concentrations in the season exceeded in all these regions.

Figure 6.

Annual trend of mean season average MDA8 concentrations from 2005 to 2019. season is defined as March–November. The lines with different colors and marks represent the model-predicted mean MDA8 concentrations in different regions. Gray line with round marks denotes China; the red line with cross marks, the PRD; the orange line with triangular marks, the NCP; and the blue line with square marks, the YRD. Note: MDA8, daily maximum 8-h average; NCP, the North China Plain; , ozone; PRD, the Pearl River Delta; YRD, the Yangtze River Delta.

Summertime (June–July–August) mean MDA8 concentrations showed an overall increasing trend during 2005–2019 in the whole of China (, ), the NCP (, ), and the YRD (, ) (Figure S14). However, similar to the season averages, no significant trend was observed for summertime mean MDA8 concentrations over 2005–2013 (, ; , ; , ; and , for the whole of China, the NCP, the YRD, and the PRD, respectively) (Figure S15). The NCP was the only region that had an overall increasing springtime pollution over 2005–2019 (, ) (Figure S16), whereas no significant trend was found for 2005–2013 across China (, ; , ; , ; and , for the whole of China, the NCP, the YRD, and the PRD, respectively) (Figure S17). As shown in Table S2, the TOAR monitoring data also did not exhibit a significant increasing trend in pollution during 2005–2013 except for the summertime pollution at the SDZ station (, ).

The predicted temporal trend was similar to the CNEMC observations over the grid cells with monitoring sites (Figure S18). Although our model may slightly underestimate the overall season MDA8 concentrations in China, the predicted summer peak concentration was almost identical with the CNEMC observations in China as a whole, as well as with the NCP, YRD, and PRD (Figure S19). Unlike our model predictions, the 2018–2019 increase in the -season mean MDA8 concentrations was not seen from the CNEMC data except for the PRD. However, there was a substantial increase in some subregions, including Central-East China (Figure 7A) and the Shandong Peninsula (locations listed on Figure S20; Figure 7B), as observed by both the CNEMC data and our predictions.

Figure 7.

The comparison of the summertime (June–July–August) mean MDA8 concentrations between our model predictions and the CNEMC monitoring data over 2014–2019 for selected regions. (A) Central-East China; and (B) the Shandong Peninsula. The blue columns represent the CNEMC observations, and the orange columns represent our model predictions. The height of the columns and the error bars represent the mean MDA8 concentrations and the standard errors. Note: CNEMC, China National Environmental Monitoring Center; MDA8, daily maximum 8-h average; NCP, the North China Plain; , ozone; PRD, the Pearl River Delta, YRD: the Yangtze River Delta.

Couple Trend between and in China

The temporal trends of population-weighted and concentrations in China and three major regions are shown in Figure 8A–D. Despite the visible inverse correlation between and levels, partial correlation analyses controlling for temperature, relative humidity, and precipitation indicated that their association varied by region. A statistically significant negative correlation was observed in the YRD region (, ), whereas in the PRD region, this correlation turned positive (, ). No significant association was observed in China as a whole (, ) and in the NCP (, ) after we controlled for temperature, relative humidity, and precipitation.

Figure 8.

Time-series comparison for population-weighted ambient and concentrations in China over 2005–2018. (A) China as a whole; (B) the NCP; (C) the YRD; and (D) the PRD. r denotes the partial correlation coefficient controlling for temperature, relative humidity, and total precipitation from MERRA-2. -Values were for the r statistics. Blue lines represent the population-weighted monthly mean MDA8 predicted by our model, and the orange lines represent the population-weighted monthly mean concentrations from Xiao et al. (2021). Note: MDA8, daily maximum 8-h average; MERRA-2, Modern-Era Retrospective Analysis for Research and Applications; NCP, the North China Plain; , ozone; PRD, the Pearl River Delta; YRD, the Yangtze River Delta.

Discussion

In the present study, we trained a random forest model to predict long-term ground-level MDA8 concentrations in China. To our best knowledge, our model’s performance was the highest among similar studies in China. For example, the 0.1° model established by Liu et al. (2020) had a spatial CV of 0.68 and a temporal CV of 0.69 for monthly mean MDA8 concentrations. Compared with their study, our model had a finer spatial resolution (0.05°) and higher spatial and temporal CV values (0.86 and 0.76, respectively). Despite a coarser spatiotemporal resolution, our model generally had a comparable performance with those established in regions with sufficient historical monitoring data. For instance, the for daily MDA8 by Di et al. (2017) ranged from 0.7 to 0.8 across the United States over 2000–2012, and the land-use regression model of Adam-Poupart et al. (2014) had an of 0.65 for daytime 8-h average concentrations in Quebec, Canada. One reason for the improved performance was the inclusion of the OMPROFOZ profile instead of the total column. The latter is a noisy proxy for surface-level because approximately 90% of the atmospheric exists in the stratosphere (Fishman and Larsen 1987). Although the gap-filled boundary-layer column was not among the top predictors for the non-spring model, it was an important predictor in springtime when ground-level pollution is greatly influenced by SI and foreign transport (Figure S12). Another reason for the higher performance of our model is that the two primary sources of ground-level (i.e., photochemical reactions and SI) were both accounted for in our study, whereas previous studies typically focused on the impact of photochemical reactions on patterns.

With the improved spatial resolution, we were able to observe the complex relationship between human activities and concentrations, as presented in Figure 5. These regions have been reported to have a high /VOC ratio so that the abundance of VOCs controls the ambient concentrations (Liu et al. 2010b). Under the VOC-limited regime, the excessive concentrations would quench molecules through titration (Jhun et al. 2015). On the other hand, concentration in regions C (Anqing City) and D (Quzhou and Jinhua City), with moderate concentrations (Figure S13), fall under the -limited regime, with a positive association between and concentrations.

Our modeling results indicated that the seasonality of concentrations in China varied across different regions and was more distinct in northern China than in the south (Figures 3 and 4). This phenomenon could be explained by three reasons. First, the absolute and relative abundance of precursors have a major impact on concentration. The observed regional and national hot spots occurred predominantly in the megacity clusters (NCP, YRD, PRD, and Sichuan Basin), which can be partially attributed to high anthropogenic emissions of and VOCs (Li et al. 2019b). Second, meteorological conditions, such as temperature and solar radiation, significantly affect surface formation (Coates et al. 2016; Schnell et al. 2009), making summer the most polluted season in most regions. However, the active monsoon activities in summer, especially in the PRD, may increase cloud cover and weaken solar radiation. formation would be restrained under such a condition (Qu et al. 2021). The heavy rainfall during the monsoon season may also lead to a reduction in surface pollution. This may help explain the absence of summer peaks in the low-latitude areas, including the PRD. The stronger seasonal variation of solar radiation in high-latitude regions may have contributed to a shorter, but more distinct, season in the NCP. Finally, the intensity of the SI events varies by region. Deep SI events occur more frequently in the high-altitude regions, with some even reaching the ground (Lin et al. 2012, 2016). Consequently, high concentrations could be observed on the Qinghai-Tibetan plateau in spring despite lower levels of precursors and relatively low temperatures.

Although our model generally captured the seasonal variation of in China, we still found a substantial underestimation of high springtime concentrations, especially at the XGLL station (Figure S7). A potential reason for this underestimation is that foreign transport was not well accounted for in our model. According to Ni et al. (2018) foreign regions contribute 40–60% to China’s springtime below the height of . This enhancement is very prominent in Southwest China owing to the massive biomass burning in South and Southeast Asia (Wang et al. 2011), which peaked in spring (Yin 2020). However, incorporating wildfire emissions from South and Southeast Asia by assuming a distance-dependent influence did not address this underestimation (Figure S8). According to Wang et al. (2011), the prevailing westerly wind and active cyclonic activity in spring facilitate transport from South Asia to China. Therefore, foreign wildfire’s impact on China’s pollution is not determined by distance itself but, rather, follows a certain trajectory. We are unable to address this trajectory in the current model because the grid cells did not extend outside China’s boundary except for the buffer.

China has experienced rapid economic growth over the past decades. Nevertheless, the contemporary increase in anthropogenic emissions also greatly exacerbated air pollution over the nation. According to Lin et al. (2017), the 95th percentile summer (June–July–August) MDA8 increased by in China over 1995–2014 and the springtime (March–April–May) median increased by concurrently. We also observed an overall increasing trend in summertime pollution during 2005–2019 in the whole of China, the NCP, and the YRD, as well as an increase in springtime pollution in the NCP (Figures S14 and S16). However, the increasing trend was not constant from year to year. As seen from our model predictions, neither the springtime nor the summer surface concentrations in China exhibited a significant increasing trend over 2005–2013 (Figures S15 and S17). These findings were mostly consistent with the TOAR monitoring data (Table S2). SDZ was the only TOAR station that observed an increasing trend in summertime pollution over 2005–2013. The increasing trend at SDZ was primarily driven by an increase from 2005–2007, whereas no significant trend was observed for 2008–2013 (, ). Xu et al. (2020) also showed that both the annual highest MDA8 and the annual fourth highest MDA8 exhibited no significant long-term increasing trend for most TOAR stations, except for SDZ. Although historical monitoring data were lacking back to the early 2000s, it can be inferred that the overall increasing trend of China’s pollution during 1995–2014 was primarily driven by an increase before 2005. Tang et al. (2009) also reported that the concentrations in Beijing increased at a rate of during 2001–2006. The differential trends of pollution between 1995–2004 and 2005–2013 were possibly attributable to meteorological conditions. According to Sun and Wang (2017), the surface temperature in Northern and Northeastern China increased during the 1990s but stabilized during 2005–2014.

A sharp decrease in concentration was observed in China during 2014–2016, especially in the PRD region. This was not an isolated event in China given that the U.S. Environmental Protection Agency also reported the second-lowest concentrations on record during the 2014–2016 average period (U.S. EPA 2017). The significant reduction in concentration in the northern hemisphere was attributable to El Niño during this period. Olsen et al. (2016) reported that the anomalous cyclonic circulation induced by El Niño events coincided with decreased tropospheric . Shen and Mickley (2017) found a negative association between the El Niño–Southern Oscillation (ENSO) and summertime air quality in the south-central states of the United States. Although the effect of the ENSO on pollution in China has yet to be fully understood, indirect evidence supported the strong impact of climate factors between 2014 and 2016. Yang et al. (2019) reported that meteorological conditions contributed to a decrease of surface in 2016 and a decrease in 2014–2015 in the PRD region.

levels in China increased rapidly after 2016 (Figure 6). This may have to do with both China’s emission control policies and meteorological conditions. The issuance of APPCAP in 2013 resulted in a dramatic nationwide decrease in emissions and levels near the end of this 5-y plan (Zheng et al. 2018). However, formation falls under the VOC-limited regime in most Chinese urban centers. Reducing emissions became an enhancer of formation in populous areas, which were then transported to other regions (Liu and Wang 2020). In addition, lower levels would modify pollution because the photolysis rates were less attenuated by aerosol light scattering and absorption (Liu and Wang 2020). Decreased concentrations could also slow down the sink of hydroperoxyl () radicals, resulting in enhanced production (Li et al. 2019b). Such complex interactions among formation, levels, and meteorological conditions are reflected in the spatially varying associations between and concentrations revealed by our partial correlation analysis. The mild weather and year-round intensive human activities cause a long season (March–November) in the YRD. During this period, levels can vary from 25 to , allowing it to affect production negatively. The PRD region was the cleanest affluent city cluster in China in terms of both and . The effect of on production was less notable, and and concentrations were both positively correlated to the emissions of their precursors (Liu and Wang 2020). Although levels in the NCP are high in winter, its cold winter strongly suppresses formation, and no significant correlation was found between and levels after controlling for meteorological conditions. The 2018–2019 increase in surface pollution in China was more likely to be driven by climatological factors, particularly the increased foehn wind frequency and the subsequent changes in temperature and relative humidity (Li et al. 2020). That explained why this round of increase was only seen in several regions by the CNEMC monitoring data (Figure 7). Again, although pollution in China is not consistently worsening from year to year as a result of the modulation of climatological factors, it generally exhibited a prominent long-term increasing trend (Lin et al. 2017; Xu et al. 2020). Controlling pollution in China remains a challenging task that requires a better understanding of this air pollutant.

Our study has a few limitations. First, the road network data was fixed over the study period. This may result in an underestimation of China’s road density, especially after 2010. Although there are alternative road length data, we prefer not to mix multiple data sets with different methodologies to avoid introducing systematic errors. Second, the precursors included in our model were limited to . Satellite retrievals of VOCs are limited and often have weak signal-to-noise ratios in the boundary layer (Zhu et al. 2020). Although our current model was able to identify regions under the VOC-limited regime of production, we will continue to explore effective indicators of ground-level VOC to improve our exposure model. Finally, the present study focused on domestic determinants of ambient . The present model may underestimate pollution attributable to foreign transport. Future users of our data set should be cautious with the springtime underestimations, especially in Southwest China. The roles of long-range transport from outside China warrants further investigation.

Conclusions

We used a data-driven modeling framework to estimate long-term, high-resolution concentrations in China. Predictors that capture the influence of photochemical reactions and SI were included in our model. This model produced reliable historical monthly mean MDA8 concentrations for China at a 0.05° resolution with little bias. This 15-y long, full-coverage national data set of ambient concentrations includes the 9 y before China’s regulatory air quality monitoring network existed. It could accelerate research on the long-term health effects in China by enabling the use of large general population cohorts established in the 2000s, such as the Chinese Longitudinal Healthy Longevity Survey (Kuang et al. 2020) and the China Health and Retirement Longitudinal Study (Zhao et al. 2014).

Click here for additional data file.

Click here for additional data file.