The Urban–Rural Heterogeneity of Air Pollution in 35 Metropolitan Regions across China

Urbanization and air pollution are major anthropogenic impacts on Earth’s environment, weather, and climate. Each has been studied extensively, but their interactions have not. Urbanization leads to a dramatic variation in the spatial distribution of air pollution (fine particles) by altering surface properties and boundary-layer micrometeorology, but it remains unclear, especially between the centers and suburbs of metropolitan regions. Here, we investigated the spatial variation, or inhomogeneity, of air quality in urban and rural areas of 35 major metropolitan regions across China using four different long-term observational datasets from both ground-based and space-borne observations during the period 2001–2015. In general, air pollution in summer in urban areas is more serious than in rural areas. However, it is more homogeneously polluted, and also more severely polluted in winter than that in summer. Four factors are found to play roles in the spatial inhomogeneity of air pollution between urban and rural areas and their seasonal differences: (1) the urban–rural difference in emissions in summer is slightly larger than in winter; (2) urban structures have a more obvious association with the spatial distribution of aerosols in summer; (3) the wind speed, topography, and different reductions in the planetary boundary layer height from clean to polluted conditions have different effects on the density of pollutants in different seasons; and (4) relative humidity can play an important role in affecting the spatial inhomogeneity of air pollution despite the large uncertainties.


Introduction
The worldwide trend in urbanization has been ongoing for decades, altering landscapes and modifying the air quality. Urbanization is a key contributor to emissions of greenhouse gases and aerosol particles and changes in land cover, which has a significant effect on regional and temporal climate changes [1][2][3][4]. Although urbanization has brought about rapid economic growth, it has adverse impacts on air quality and human health [5][6][7][8][9]. Rapid urbanization can lead to the destruction of cultivated land, overcrowding, complex infrastructures, and extreme weather disasters [10].
Reduced vegetation in urban areas, urban construction materials and structures, and anthropogenic heat emissions are among the major factors leading to the urban heat island (UHI) [11][12][13]. Urban regions are generally warmer than their rural surroundings because of the surfaces changing from permeable and moist to impermeable and dry. The effect is proportional to the urban size and varies considerably due to the great differences in numerous influential factors [14,15]. The UHI can alter The spatial resolution of Landsat data is 30m, only summertime (June, July, and August) images before 2000 or in 2000, in 2010, and in 2015 are used to extract urban impervious surfaces and outline urban contours.
Only satellite-based AOD was used here. The MODIS Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm is an advanced algorithm which uses time series analyses and a combination of pixel-and image-based processing to improve the accuracies of cloud detection, aerosol retrievals, and atmospheric correction [35,36]. Its spatial and temporal resolutions are 1 km and 1 day, respectively. Linear interpolation is used to fill the missing values based on existing values when the data coverage is greater than 30%. For each city, five zones were selected based on the distance to the urban geometric center: Zone 1: 0-10 km, Zone 2: 11-20 km, Zone 3: 21-30 km, Zone 4: 31-40 km, and Zone 5: 41-50 km. The average AOD for each zone was then calculated.
The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is a two-wavelength polarization lidar that performs global profiling of aerosols and clouds in the troposphere and lower stratosphere. CALIOP is the primary instrument onboard the CALIPSO satellite, which has flown with the NASA A-train constellation of satellites since May 2006 [34].
The global emission inventory dataset called PKU-FUEL produced by Peking University was used to study emission differences between urban and rural areas. The inventory includes CO2, CO, PM2.5, PM10, BC, OC, SO2, NOx, NH3, total suspended particles, and polycyclic aromatic hydrocarbons. Satellite datasets used in this study include those from the Land Satellite Thematic Mapper/ Enhanced Thematic Mapper (Landsat), the Moderate Resolution Imaging Spectroradiometer [MODIS; aerosol optical depth (AOD)], and the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) satellite [34]. Ground-based observation data include meteorological stations, hourly fine-particulate matter (PM 2.5 ) concentration measurements, and daily sounding data.
The spatial resolution of Landsat data is 30m, only summertime (June, July, and August) images before 2000 or in 2000, in 2010, and in 2015 are used to extract urban impervious surfaces and outline urban contours.
Only satellite-based AOD was used here. The MODIS Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm is an advanced algorithm which uses time series analyses and a combination of pixel-and image-based processing to improve the accuracies of cloud detection, aerosol retrievals, and atmospheric correction [35,36]. Its spatial and temporal resolutions are 1 km and 1 day, respectively. Linear interpolation is used to fill the missing values based on existing values when the data coverage is greater than 30%. For each city, five zones were selected based on the distance to the urban geometric center: Zone 1: 0-10 km, Zone 2: 11-20 km, Zone 3: 21-30 km, Zone 4: 31-40 km, and Zone 5: 41-50 km. The average AOD for each zone was then calculated.
The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is a two-wavelength polarization lidar that performs global profiling of aerosols and clouds in the troposphere and lower stratosphere. CALIOP is the primary instrument onboard the CALIPSO satellite, which has flown with the NASA A-train constellation of satellites since May 2006 [34].
The global emission inventory dataset called PKU-FUEL produced by Peking University was used to study emission differences between urban and rural areas. The inventory includes CO 2 , CO, PM 2.5 , PM 10 , BC, OC, SO 2 , NO x , NH 3 , total suspended particles, and polycyclic aromatic hydrocarbons. This global emission inventory has been developed using a bottom-up approach with a 0.1 • × 0.1 • (about 10 km) spatial resolution and a monthly temporal resolution, covering the period from 1960 to Remote Sens. 2020, 12, 2320 4 of 22 2014 [37]. For each city, one urban window (0.1 • × 0.1 • ) inside urban contour in 2000 and four rural windows (0.1 • × 0.1 • ) that are not affected by urban expansion and outside urban contour in 2015.
Two sets of meteorological stations and PM 2.5 stations were selected in each city: one set located in the urban area of the city and the other set located in the rural area near the city. Meteorological parameters such as visibility, surface wind speed, temperature, and precipitation were observed at three-hourly intervals at each meteorological station. Hourly PM 2.5 concentrations were measured at each PM 2.5 measurement site (shown in Figure S1) from the Chinese national monitoring network. Hourly PM 2.5 data are published in real time by the Ministry of Ecology and Environment of China (http://106. 37.208.233:20035). The PM 2.5 concentration is measured by tapered element oscillating microbalance (TEOM) or beta attenuation monitor. High correlation between the satellite-based AOD and ground-based PM 2.5 observations were reported in previous studies [38][39][40]. One urban and one rural meteorological station was selected for each city. Meteorological (1100 and 1400 BJT) and PM 2.5 concentration (1300 and 1400 BJT) observed at near 1330 Beijing time (BJT) were linearly interpolated to the MODIS imaging time (1330 BJT) for matching purposes. According to the evaluation standard of the Ministry of Ecology and Environment of China, urban PM 2.5 concentration are classified into three levels of air quality: light pollution (PM 2.5 < 50 µg/m 3 ), moderate pollution (100 µg/m 3 < PM 2.5 < 150 µg/m 3 ), and heavy pollution (200 µg/m 3 < PM 2.5 ). The classification is based on 2-hour means of PM 2.5 concentration, different from other standards using annual mean or 24-hour mean (e.g., the standards of WHO, US, EU).
Four sounding stations, representative of their respective regions, were selected, i.e., Shenyang, Beijing, Xi'an, and Nanjing. The high-resolution radiosonde network of the L-band sounding system, developed by the China Meteorological Administration in 2011, provides fine-resolution profiles of temperature, pressure, RH, wind speed, and wind direction twice a day at 0800 BJT and 2000 BJT. The sounding quality was strictly controlled and adequate to characterize PBL features in China [41,42]. Unless noted otherwise, 2000 BJT soundings were used in this study, given our focus on the daytime effect.

Urban Impervious Surfaces and Urban Contours
Many previous studies have used nighttime stable-light data to extract urban areas, but their low spatial resolution limits the extraction accuracy. So here, the Landsat data were used to accurately identify urban impervious surfaces and outline urban contours. The difference between the normalized difference build-up index (NDBI) and the soil-adjusted vegetation index (SAVI), i.e., NDBI − SAVI, is used to extract urban impervious surfaces because this difference can effectively differentiate urban impervious surfaces from other land-use types [43]: where L is the soil adjustment factor whose value is 0.5, and ρ n is the nth Landsat reflectance band. Urban impervious surfaces were extracted using different thresholds ranging from 0.1 to 0.3, then urban boundaries were determined based on the urban impervious surfaces. Google Earth and a land-use map with a 1:100,000 scale from the Resource and Environment Science Data Center of the Chinese Academy of Sciences verified the results. The physical contours of urban areas were extracted based on the difference in the underlying surfaces of urban and rural areas. Urban contours in three periods were outlined at each city, which were used to distinguish urban and rural areas. Figure 2 shows an example diagram of Nanjing.

Urban Shape
The Boyce-Clark shape index (SBC) has been used to reveal the urban morphology through comparisons with a standard shape [44]. The principle behind this index is to compare the urban shape with a standard circle, then calculate a relative shape index. This is also called the radius shape index on account of radii being involved in the algorithm, expressed as follows: where is the length from the vantage point to the boundary, and n is the number of radii with an equal-angle difference. The vantage point of urban areas can be the central business district or the centroid of the urban contour. For example, when n is equal to 16, the angle between adjacent radii is 22.5°, and when n is equal to 32, the angle between adjacent radii is 11.25°. The minimum SBC is 0, which is the SBC of a standard circle. It represents the highest land-use efficiency and the most compact urban area. Larger SBC values generally mean lower land-use efficiencies and less compact urban areas. Urban contours of 35 cities extracted from Landsat in 2015 were used to calculate SBC here.

Research Windows for Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)
CALIPSO data can provide vertical variation information, unlike MODIS data. Two metropolitan clusters were selected as research areas for CALIPSO: the Greater Beijing Metropolitan Area (GBMA) and the Yangtze River Delta (YRD). Both GBMA and YRD have experienced dramatic economic development and rapid urbanization over the past three decades. The GBMA includes Beijing, Tianjin, and some industrial cities in Hebei province, and the YRD region includes Shanghai, Suzhou, Wuxi, Changzhou, Zhenjiang, Yangzhou, Nanjing, and five other prosperous cities [45,46]. Each research area has one urban research window and one rural research window.

Calculation of the Planetary Boundary Layer Height
Five cities with sounding stations that are part of the China Meteorological Administration's radiosonde network of L-band sounding systems were selected: Beijing, Shenyang, Chengdu, Xi'an, and Nanjing. Radiosonde profiles from 1 January 2013 to 31 December 2015 were analyzed. The bulk Richardson number ( ) has been used to estimate the PBLH [41,47,48], which is also done here. The is defined as the ratio of turbulence associated with buoyancy to that induced by mechanical shear and is expressed as where denotes the height above the ground, denotes the surface, is the acceleration due to gravity, is the virtual potential temperature, and are the components of wind speed, and *

Urban Shape
The Boyce-Clark shape index (SBC) has been used to reveal the urban morphology through comparisons with a standard shape [44]. The principle behind this index is to compare the urban shape with a standard circle, then calculate a relative shape index. This is also called the radius shape index on account of radii being involved in the algorithm, expressed as follows: where r i is the length from the vantage point to the boundary, and n is the number of radii with an equal-angle difference. The vantage point of urban areas can be the central business district or the centroid of the urban contour. For example, when n is equal to 16, the angle between adjacent radii is 22.5 • , and when n is equal to 32, the angle between adjacent radii is 11.25 • . The minimum SBC is 0, which is the SBC of a standard circle. It represents the highest land-use efficiency and the most compact urban area. Larger SBC values generally mean lower land-use efficiencies and less compact urban areas. Urban contours of 35 cities extracted from Landsat in 2015 were used to calculate SBC here.

Research Windows for Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)
CALIPSO data can provide vertical variation information, unlike MODIS data. Two metropolitan clusters were selected as research areas for CALIPSO: the Greater Beijing Metropolitan Area (GBMA) and the Yangtze River Delta (YRD). Both GBMA and YRD have experienced dramatic economic development and rapid urbanization over the past three decades. The GBMA includes Beijing, Tianjin, and some industrial cities in Hebei province, and the YRD region includes Shanghai, Suzhou, Wuxi, Changzhou, Zhenjiang, Yangzhou, Nanjing, and five other prosperous cities [45,46]. Each research area has one urban research window and one rural research window.

Calculation of the Planetary Boundary Layer Height
Five cities with sounding stations that are part of the China Meteorological Administration's radiosonde network of L-band sounding systems were selected: Beijing, Shenyang, Chengdu, Xi'an, and Nanjing. Radiosonde profiles from 1 January 2013 to 31 December 2015 were analyzed. The bulk Richardson number (R i ) has been used to estimate the PBLH [41,47,48], which is also done here. The R i is defined as the ratio of turbulence associated with buoyancy to that induced by mechanical shear and is expressed as where z denotes the height above the ground, s denotes the surface, g is the acceleration due to gravity, θ v is the virtual potential temperature, u and v are the components of wind speed, and u * is the surface friction velocity. Here, u * can be ignored because its magnitude is small [48]. Previous theoretical and laboratory studies [49] have suggested that when R i is smaller than the critical value (0.25), the laminar flow becomes unstable. The lowest level z at which interpolated R i crosses the critical value of 0.25 is thus referred to as the PBLH in this study, similar to the criterion used by others [41]. Figure 3a shows detailed mean visibilities of each city in urban and rural areas in summer and winter under polluted conditions. Figure 3b shows the overall mean visibilities of all cities. Only visibilities associated with RH levels less than 85% were analyzed to eliminate the influence of water vapor. Urban visibilities of less than 10 km are considered as grossly polluted in summer (June, July, and August) and winter (December, January, and February).

Visibility Differences between Urban and Rural Areas
is the surface friction velocity. Here, * can be ignored because its magnitude is small [48]. Previous theoretical and laboratory studies [49] have suggested that when is smaller than the critical value (0.25), the laminar flow becomes unstable. The lowest level at which interpolated crosses the critical value of 0.25 is thus referred to as the PBLH in this study, similar to the criterion used by others [41]. Figure 3a shows detailed mean visibilities of each city in urban and rural areas in summer and winter under polluted conditions. Figure 3b shows the overall mean visibilities of all cities. Only visibilities associated with RH levels less than 85% were analyzed to eliminate the influence of water vapor. Urban visibilities of less than 10 km are considered as grossly polluted in summer (June, July, and August) and winter (December, January, and February).

Visibility Differences between Urban and Rural Areas
Averaging over all cities, the wintertime urban mean visibility under polluted conditions is ~5.8 km (standard deviation, or std, = 1.1 km), and the wintertime rural mean visibility is ~10.7 km (std = 5.3 km). In the summertime, the urban mean visibility is ~7.0 km (std = 1.9 km), and the rural mean visibility is ~16.4 km (std = 6.8 km). The overall differences in visibility between urban and rural areas are ~9.4 km and ~4.9 km for summertime and wintertime, respectively.  Figure 4 shows the mean difference in AOD (∆ = − ) between the urban and rural areas of each city and all cities. Summertime ∆ s are larger than those in winter in all cities except for Chongqing. Figure 5 shows the variation trends of mean AOD as a function of distance from the urban geometrical center of each city in winter and summer. As the distance from the urban geometrical center increases, summertime AODs decrease more rapidly than wintertime AODs. The range of the decrease in summer is greater than that in winter.  Averaging over all cities, the wintertime urban mean visibility under polluted conditions is 5.8 km (standard deviation, or std, = 1.1 km), and the wintertime rural mean visibility is~10.7 km (std = 5.3 km). In the summertime, the urban mean visibility is~7.0 km (std = 1.9 km), and the rural mean visibility is~16.4 km (std = 6.8 km). The overall differences in visibility between urban and rural areas are~9.4 km and~4.9 km for summertime and wintertime, respectively. Figure 4 shows the mean difference in AOD (∆AOD = AOD urban − AOD rural ) between the urban and rural areas of each city and all cities. Summertime ∆AODs are larger than those in winter in all cities except for Chongqing. Figure 5 shows the variation trends of mean AOD as a function of distance from the urban geometrical center of each city in winter and summer. As the distance from the urban geometrical center increases, summertime AODs decrease more rapidly than wintertime AODs. The range of the decrease in summer is greater than that in winter. Figures 4 and 5 indicate that the ∆AOD in summer is larger than that in winter. The overall mean ∆AODs are 0.175 in summer and 0.07 in winter, and the relative difference between summer and winter is 0.105. Moreover, the AOD variations of different cities in Figure 5 are inconsistent, because of the joint complex effects of meteorological conditions, topography, and aerosol types.

Aerosol Optical Depth (AOD) Differences between Urban and Rural Areas
The spatial difference in CALIPSO aerosol distributions from urban to rural areas was also analyzed. Figure 6 shows the frequency of occurrence of aerosols at different altitudes over the GBMA and YRD regions from 2006 to 2015. In both regions, the occurrence frequency differs more in summer than in winter, consistent with the column total AOD from MAIAC. Aerosols reach higher into the atmosphere in summer than in winter, i.e., higher than 3 km (GBMA) and 2.5 km (YRD) in summer, and lower than 2 km over both regions in winter. The occurrence frequency as determined by CALIPSO at lower altitudes (<1.5 km) is larger in urban areas than in rural areas in summer (Figure 6a,c). This is not as obvious in winter (Figure 6b,d). In rural areas, larger occurrence frequencies in summer are located at higher altitudes (0.5-2.5 km), whereas they are at lower altitudes (<1 km) in winter. The rural frequency of occurrence is close to, or even larger than, the urban frequency of occurrence at higher altitudes in summer. Urban-rural circulation may explain this. Because of higher temperatures, urban areas act as convergence zones, bringing in near-surface pollutants from rural areas that are then lifted by updrafts to higher altitudes and dispersed to rural surrounding areas, especially in summer. This is more clearly seen at the border between urban and rural areas (Figure 6a,c) This process diminishes the vertical gradient of pollution in the boundary areas of a city. Since updrafts are stronger in summer than in winter, pollutants can be carried to higher altitudes, leading to more aerosols at higher altitudes in summer than in winter.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 22 that the ∆ in summer is larger than that in winter. The overall mean ∆ s are 0.175 in summer and 0.07 in winter, and the relative difference between summer and winter is 0.105. Moreover, the AOD variations of different cities in Figure 5 are inconsistent, because of the joint complex effects of meteorological conditions, topography, and aerosol types.   in winter and summer, respectively. The overall mean ∆ calculated using data in winter and summer are shown as black and green lines, respectively. The city marked with an asterisk (Chongqing) shows results opposite to those of the other cities.

PM 2.5 Differences between Urban and Rural Areas
Figure 7a-c shows detailed PM 2.5 differences for each city. Figure 7d shows summertime and wintertime overall PM 2.5 differences (∆PM 2.5 = PM 2.5−urban − PM 2.5−rural ) at each pollution level of all cities using data from 2013 to 2015. Summertime PM 2.5 differences are generally larger than wintertime PM 2.5 differences, with greater increases in magnitude than wintertime PM 2.5 differences as PM 2.5 increases.
The results of Section 3 indicate that the SIAP is found to exhibit a pronounced seasonality that is not recognized before, namely, that summertime SIAP > wintertime SIAP, it can reveal the linkage of air pollution to the urban heat island, clouds, and precipitation, among others. Combining multi-source datasets shows the 3-D distribution of air pollution between urban and rural areas that have rarely been investigated before (mostly 2-D), and this is important for validating future model simulations.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 22   ) at each pollution level of all cities using data from 2013 to 2015. Summertime . differences are generally larger than wintertime . differences, with greater increases in magnitude than wintertime . differences as . increases. The results of Section 3 indicate that the SIAP is found to exhibit a pronounced seasonality that is not recognized before, namely, that summertime SIAP > wintertime SIAP, it can reveal the linkage of air pollution to the urban heat island, clouds, and precipitation, among others. Combining multisource datasets shows the 3-D distribution of air pollution between urban and rural areas that have rarely been investigated before (mostly 2-D), and this is important for validating future model simulations. Figure 8 presents the spatial distribution of total PM2.5 emissions in the urban and rural areas of each city in summer and winter. For all cities, urban PM2.5 emissions are larger than rural emissions in both winter and summer. Summertime PM2.5 emissions from most cities located south of 32° are slightly larger than in winter. The opposite is the case for cities located north of 32° where emissions from heating in winter may make a difference. For most cities, urban-rural differences in PM2.5 emissions in summer are slightly larger than in winter. Figure 9 shows the urban-rural emission differences of precursors (including NOx, SO2, and NH3), which promote the formation of particulate matter. Results also show that the urban-rural differences in precursors in summer are larger than in winter in most cities. Overall, the seasonal differences are not obvious, suggesting that the emission difference is a factor but not the only one causing the different spatial inhomogeneities of air pollution in different seasons. It should be noted that the AOD differences in a few cities are inconsistent with emission differences, because AOD is not only affected by PM2.5, but also by PBL height, RH, and larger particles (e.g., PM10, dust).  emissions from heating in winter may make a difference. For most cities, urban-rural differences in PM 2.5 emissions in summer are slightly larger than in winter. Figure 9 shows the urban-rural emission differences of precursors (including NO x , SO 2 , and NH 3 ), which promote the formation of particulate matter. Results also show that the urban-rural differences in precursors in summer are larger than in winter in most cities. Overall, the seasonal differences are not obvious, suggesting that the emission difference is a factor but not the only one causing the different spatial inhomogeneities of air pollution in different seasons. It should be noted that the AOD differences in a few cities are inconsistent with emission differences, because AOD is not only affected by PM 2.5 , but also by PBL height, RH, and larger particles (e.g., PM 10 , dust).    Figures 10a,b show the relationship between SBC and ∆ in summer and winter, respectively. In summer, ∆ decreases as SBC increases. In winter, there is no clear relationship between them. Figure 10c shows the extent of changes in AOD every 10 km as a function of SBC. Results show that the extent of decreases in AOD reduces with increasing SBC in summer. The extent of decreases in AOD with increasing SBC barely changes in winter. Figure 10d shows the relationship between PM2.5 concentration and building density. PM2.5 concentrations first increase then decrease as the building density increases in both summer and winter, but the correlation in summer is more significant than in winter. The above results indicate that urban structure has a more obvious relationship with the spatial inhomogeneity of air pollution in summer than in winter. Low SBC values indicate simpler shapes and more compact urban areas, which serve as a barrier against pollutants being transported elsewhere, causing greater differences between the pollution levels of  Figure 10a,b show the relationship between SBC and ∆AOD in summer and winter, respectively. In summer, ∆AOD decreases as SBC increases. In winter, there is no clear relationship between them. Figure 10c shows the extent of changes in AOD every 10 km as a function of SBC. Results show that the extent of decreases in AOD reduces with increasing SBC in summer. The extent of decreases in AOD with increasing SBC barely changes in winter. Figure 10d shows the relationship between PM 2.5 concentration and building density. PM 2.5 concentrations first increase then decrease as the building density increases in both summer and winter, but the correlation in summer is more significant than in winter. The above results indicate that urban structure has a more obvious relationship with the spatial inhomogeneity of air pollution in summer than in winter. Low SBC values indicate simpler shapes and more compact urban areas, which serve as a barrier against pollutants being transported elsewhere, causing greater differences between the pollution levels of urban and rural areas. High SBC values indicate more complex shapes and less compact urban areas, which favors the dispersion of pollutants, leading to lower differences in the pollution levels between urban and rural areas [50,51]. The urban structure affects the spatial distribution of air pollution in the summer because the compactness of urban areas plays an important role in the dispersion of pollutants in the more active PBL. In winter, the PBL is very stable, so the influence of the compactness of urban areas likely weakens, lessening the effect.

Urban Structure Effect
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 urban and rural areas. High SBC values indicate more complex shapes and less compact urban areas, which favors the dispersion of pollutants, leading to lower differences in the pollution levels between urban and rural areas [50,51]. The urban structure affects the spatial distribution of air pollution in the summer because the compactness of urban areas plays an important role in the dispersion of pollutants in the more active PBL. In winter, the PBL is very stable, so the influence of the compactness of urban areas likely weakens, lessening the effect.

Meteorological and Topographical Effect
To gain further insight into the impact of the PBL on air quality, their relationship is investigated more rigorously. Note that pollutants mainly exist within the PBL. A lower (higher) PBLH means less (more) space available for aerosols, so changes in the PBLH will affect the spatial inhomogeneity of air pollution to varying degrees [21,52]. Figure 11a-d show the mean PBLHs in summer and winter at four cities at three pollution levels. Figure 11e shows PBLH differences between clean and heavy polluted conditions in the four cities in summer and winter. The PBLH decreases persistently as the air quality changes from clean, moderately polluted to heavily polluted, and the extent of changes in PBLH is more in summer (~0.2-0.9 km) than in winter (<0.4 km). In absolute terms, the PBLH decreases in summer (~0.54-0.99 km) are more than twice that in winter (<0.35 km) at these four cities when air quality changes from clean to polluted conditions. In relative terms, the PBLH decreases 51%, 50%, 52%, and 35% in summer at the four cities, respectively, but 21%, 30%, 38%, and 26% in winter. A greater PBLH reduction leads to greater PM 2.5 concentrations within the PBL. The greater reduction in the PBLH in summer enhances the spatial inhomogeneity of air pollution more significantly than that in winter.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 22 To gain further insight into the impact of the PBL on air quality, their relationship is investigated more rigorously. Note that pollutants mainly exist within the PBL. A lower (higher) PBLH means less (more) space available for aerosols, so changes in the PBLH will affect the spatial inhomogeneity of air pollution to varying degrees [21,52]. Figure 11a-d show the mean PBLHs in summer and winter at four cities at three pollution levels. Figure 11e shows PBLH differences between clean and heavy polluted conditions in the four cities in summer and winter. The PBLH decreases persistently as the air quality changes from clean, moderately polluted to heavily polluted, and the extent of changes in PBLH is more in summer (~0.2-0.9 km) than in winter (<0.4 km). In absolute terms, the PBLH decreases in summer (~0.54-0.99 km) are more than twice that in winter (<0.35 km) at these four cities when air quality changes from clean to polluted conditions. In relative terms, the PBLH decreases 51%, 50%, 52%, and 35% in summer at the four cities, respectively, but 21%, 30%, 38%, and 26% in winter. A greater PBLH reduction leads to greater PM2.5 concentrations within the PBL. The greater reduction in the PBLH in summer enhances the spatial inhomogeneity of air pollution more significantly than that in winter.  Figure 12 shows the variation of ΔPM2.5 with wind speed. In Beijing, ΔPM2.5 increases with increasing wind speed in summer, but decreases in winter. The prevailing southeast wind in summer carries pollutants from surrounding area to urban area, and the northwest mountain areas prevent diffusion of urban pollutants. This process facilitates the accumulation of pollutants in urban areas and increase urban-rural pollution difference. However, in winter, the prevailing northwest wind from mountain area is cleaner and pushes urban pollutants to surrounding areas, reducing urbanrural pollution differences. The prevailing wind directions in Xi'an are same as that in Beijing in  Figure 12 shows the variation of ∆PM 2.5 with wind speed. In Beijing, ∆PM 2.5 increases with increasing wind speed in summer, but decreases in winter. The prevailing southeast wind in summer carries pollutants from surrounding area to urban area, and the northwest mountain areas prevent diffusion of urban pollutants. This process facilitates the accumulation of pollutants in urban areas and increase urban-rural pollution difference. However, in winter, the prevailing northwest wind from mountain area is cleaner and pushes urban pollutants to surrounding areas, reducing urban-rural pollution differences. The prevailing wind directions in Xi'an are same as that in Beijing in summer and winter. There are also mountains near Xi'an, but they are situated south of the urban area. So the combined effects of wind and topography in Xi'an is opposite to that in Beijing. There is no continuous mountain area around Nanjing. Winds from any directions help disperse urban pollutants and decrease urban-rural pollution differences. In Shenyang, the prevailing winds in summer and winter are southeast wind and northeast wind, respectively. The cleaner winds from eastern mountain area help remove urban pollutants. The above results indicate that the combination of wind and topography can affect the spatial inhomogeneity of air pollution.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 22 summer and winter. There are also mountains near Xi'an, but they are situated south of the urban area. So the combined effects of wind and topography in Xi'an is opposite to that in Beijing. There is no continuous mountain area around Nanjing. Winds from any directions help disperse urban pollutants and decrease urban-rural pollution differences. In Shenyang, the prevailing winds in summer and winter are southeast wind and northeast wind, respectively. The cleaner winds from eastern mountain area help remove urban pollutants. The above results indicate that the combination of wind and topography can affect the spatial inhomogeneity of air pollution.  Figure 13 shows the ΔAODs of different terrain cities. The results show that the ΔAODs in summer are larger than in winter in all such cities. Comparing different categories of cities in Figure  13, ΔAOD has smaller changes in summer than in winter. This indicates that the effect of topography in winter is more obvious than that in summer.  Figure 13 shows the ∆AODs of different terrain cities. The results show that the ∆AODs in summer are larger than in winter in all such cities. Comparing different categories of cities in Figure 13, ∆AOD has smaller changes in summer than in winter. This indicates that the effect of topography in winter is more obvious than that in summer.

Humidity Effects
Humidity facilitates both particle growth and new particle formation, having a stronger effect when RH levels are high [53][54][55]. New particle formation is highly localized and extremely tiny (from nano-to tens of nano-meter). This process tends to have little effect on AOD and PM 2.5 that are most strongly influenced by particles >0.1 m. However, particle growth can largely enhance the scattering efficiency and increase particle size and AOD, also affecting nucleation and coagulation, and thus dry and wet depositions, which would influence the mass concentration. Humidity exerts a considerable influence on aerosols and generally increases both AOD and PM 2.5 [54,56,57]. Due to the spatial variation in RH, humidity can influence the spatial inhomogeneity of air pollution in terms of mass concentration, particle diameter, and optical depth [55,[58][59][60].
Nanjing, Shenyang, and Xi'an in summer (green curves) and winter (black curves). The ΔPM2.5 concentration is the difference in PM2.5 concentration between urban and rural areas. Figure 13 shows the ΔAODs of different terrain cities. The results show that the ΔAODs in summer are larger than in winter in all such cities. Comparing different categories of cities in Figure  13, ΔAOD has smaller changes in summer than in winter. This indicates that the effect of topography in winter is more obvious than that in summer.  Table 1 summarizes the correlation coefficients between ∆AOD and ∆RH for all cities. Twenty-eight cities have positive relationships between ∆AOD and ∆RH in summer, while only eleven cities have positive relationships in winter. Figure 14 shows the variations in ∆PM 2.5 concentration as a function of RH in Beijing, Nanjing, Shenyang, and Xi'an. Results indicate that ∆PM 2.5 continuously increases as RH increases in summer in all cities, while the variations in ∆PM 2.5 differ with different values of RH in winter. In winter, the ∆PM 2.5 variations in Beijing and Xi'an are different from other cities. ∆PM 2.5 changes at high RH in Xi'an are opposite to those in Beijing, which may be caused by different aerosol compositions and properties between urban and rural areas. Urban particles grow more than rural particles under high RH in Beijing, which is contrary to Xi'an. This may lead to different relationship between ∆AOD and ∆RH in different seasons in Table 1 in a few cities. Figure 15 shows the mean RH in each city in summer and winter. Twenty-seven cities have higher RHs in summer than in winter. The other eight cities, mainly located in the northernmost and southernmost parts of China, have lower RHs in summer than in winter. Previous studies have shown that urban aerosols are more hygroscopic than rural aerosols, especially in summer, because of the higher proportion of inorganic aerosols and higher RH in summer (shown in Figure 15). The spatial inhomogeneity of AOD or pollutants in summer is thus more sensitive to the effects of RH [61][62][63][64][65]. Furthermore, the differences in aerosol hygroscopicity can affect the spatial inhomogeneity of AOD between urban and rural areas. Differences in RH affect the process of particle growth, serving as a potential factor explaining the spatial differences in aerosol loading. Note that the effects of humidity on aerosols are still highly uncertain. In summary, the above analyses suggest that the effects on seasonal SIAP are complex, likely driven by multi-factors that vary with season, including emission, urban structure, PBL, and relative humidity.

Discussion for Potential Factors
The interaction mechanism between urbanization and air pollution is very complicated. Here, the four potential factors we proposed provide new insights into this subject. Meanwhile, they raise more questions than what can be addressed in this study, due to the limitations of currently available observational data. For the sake of a future study, the limitations of the approach and uncertainties or ambiguities in the interpretation of the results are stated here. (1) For the emission effect, the spatial resolution of the emission data is about 10 km. While it is higher than some related products, it is still not fine enough to resolve any variations on urban-rural scales for which higher resolution data would be needed.
(2) For the urban structure effect, the relationship between SBC and ∆AOD is faint. It is thus impossible to solely extract the urban effects, as other factors disturb the signal of urban structures. Resolving this issue entails model simulations (e.g., LES-Large eddy simulation or WRF-Chem simulation) [66]. (3) While the influences of PBLH, wind speed, and topography are analyzed here, quantification of their respective importance is still wanting. (4) For the humidity effect, a handful of studies have been conducted, there still exit large uncertainties, especially concerning urban-rural difference. Understanding the effects of humidity merit further comprehensive observations of aerosol chemical composition between urban and rural areas.   In summary, the above analyses suggest that the effects on seasonal SIAP are complex, likely driven by multi-factors that vary with season, including emission, urban structure, PBL, and relative humidity.

Discussion for Potential Factors
The interaction mechanism between urbanization and air pollution is very complicated. Here, the four potential factors we proposed provide new insights into this subject. Meanwhile, they raise more questions than what can be addressed in this study, due to the limitations of currently available observational data. For the sake of a future study, the limitations of the approach and uncertainties or ambiguities in the interpretation of the results are stated here. (1) For the emission effect, the spatial resolution of the emission data is about 10 km. While it is higher than some related products, it is still not fine enough to resolve any variations on urban-rural scales for which higher resolution data would be needed. (2) For the urban structure effect, the relationship between SBC and ΔAOD is faint. It is thus impossible to solely extract the urban effects, as other factors disturb the signal of urban structures. Resolving this issue entails model simulations (e.g., LES-Large eddy simulation or WRF-Chem simulation) [66]. (3) While the influences of PBLH, wind speed, and topography are analyzed here, quantification of their respective importance is still wanting. (4) For the humidity effect, a handful of studies have been conducted, there still exit large uncertainties, especially concerning urban-rural difference. Understanding the effects of humidity merit further comprehensive observations of aerosol chemical composition between urban and rural areas.

Conclusions
Satellite data and ground-based observation data from multiple sources were used to analyze the spatial inhomogeneity of air pollution between the urban and rural areas of 35 cities in China in different seasons. The results obtained from these multi-source datasets are consistent and provide observational analysis support for previous model simulation work [30]. Results show that the spatial inhomogeneity of air pollution between urban and rural areas is larger in summer than in winter. Urban pollution is often more severe than rural pollution in summer. However, in winter, both urban and rural pollution are severe, and rural pollution is even more severe than urban pollution in few cities.
Emissions, urban structures, meteorology and topography, and humidity are the potential reasons for why the spatial inhomogeneity of air pollution is larger in summer: (1) For most cities, the urban-rural differences in PM 2.5 and precursor-gas emissions in summer are larger than in winter but not obviously. This indicates that besides the emission effect, there are also other factors causing different spatial inhomogeneities of air pollution in different seasons.
(2) The effect of the urban structure is more significant in summer than in winter. The compactness of urban areas has a more obvious effect on the dispersion of pollutants in summer because of an active planetary boundary layer, but not in winter.
(3) The reduction in planetary boundary layer height over urban areas is larger in summer than in winter from clean to heavy polluted conditions. This causes more significant influence on urban pollution concentrations in summer than in winter, as is its impact on the spatial inhomogeneity of air pollution. Moreover, both wind speed and topography affect the spatial inhomogeneity of air pollution.
(4) For most cities, higher relative humidity can cause the larger spatial inhomogeneities of air pollution in summer and winter due to particle growth.
Note that these explanations are not independent, i.e., there are interactions between them, influencing each other. However, their respective importance cannot be quantitatively assessed by analyzing observation data alone. Numerical simulations will be conducted to gain further insights. This study mainly focused on physical factors. Synchronous experiments of urban and rural areas are needed in the future to gain a deep insight into related problems. The chemical transformation of pollutants and socioeconomic factors also need to be investigated to understand their impacts on the spatial inhomogeneity of air pollution.