Aerosol optical depth regime over megacities of the world

Currently, 55 % of the world's population resides in urban areas and this number is projected to increase to 70 % by 2050. Urban agglomerations with a population over 10 million, characterized as megacities, are expected to be more than 100 by 2100. Such large concentrations of population could boost creativity and economic progress, but also raises several environmental challenges such as air quality degradation. In this study, we investigate the spatial and temporal variability of urban aerosol state of 81 cities with a population over 5 million, relying on daily satellite-based aerosol optical depth (AOD) retrievals, derived at fine spatial resolution ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$), over an 18-year period spanning from 2003 to 2020.

According to our results, the lowest long-term mean AOD values worldwide were found in European and American cities (from 0.08 to 0.20). For almost all African and Asian cities, mean AOD ranged from 0.25 up to 0.90, but a considerable dust aerosol contribution (up to 70 %) was found for some of them with associated mean dust optical depth (DOD) values reaching up to 0.4. Mostly Chinese and Indian cities tend to have higher mean AOD values in the areas surrounding their center, while the opposite was found for most of the cities in the rest of the world. High intraannual AOD variability was revealed for the eastern American cities, while lower values were found in Chinese, eastern Indian and the eastern Mediterranean cities. During the study period, statistically significant negative AOD decadal trends were found for East Asian, European and North American cities, with the greatest decrease of −0.1 to −0.3 per decade recorded for the Chinese cities, in which the maximum mean AODs (0.45–0.91) are observed. In most of the US cities, where low mean AOD <0.17 was recorded, considerable declining AOD trends were found (−30 % to −50 % per decade). For the rest of Asian, African and South American cities, statistically significant AOD increase was found, with the greatest values of +0.07 to +0.16 per decade recorded for Indian cities. In Bengaluru (India), it is reported the lowest mean AOD value (0.2) and the maximum AOD increase (+69 %), which may be partially attributed to the population growth over the study period. The agreement of the satellite-derived AOD trends against those obtained from ground-based AERONET measurements was examined. For ground-based stations within the geographical limits of the contiguous urban area of the examined cities, a 0.93 correlation for the long-term means of AOD was found and ∼75 % of the derived trends agreed in sign. It was found that the spatial homogeneity within the examined satellite domain and the location of the surface station were key factors that determined their agreement.

The present study highlights the vital and essential contribution of spaceborne products to monitor aerosol burden over megacities of the planet towards fulfilling the United Nations Sustainable Development Goal of “sustainable cities and communities”, dealing with urban air quality.

1 Introduction

Since 2018, about 55 % of the world's population has been living in urban areas and the current projections suggest that this percentage will rise to 60 % by 2030 and to 70 % by 2050 (UN, 2019b). As a consequence, the population growth will impose difficulties in implementing the United Nations (UN) Agenda for Sustainable Development Goals (SDGs), including SDG 11 related to “making cities and human settlements inclusive, safe, resilient and sustainable”. Megacities are defined as the cities with more than 10 million inhabitants. According to UN (2018a), there are 33 megacities worldwide and this number is expected to increase to 43 in 2030. Despite the fact that ∼7 % of the global population in 2018 resides in those megacities (UN, 2018a), their continuous increase both in size and number bring them in the spotlight, as they will accommodate an increased share of the world's population (projected to ∼10 % by 2030). This population growth of cities, and especially those with more than 5 million inhabitants, raises urgent and critical environmental issues, like air quality (WHO, 2021). The worst pollutant affecting the megacities is the suspended particulate matter or aerosols. This is of particular concern, as high levels of aerosols are known to be related to increased morbidity and mortality rates, and in many of the megacities in developing countries healthcare for acute cases is less proficient than in developed countries. Particularly, in 2019, cities' air pollution constituted the 4th leading risk factor for early death at a global scale (HEI, 2020).

Many countries worldwide have enforced policies in recent decades to reduce anthropogenic aerosol emissions in urban areas towards mitigating the aerosol adverse health effects. These policies include the transition of thermal engines technologies, pollution fees on big industries and schemes for air quality control. The exact policies and the level of implementation largely differ between countries, but it is the main cause of decrease in aerosol concentration in some megacities. WHA69 (WHO, 2016) set a global roadmap and targets for reduction of air pollution related deaths, urban air quality and clean energy. In the United States (US), the first regulations were legislated in 1970, and Environmental Protection Agency (EPA) reported a drop of ∼40 % in aerosol related air pollution in major US cities in the last half century (DeMocker, 2003). In the European Union (EU) countries, a decrease of ∼29 % in aerosols has been reported since 1970 as a result of policy measures and technological improvements (Turnock et al., 2016). In South and Central American countries, a poor emission control has been reported, resulting in 150 million people living in urban areas with poor air quality (UNEP, 2016). Jin et al. (2016) summarized the policies and their implementation over the last three decades for China, and as a consequence China's anthropogenic emissions markedly declined in the last decade (Zheng et al., 2018). In 2019, India launched the National Clean Air Programme, with the objective to reduce in 2024 the particulate air pollution by 20 %–30 % with respect to 2017 levels, in 122 cities (Ganguly et al., 2020). The regulations for air pollution for many African cities are weak or non-existent (Abera et al., 2020), which will be challenging for their future air quality, since many African cities have an unprecedented population growth. So, it is of vital importance to understand the trends of aerosol loads and their spatial variability over those great population agglomerations, this way assessing the effectiveness of air pollution emission regulations.

Aerosol-related air pollution can be quantified, in optical terms, through aerosol optical depth (AOD), which is the most comprehensive variable for assessing the aerosol load of the atmospheric column. Ground-based sun photometers have been deployed during the last 20–25 years, providing long-term AOD measurements at established global/regional networks such as AERONET (Holben et al., 1998), GAW-PFR (Kazadzis et al., 2018) and SKYNET (Nakajima et al., 2020). Although the most precise method for monitoring AOD is provided by surface sun photometers, these measurements are scarce, lacking full spatial and temporal coverage. On the contrary, satellite remote sensing is a powerful tool for monitoring AOD around the globe (Kaufman et al., 2002) at considerable accuracy, almost on a daily basis and at relatively fine spatial resolution. Despite the fact that the quality of spaceborne AOD over urban surfaces depends strongly on the limitations of the retrieval algorithms (Gupta et al., 2016), it is an aerosol parameter available worldwide at high spatial and temporal resolution.

There is a large number of studies dealing with AOD trends at local (e.g., Raptis et al., 2020; Vohra et al., 2021), regional (e.g., Che et al., 2019; Cherian and Quaas, 2020; Zhao et al., 2017) and global scale (e.g., Buchholz et al., 2021; Gupta et al., 2022; Logothetis et al., 2021; Wei et al., 2019a). AOD trends (either from spaceborne or from ground-based observations) in the last two decades indicate robust regional patterns, with decreasing aerosol loads over Europe and US since 2000 and reversed to decreasing since 2010 for East Asia and continuously increasing over southern Asia, which is also supported by in situ particulate concentration measurements (Gulev et al., 2021 and references therein). Nevertheless, limited studies have focused on aerosol regime over megacities, which are the major anthropogenic aerosol sources worldwide, and hence studies focused on them provide the direct link between aerosol load and emissions variability. In order to conduct studies like this, it is necessary to have available a dense surface-based network or to take advantage of the satellite remote sensing coverage capabilities. Alpert et al. (2012), relying on spaceborne aerosol measurements (MODIS and MISR), investigated the AOD tendencies over the 189 largest world cities for the 8-year period between 2002–2010. According to their findings, increasing AOD trends were found for the largest cities in the Indian subcontinent, the Middle East and North China, whereas opposite trends were evident in European, northeast US and southeast Asian megacities. The aim of this study is to use up-to-date and state-of-the art spaceborne aerosol retrievals (AOD and dust optical depth – DOD) to study local and transported aerosols that can affect air quality in large cities around the globe. Based on the idea of Alpert et al. (2012), we used a finer spatial resolution (0.1^∘) aerosol product that is able to reveal aerosol air pollution disparities within megacities and will allow for a more detailed AOD spatial analysis at urban scale and for a period from 2003 to 2020, in order to have more robust trends. Analyzing the aerosol spatiotemporal variability over megacities is of great importance for air quality and human health.

The present study aims to investigate the aerosol regime of megacities by addressing the following scientific questions:

• What is the spatial variability of AOD in megacities at different geographical locations?

• What is the temporal variability?

• Are ground-based AOD measuring stations adequate for reflecting the spatiotemporal changes of AOD in megacity scales?

2 Data and methods

2.1 Megacities information

Megacities are defined as cities with a population of more than 10 million, but in our analysis, we include cities with populations between 5 and 10 million as well, due to their potential to become megacities in the coming decades, as estimated by UN projections (UN, 2018a). Table 1 summarizes the number of cities according to their population by 2018 and future projections. According to UN projections, 10 cities with populations between 5 and 10 million are expected to become megacities in the near future. Here, we focus on 81 cities, with more than 5 million inhabitants as reported in 2018, which are listed in Table A1 (UN, 2018a, 2019a) and their geographical location is depicted in Fig. 1.

It should be mentioned that the UN data correspond mostly (65 %) to urban agglomerations, for which the city's boundaries are defined as the extent of the contiguous urban area (or built-up area). A small part refers to the metropolitan areas (25 %) and even smaller (10 %) to the city proper. For a major metropolitan area (MMA), more than one city is close together, and in most cases with no distinct limits. For example, in the Kinki MMA, Osaka is the city with the highest population and the Los Angeles metropolitan area includes the Long Beach and Santa Ana population. There are five MMAs included in this study, and they are shown in bold in Table A1. One last remark concerning population data is about Guangzhou, Shenzhen, Hong Kong, Dongguan and Foshan (shown in italics in Table A1). They are all located in the Guangdong–Hong Kong–Macau Greater Bay Area. This area, which is confined in a domain narrower than $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, also includes five other cities with population more than 1 million, forming a megalopolis, which is the largest urban agglomeration on the planet. The cities in this particular region and the MMAs will be exempted from the spatial gradient analysis.

Table 1 Number of cities according to their population (adopted by UN, 2018a).

Figure 1 Map of the location of the 81 cities with populations greater than 5 million (black circles). Red circles denote the 10 megacities with the highest population.

2.2 Spaceborne aerosol optical depth (AOD) and dust optical depth (DOD)

Daily retrievals of AOD at 550 nm from the MODerate resolution Imaging Spectroradiometer onboard the Aqua satellite (MODIS-Aqua) were used for this analysis. Specifically, quality assured Collection 6.1 MODIS – Aqua Level 2 AOD retrievals (Levy et al., 2013; Wei et al., 2019b) from the Dark Target (DT) and Deep Blue (DB) retrieval algorithms were merged according to Sayer et al. (2014) into one dataset and presented by Gkikas et al. (2021). AOD values from DT retrieval algorithm above ocean and land and DB above land were selected based on quality assurance (QA) of the retrievals. Then, the AOD values from both retrieval algorithms were merged by applying additional quality filtering criteria (see Sect. 2.1 in Gkikas et al., 2021) and stored in a scientific dataset. Based on this MODIS AOD combined product and on MERRA-2 (Randles et al., 2017; Gelaro et al., 2017) reanalysis data of dust optical depth (DOD) to AOD ratio, the MIDAS (ModIs Dust AeroSol; Gkikas et al., 2021, 2022) dataset was developed providing columnar mid-visible (550 nm) DOD and the corresponding AOD, on a daily basis, at fine spatial resolution ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$) and at global scale over the period 2003–2020. Regarding the DOD uncertainties, there were assessed in Gkikas et al. (2021) by considering the AOD and Merra-2 dust fraction uncertainties, as they were calculated using AERONET retrievals (Giles et al., 2019) and LIVAS database (Amiridis et al., 2015) as references, respectively.

2.3 Ground-based (GB) measurements

AERONET ground-based measurements of AOD for cities that have available long-term time series in the 2003–2020 period were also analyzed. The level 2, version 3 daily mean product of AOD at 500 nm was collected for stations with at least 8 years of data within the study period. Using the Ångström exponent, 440–870 nm, the AOD values were interpolated at 550 nm. Table B1 gives the details of the 27 AERONET stations utilized for the comparison. Statistics and trends of ground-based measurements were calculated with the same methodology and data availability criteria applied to the satellite data (see Sect. 2.5 and 2.6).

2.4 Spatial features

Initially, the geographical distribution of long-term (2003–2020) mean annual and seasonal AOD was derived for a square area spanning of $\pm \mathrm{1}{}^{\circ }$ around the city center (as defined by the UN 2018b database). The goal of this analysis was to investigate the AOD variability inside the urban areas, as retrieved by spaceborne data. In addition, visual inspection of aerosol distributions over $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ areas and the identified patterns was used for the first qualitative classification of megacities. The main advantage of the MIDAS AOD dataset is the high spatial resolution (of 0.1^∘) and the daily data availability. However, there were challenges during averaging regarding the data availability, and for this reason temporal availability criteria were applied. Every seasonal mean value was calculated when at least 9 d were available (1st threshold applied), whereas the annual means were computed when at least three seasonal means were valid (2nd threshold applied). The choice of the 9 d was decided after conducting a sensitivity analysis (not shown here) based on data availability limitations for calculating AOD trends. Based on the filtered year-by-year seasonal and annual means, their respective long-term averages were calculated. In Fig. 2, the geographical distributions of the annual (Fig. 2b) and seasonal (Fig. 2d) long-term averages are illustrated for the megacity of Tokyo. Comparing the AOD field $\pm \mathrm{1}{}^{\circ }$ around the city's center (black circle) with the corresponding Google Earth satellite map (Fig. 2a), higher values of mean AOD can been seen over the urban agglomeration of Tokyo, with respect to the surrounding area, throughout the year (Fig. 2b, d). The blank boxes are pixels that do not fulfill the applied criteria for data availability or data are not available. This example case demonstrates the applied methodology as it will be explained in the following sections.

Figure 2 (a) Geographical limits of the study area for the megacity of Tokyo on © Google Earth maps, © Google Maps. (b) Geographical distribution of the long-term (2003–2020) annual averaged AOD for the broader area of Tokyo megacity. (c) Spatial representation equal lat/long grid of AOD data and the sectors under investigation with different colors. (d) Geographical distribution of the seasonally averaged AOD for the broader area of Tokyo megacity. Blank pixels do not fulfil the data availability criteria. The black circle denotes the pixel of the megacity center.

In order to investigate further the spatial AOD variability, we quantified the AOD gradients between the city's center and the surrounding areas. The square area at each city was divided into six different sectors (Fig. 2c). Sector 1 (or S1) is the 4-pixel area that encompasses the city center. The daily AOD time series has been constructed by calculating the AOD median for each sector under investigation. Median was selected rather than mean for AOD aggregation, as a non-parametric statistic of central tendency, accounting for the not normally distributed AOD data at those spatial and temporal scales (Sayer and Knobelspiesse, 2019). This aggregation process increased the data availability compared to the single pixel approach. Following the previous procedure with filtered seasonal and annual means, the long-term average AOD for every sector was derived. An example for Tokyo is given in Fig. 3. For this megacity, a mean AOD value ∼0.35 was found over the city's center and decreasing mean AOD values moving towards the outer sections. Although a uniform approach like this ignores the effect of topography (mainly mountains and sea) that breaks the symmetry around many cities, it could be considered an indicator of the spatial distribution of AOD. It should also be highlighted that the “city center” area could have different characteristics and be wider than the area assumed in this approach. These assumptions are used only to make the results comparable, and studies focusing on a small number of cities should make a specific division for each case. In order to give a single number that will describe the spatial gradient of AOD field from megacities center to the surrounding area, linear regression was performed to the sectoral annual AOD averages in order to calculate the AOD changes per 0.1^∘ along with their statistical significance. The results then were expressed in Δ(AOD) per 1^∘ and finally were converted in Δ(AOD) per 100 km (see next paragraph) and were used for a categorization of the megacities.

Figure 3 The long-term average of AOD (red points) for Tokyo megacity for the six different sectors considered around every megacity center. The shaded area denotes 1σ. The solid red line is the regression line resulted from the six points.

Another source of uncertainty for the unified pixel-based approach that was followed for all cities would be the east–west (E–W) direction size of the pixel. While the pixel size in the north–south (N–S) direction is the same (∼110 km for 10 pixels $=\mathrm{1}{}^{\circ }$), this is not the case in the E–W direction. More than 50 % of the cities are located at absolute latitudes 15–45^∘ (left panel Fig. 4) with a median size of 10 pixels ($=\mathrm{1}{}^{\circ }$) in E–W direction of 95 km (central panel Fig. 4) and with majority of the cities (over 75 %) within the ±10 km limit. Even for the high latitude cities ($\sim \mathrm{60}{}^{\circ }$), the E–W size of the 10 pixels ($=\mathrm{1}{}^{\circ }$) is ∼65 km ($\sim \mathrm{2}/\mathrm{3}$ of the median). This difference for the high latitudes cities has been considered in the subsequent two-domain analysis (see Sect. 2.5), by including additional pixels in the analysis. For the six-sector analysis, we expressed the single measure of the AOD gradients in terms of distance in order to account for the difference in the east–west direction. We derived the equivalent distance of the 1^∘ for every city by taking the mean of 110 km size in the N–S direction and the distance in E–W direction (central panel Fig. 4) and subsequently we converted the single measure of the AOD gradients from Δ(AOD) per 1^∘ in terms of distance in Δ(AOD) per 100 km. Finally, the right panel of Fig. 4 shows the distribution of cities elevation, a parameter that has been also considered in the interpretation of the results.

Figure 4 The distribution of the absolute values of megacities latitude (a). The distribution of the equivalent length of 1^∘ east–west (E–W) direction in km (b). The distribution of megacities elevation in m (c).

2.5 Temporal variability and long-term averages

In order to investigate the temporal variability of AOD and derive the long-term mean values, an area of $\mathrm{0.4}{}^{\circ }\times \mathrm{0.4}{}^{\circ }$ (4×4 pixels, S1 and S2 combined) was considered as representative of the urban agglomerations, denoted as urban domain hereafter (Fig. 5). In the same manner as in six-sector analysis (Sect. 2.4), the daily AOD time series has been constructed for the urban domain by calculating the AOD median. Filtered (the same temporal criteria as in Sect. 2.4) annual mean AOD values were calculated and, based on them, the long-term (2003–2020) annual mean values were derived for the urban domain. The AOD interannual variability and decadal trends for the urban domain were derived, and the methodological details are described in the following subsection. A 4×6 pixels area (instead of 4×4) was considered for high latitude cities (London, Moscow and Saint Petersburg) in order to have comparable areas for all cities. The differences in annual mean AOD were 1 %–5 % compared to the same results without making this correction. Despite the small increase in data availability (∼10 %) when this correction was applied, the trend calculations remained stable. In general, the uncertainties introduced due to the pixel size are minor compared to those associated with the data availability.

Figure 5 Spatial representation of equal lat/long grid of the AOD data and the urban domain (red) and surrounding area domain (brown) under investigation.

2.6 AOD interannual variability and trends

The linear trends were calculated by using simple linear regression for the filtered annual mean values when at least 10 years (out of 18) were available. The statistical significance of the trends was assessed by performing the t-test. In Fig. 6, we are presenting an example for Tokyo. The red and brown curves correspond to the annual means for the urban and surrounding domains, respectively, the shaded areas represent the standard deviation, and the same colored lines denote the calculated linear trends. For this interannual time series, declining AOD trends are evident both for the urban and the surrounding domains of the Tokyo megacity.

Figure 6 Annual mean AOD values for the Tokyo megacity (red points) and the surrounding area (brown points) and the corresponding standard deviations (shaded areas). The linear trends are depicted with the straight lines.

2.7 AOD intraannual variability

Monthly mean AOD values were calculated when 3 d were available. From those filtered monthly AODs, the long-term monthly mean AOD values were derived, in order to assess the intraannual variability of aerosol loads for every city's urban domain. As an example, the intraannual variability for the megacity of Tokyo is given in Fig. 7. In order to give a single measure of this intraannual variability of AOD, the temporal coefficient of variation (CV) was derived. The mean and the standard deviation (SD) were calculated for cities with at least nine monthly values in order to calculate the temporal CV using the following formula:

$$\begin{array}{}\text{(1)}& \mathrm{CV}={\displaystyle \frac{\mathrm{SD}}{\mathrm{mean}}}\mathrm{100}\mathit{\%}.\end{array}$$

Figure 7 (a) AOD intraannual variability for Tokyo megacity (red line) and the surrounding area (brown line). (b) Normalized monthly mean AOD values with respect to June mean AOD value.

2.8 Spatiotemporal variability

In order to investigate the spatiotemporal AOD variability, the long-term mean AOD values and the linear trends were also derived for a second domain corresponding to the surrounding area of the urban agglomerations (brownish cells in Fig. 5). This domain (surrounding domain hereafter, Fig. 5) is the remainder from a broader 8×8 pixel domain around the city's center after eliminating the urban domain (4×4 pixels). Again, for high latitude cities, extra pixels were considered (12×8 pixels) adjusting for different pixel area, but even smaller differences were observed (1 %–4 %), with respect to those found for the urban domain considering this change.

3 Results and discussion

3.1 Megacities' spatial AOD characteristics

In this section, megacities' AOD spatial variability has been investigated and used as a proxy for cities' classification.

3.1.1 Geographical distribution

As already mentioned in Sect. 2.4, the aim of deriving the geographical distribution of the 18-year period mean annual and seasonal AOD at the original spatial resolution (0.1^∘) of satellite retrievals over $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ area around every cities center was to investigate the AOD variability inside and around cities by the visual inspection of those geographical distributions. An example of the analysis performed is given for Tokyo megacity in Sect. 2.4. The same approach was followed for all the examined cities (not shown). The first important result of this analysis was that the comparison of the AOD fields ($\pm \mathrm{1}{}^{\circ }$ around cities' centers) with the corresponding Google Earth maps confirms that the selection of the 4×4 urban domain is representative of cities' urban agglomerations and the remainder up to 8×8 pixel area representative of cities' surrounding areas. Additionally, by examining the spatial AOD heterogeneities between the urban and surrounding domains for all megacities, interesting features were found, which led to their first qualitative classification.

Out of the 81 cities, 53 % are coastal and 47 % are inland. This classification is important, as coastal cities are related with low data availability – due to satellite retrieval algorithm restrictions – and most cities that are not present in the following results are from this group because they did not fulfill the criteria of data availability. One group consists of high-elevated (>2000 m) cities such as Mexico City and Bogota. Those cities are related with low data availability, and for the available pixels the AOD values are very low, which is not consistent with the GB measurements (see Sect. 3.5), so they were excluded from the spatial gradient analysis. Another group contains inland cities situated nearby large deserts where dust is regularly dominating the aerosol mixture (Cairo, Tehran, Riyadh, Baghdad and Khartoum).

Cities with few blank pixels (<50 %) within the defined limits are clustered in three distinct groups. There are cities where the AOD fades from the central sector towards the surrounding area (∼60 % for coastal and ∼35 % for inland). This is always the case for the MMA group of megacities. For a number of cities (∼40 % for coastal and for inland), an opposite gradient is found, attributed to the higher anthropogenic activities (e.g., industrial zones) in the outer domain. Finally, in the third category are grouped only inland cities (∼25 %) where a uniform distribution of mean AOD exhibits a weak variation in spatial terms. This is a qualitative classification and it is based on the more uncertain results (based on the data at the original spatial resolution at 0.1^∘). In order to give a quantitative measure of this spatial variability and with lower uncertainty compared to the singe pixel approach, the six-sector and two-domain analysis was performed (see Sect. 3.1.2).

3.1.2 Long-term means and spatial variability

One of the objectives of this analysis was to take advantage of the high spatial resolution (0.1^∘) of MODIS AOD product in order to assess inequities in aerosol air pollution exposure levels within those high-risk communities. The aggregation of the AOD for the different sectors and domains that have been applied (see Sect. 2.4 and 2.5) increase the confidence in the results relative to the single pixel approach. It must be pointed out here that the combination of the local topography and the location of the anthropogenic activity outside the examined cities (industrial zones, nearby smaller cities) distorts the spatial AOD distribution, highlighting the possible non-homogeneity of the city sectors or domain of the surrounding area we have used.

Spatial (six-sector) analysis

In order to give a single measure of the AOD gradients (i.e., sharp or smooth) in the vicinity of the megacities, a linear fit was applied to the six sectoral 18-year mean AOD values. The slopes are expressed as Δ(AOD) per 100 km and along with their statistical significance are presented in Fig. 8. Negative values indicate that greater values of AOD were found in the city center and are denoted with blue color. Positive values are related with higher AOD values moving away from cities center and are denoted with red color. Due to low data availability, 24 cities were discarded from the analysis. The MMAs and cities with high elevation (>2000 m) were also excluded from this spatial gradient analysis. In Appendix C for every city, the mean AOD values for six sectors are given, normalized with the 18-year mean value of S1 (Fig. C1).

For the majority of megacities (∼65 %), lower mean AOD values resulted for the surrounding areas with respect to the city's center (i.e., negative Δ(AOD)). This percentage is raised to ∼75 % when only the statistically significant results are considered. The largest values of negative AOD gradients were found for Xian (China), Alexandria (Egypt) and Santiago (Chile), attributed to complex topography (high mountains and coastal areas). On the other hand, only 19 (∼35 %) cities have larger mean AOD values in the surrounding sectors (i.e., positive Δ(AOD)) and most of them are located in China and India. The reason that higher AODs were recorded in the surrounding sectors for most of those cities is that they are enclosing smaller-scale cities (MMAs have been already excluded from the analysis), which, however, host significant sources of anthropogenic aerosols. For example, Shanghai is the Chinese city with the largest positive AOD gradient attributed to the influence from the nearby megacity Suzhou at the west and the Pacific Ocean at the east part of the city.

Figure 8 Spatial AOD gradients per 100 km (Δ(AOD) per 100 km) from city center to surroundings. The circle radius is proportional to the statistical significance of the trends. Three cities with greatest negative AOD gradient are Xian, Alexandria and Santiago and the three cities with greatest positive AOD gradient are Shanghai, Jinan and Atlanta.

Two-domain spatial analysis

The geographical distribution analysis (Sect. 3.1.1) revealed that the 4×4 pixels area is a domain well fitted to the cities' urban agglomeration (urban domain). The remainder between the urban domain and an area up to 8×8 pixels (surrounding domain) was found also to be representative for the surrounding area of all cities. The results for this two-domain analysis are presented in this section.

The geographical distribution of the long-term mean AOD values for the urban domain is presented in Fig. 9. Four cities (BOGO, HCMC, KUAL, SING) are not included due to low data availability. Low mean AOD values were found for European and American cities in contrast to Asian and African cities. Specifically, the larger mean AOD values (>0.5) were recorded for Chinese cities, with Indian cities following, since both areas are densely populated and with high industrial activity. Many megacities (mostly Asian) lying in the proximity of great deserts are also influenced by natural aerosols of desert dust (Proestakis et al., 2018; Gkikas et al., 2022). In order to investigate further this feature and quantify this influence, the long-term mean DOD and the DOD to AOD ratio – as a metric of dust contribution in the aerosol mixture – were derived using the MIDAS DOD and AOD products and the obtained results are presented in Fig. 10. Significant dust contribution (∼20 %–40 %) was found for the western Chinese cities, with mean DOD values ranging from 0.1 to 0.2, that were influenced by the nearby arid and semi-arid regions (namely the great Gobi Desert and the Taklamakan Desert). For Indian and Pakistan cities around the Thar Desert (India–Pakistan borders), even bigger dust contribution was found ∼40 %–50 % and mean DOD values ranging from 0.15 to 0.35. The greatest dust contribution (up to 60 %) with significant mean DOD values up to ∼0.4 was found for the African cities of Khartoum and Cairo and for all Middle East cities (Riyadh, Baghdad and Tehran) that are influenced by the great deserts of North Africa and Middle East.

Figure 9 Long-term (2003–2020) mean AOD for megacities urban domain ($\sim \mathrm{40}\mathrm{km}\times \mathrm{40}\mathrm{km}$ around cities center). The circle radius is in proportion to cities' population. Three cities with greatest mean AOD are those of Chengdu, Foshan and Chongqing and the three cities with lowest mean AOD values are those of Atlanta, Houston and Belo Horizonte.

Figure 10 Long-term (2003–2020) mean DOD for megacities urban domain ($\sim \mathrm{40}\mathrm{km}\times \mathrm{40}\mathrm{km}$ around cities centre). The circle radius is in proportion to the DOD to AOD ratio. Five cities with greatest dust contribution (mean DOD $/$ AOD ratio) are Khartoum, Riyadh, Baghdad, Cairo and Tehran.

In order to contrast the mean AODs for the urban domain against those of the surrounding area, we have reproduced a global scatterplot with their matchups, presenting also ancillary information (i.e., continent, population, coastal/inland) (Fig. 11). Points residing on top of the one-by-one line indicate cities with homogeneous spatial AOD distributions (inside and surrounding area), whereas above/below the equality line, AODs are higher in the surrounding and urban domain, respectively.

Almost 60 % of the cities have greater mean AOD values over the urban domain and 40 % have greater values over the surrounding domain. There are more cities with differences exceeding 10 % when higher values of AOD are recorded over the urban domain compared to the opposite case. In general, the obtained results are in line with those of spatial gradients (sectoral analysis, Fig. 8), while almost all cities (apart from four) are included here due to the greater data availability for the two-domain approach. Shenzhen, Osaka and Jakarta show the biggest decrease of 0.15–0.25 (∼25 %–35 %) in a range of few kilometers from their center. On the other hand, Chinese cities (Shanghai, Suzhou and Qingdao) and Atlanta have greater mean AOD values in the area surrounding the center.

Additionally, the classification of the cities according to their geographical location (in line with Fig. 9) came up, according to the obtained results for both domains. All European and American cities yield mean AODs ranging from 0.08 to 0.20, in contrast to African and Asian cities in which the corresponding levels range from 0.25 to 0.9 (apart from Bangalore (0.2), Dar es Salaam (0.17) and Johannesburg (0.09)).

Figure 11 Long-term (2003–2020) mean AOD for megacities' urban domain ($\sim \mathrm{40}\mathrm{km}\times \mathrm{40}\mathrm{km}$ around cities' centers) versus the surrounding area domain. The name of the cities with difference in mean AOD greater than ±0.04 between the two domains is depicted. The dotted black line is the identity line and the dashed red lines denote ±10 % difference in mean AOD between the two domains.

3.2 Temporal variability

3.2.1 AOD trends

The geographical distribution of the AOD changes per decade (Fig. 12, after excluding 11 cities which did not fulfill the temporal criteria) revealed pronounced regional features. The resulted AOD trends for the megacities are reflecting the changes in the anthropogenic emissions, associated with the air quality regulations implemented throughout the years, in the same manner that previous studies have shown the connection of the satellite-observed AOD trends with the implemented mitigation policies in regional scales (Gupta et al., 2022; Zhao et al., 2017)

For all European and US/Canadian cities, decreasing AOD values were found (up to ∼0.03 and ∼0.07 per decade, respectively), in accordance with the AOD decrease in western Europe and eastern North America that have already been reported in the literature, and was associated with a series of air quality control measures that have been implemented (Gupta et al., 2022; Zhao et al., 2017). While negative trends have been found in this study for Los Angeles, other recent studies report small positive trends for the western United States of America that are associated with reduced precipitation and increased fire activity over the area (Gupta et al., 2022; Cherian and Quaas, 2020). However, the latter studies refer to much wider areas and not to the city of Los Angeles, which indicates that the trends in the city may be related to different mechanisms (i.e., reduction of anthropogenic emissions dominates over the increase due to increasing dust and smoke events), relative to the trends over the wider region of the western US. Statistically significant negative AOD trends were derived for the eastern Asian megacities, with the highest values (in absolute terms) up to ∼0.3 per decade being evident for the Chinese megacities. This result (net negative AOD trend for the whole study period for Chinese cities) is in agreement with recent studies that reported AOD decrease for eastern China (Gupta et al., 2022), which was associated with the implemented emission control policies. Specifically, for China, while up to 2010 AOD was increasing (e.g., Hsu et al., 2012) due to the rapid economic and industrial development of the country, after 2011 declining AOD trends have been recorded for the central and eastern sectors of the country, related with the reduction in anthropogenic aerosol emissions due to the implementation of emission control measures (Zhao et al., 2017; Sogacheva et al., 2018). According to Sogacheva et al. (2018), the gradual AOD decline after 2011 is more prominent for the highly populated and industrialized southeast China regions, being in agreement with our results.

On the contrary, strong positive AOD trends (ranging from 0.07 up to 0.16 per decade) were found for all megacities in the Indian subcontinent, reflecting the increased industrial and financial development during this period, and is in agreement with previous studies (e.g., Buchholz et al., 2021). In a recent study by Samset et al. (2019), they have shown that the climate implications of this dipole pattern of positive AOD trends over southern Asia and strong negative values over eastern Asia observed since 2010 might be strong not only on a regional scale but also for areas away from the sources.

Positive AOD trends were also found for the Middle East megacities, ranging from 0.03 to 0.1 per decade. This finding is in line with the AOD trends that resulted from the analyses of retrievals from different satellite sensors (Che et al., 2019), although negative trends have been reported in the study of Gupta et al. (2022) who analyzed AOD from CALIOP. For the Latin American megacities, the highest AOD trends are positive (∼0.03 per decade for Buenos Aires and Lima), reflecting the poor emission control measures over this geographical domain. In recent studies (Gupta et al., 2022; Buchholz et al., 2021), a reduction in AOD was reported for South America, attributed to the decline in forest fires, which might explain the negative AOD trends for the rest of Latin American cities. Sub-Saharan Africa is an area for which AOD decrease has been reported in previous studies, which is opposite to our findings (+0.06 per decade for Luanda and +0.03 per decade for Dar es Salaam). For this area, an increase in AOD was also found by the recent study of Gupta et al. (2022), which was associated with the increase in biomass burning of agricultural activity in the dry season over the area. However, since our study is focused on megacities, our findings may reflect the increasing urbanization in combination with the limited air quality regulations over the area, but further investigation is needed here.

Figure 12 Linear trend of AOD per decade for urban domain ($\sim \mathrm{40}\mathrm{km}\times \mathrm{40}\mathrm{km}$ around cities' centers). The circle radius is in proportion to the statistical significance of the trends. Three cities with greatest increase of AOD are Hyderabad, Kolkata and Bangalore and the three cities with greatest decrease of AOD are Chengdu, Chongqing and Xian.

The comparison of the AOD trends against the long-term mean AOD of the urban domain revealed an interesting clustering of cities, according to their geographical location (Fig. 13). It has to be pointed out here that only the statistically significant trends were included in Fig. 13, thus ensuring robust and meaningful results.

In all European and North American cities, low AOD climatological values and decreasing trends were found. Note that in some US cities like Atlanta, Washington, Philadelphia and New York, apart from their relatively low mean AOD values (<0.17), considerable negative trends (ranging from −53 % to −28 %) were recorded.

All Indian cities have results with positive AOD trends regardless of their mean AOD levels, which span from low to high values. Among them, Kolkata is the city where extremely high mean AODs (0.70) and large positive trends (+22 % per decade and the higher in absolute values $+\sim \mathrm{0.16}$ per decade) were revealed. Of particular interest is Bangalore, which has relatively low mean AOD value (0.20), but the maximum positive AOD trend is almost +69 % per decade. Bangalore's population increased from 6 million to 12 million during the last two decades, which is one of the biggest population increases for cities at this scale and thus it can explain to a large degree the extreme increase of AOD. Meanwhile, financial development in the area is linked more with new technologies and not heavy industry, which partially answers the relatively low mean.

All Chinese cities are subjected to AOD decrease (up to ∼30 %), while at the same time they have the highest mean AOD values ranging from 0.45 to 0.91. Chengdu is the city with the biggest mean AOD value along with the highest AOD decrease of ∼30 % per decade (or ∼0.3 per decade).

Figure 13 AOD decadal changes expressed in percentages (Δ(AOD) % per decade) versus the long-term mean AOD for the urban domain. Only statistically significant AOD trends are presented (P<0.1). Cities with absolute values of Δ(AOD) % per decade greater than 25 % are denoted in the figure.

3.2.2 Intraannual variability

AOD exhibits strong intraannual variability (Zhao et al., 2018), which is quite important to be analyzed at megacities' scales in order to improve our understanding regarding the aerosol-related health effects. The normalized (using the mean AOD of June as reference ) monthly mean values for all cities are provided in the Appendix D, Fig. D1, classified for the different geographical domains. In order to investigate the intraannual variability of AOD at city level, the temporal CV (Sect. 2.7) was calculated. The geographical distribution of temporal CV is illustrated in Fig. 14, only for cities complying with the defined temporal criteria (four cities were omitted). Low CV values (<20 %) were found for eastern India, most of the Chinese cities and for the eastern Mediterranean as well. The highest CV values >70 % were recorded for three US cities. In order to investigate those results further, the relationship between temporal CV and the long-term mean AOD values for the urban domain was examined (Fig. 15).

A non-linear relationship was found, with the gradual decrease of the mean AOD to be related with increasing CV levels. The CV values can reach at their maximum levels (>50 %) only in the group of cities with low mean AOD values (<0.25, with only exception being Luanda) and all are American cities on the East coasts. For North American cities, the monthly mean AOD values during the boreal winter are decreasing close to zero and maximized during summer, this resulted in this high intraannual variability (Fig. D1b) especially for southeastern US cities. These high values of intraannual variability for southeastern US cities, and most notably in Atlanta, are explained by the high AOD values over the area during summer due to secondary aerosol formation by biogenic volatile organic compound (BVOC) oxidation (Goldstein et al., 2009). For most southern American cities, rather stable values were found throughout the year (Fig. D1c). For moderate AOD levels (0.25–0.50), a limit of 50 % was found for the CV values. For the cities with high aerosol loads (>0.50), which are mostly Chinese cities, the CV values are confined to the limit 10 %–40 %.

Figure 14 Temporal coefficient of variation (CV) of monthly mean AOD (intraannual variability) expressed in percentages for megacities' urban domains ($\sim \mathrm{40}\mathrm{km}\times \mathrm{40}\mathrm{km}$ around cities' center). The circle radius is proportional to the cities' population. Three cities with the highest temporal CV are Atlanta, Houston and Washington, DC and the three cities with lowest values are Kolkata, Chongqing and Pune.

Figure 15 Temporal coefficient of variation expressed in percentages (CV %) against the long-term mean AOD for the urban domain. The size and the color of the points denote cities' population and geographical domain, respectively. The shape of the points denotes costal or inland cities.

3.3 Spatiotemporal variability

Towards assessing the AOD spatiotemporal variability, the decadal AOD linear trends for the urban and the surrounding domains were compared (Fig. 16). We are presenting only the points in the scatterplot when the AOD trends are statistically significant (p<0.1), both for the surrounding and urban domains (43 pairs). Points residing above/below the identity line correspond to AOD trend (negative or positive) that is higher in the surrounding/urban domain, respectively.

All cities in India exhibit an increase of AOD both in the urban and surrounding domains. The largest difference was found in Hyderabad where AOD trend is 10 % larger for the urban domain, indicating an increase of the anthropogenic activity within the boundaries of the urban agglomeration. For cities where AOD decreases were found, two groups are shaped. The first one, that has greater AOD decrease for the urban domain, consists mostly of East Asian cities. Chongqing is the Chinese megacity with the largest difference (∼20 %), reflecting the adaptation of effective measures towards reducing city's aerosols air pollution (Gao et al., 2021). The second group, which consists mostly of the North American cities, has greater AOD decrease for the surrounding domain.

Figure 16 AOD decadal changes expressed in percentages (Δ(AOD) % per decade) for the surrounding domain versus the urban domain. Only statistically significant AOD trends are presented (P<0.1). The dotted line is the identity line. Cities with absolute differences of Δ(AOD) % per decade greater than 10 % between the two domains are denoted in the figure.

3.4 Evaluation with GB measurements (AERONET)

For cities with available long-term time series of AOD from GB stations of the AERONET network, an evaluation of the satellite AOD averages and trends was performed. Table B1 summarizes the information of the AERONET stations utilized. This comparison is separated by the GB stations located within the 4×4 pixel area of the urban domain and those residing in the surrounding domain of a city. The daily averaged AOD product was used from the AERONET stations to derive long -term mean AOD and linear trends, using the same approach as followed with the satellite data. The diurnal AOD variability of the AERONET data was not considered and this may slightly affect the long-term mean AOD (Smirnov et al., 2002), but it plays a minor role for the calculation of trends.

Regarding the long-term mean AOD (Fig. 17), for the urban domain there is a good agreement between the satellite-derived and the GB values (correlation coefficient R∼0.93). The largest deviations (expressed in percentage terms) between spaceborne and ground-based AODs are recorded in general for weak-load cities, and those deviations are maximized in the high-altitude city (>2000 m) of Mexico and Osaka (limited satellite data availability). For the surrounding domain, the studied cases are few, so the R is not representative. The satellite derived mean AOD is overestimated, but in this case the area that was utilized to derive the satellite results is considerably large and it is unlikely that the point measurement statistics would coincide with satellite statistics.

Figure 17 (a) Satellite-derived (red line) and GB (blue line) long-term mean AODs and their relative difference expressed in percentage (black line) for the urban domain of cities. (b) The same for the surrounding domain (brown line for the satellite-derived values).

A relatively good agreement was also found for the decadal linear trends (Fig. 18), with R∼0.79 for the urban domain. For 75 % of cases, GB and satellite-derived trends have the same sign for both domains. In Fig. 18, the results are presented only when both spaceborne and ground-based AOD trends are statistically significant. For the urban domain, all the statistically significant trends have the same sign. For the surrounding domain, the trends' signs differ only for LOSA. This difference can be attributed to the shorter time period of measurements at the LOSA AERONET station (8 years). Additionally, in a recent study (Wei et al., 2019c), western North America was found to be the area with fewer sites with the same signs between MODIS- and AERONET-derived AOD trends. Nevertheless, our analysis does not focus on a validation of satellite AOD against AERONET, thus not an exact colocation of the AERONET stations with satellite pixels was attained. In order to further investigate the obtained differences, the spatial CV was calculated for the urban and surrounding domain, as the ratio of standard deviation and mean values of the pixel long-term mean AODs derived in geographical domain analysis (see Sect. 2.4). The results are presented in the lower panel of Fig. 18 along with the percentage of the available pixels inside the domain under investigation. For LOSA, the greatest value of spatial CV was found (>50 %), which was anticipated, since the spatial extent of Los Angeles metropolitan area is high (including Long Beach and Santa Ana), highlighting the importance of the selection of the GB location.

An additional interesting feature is that the absolute values of all satellite-derived statistically significant trends (Fig. 18) have lower magnitude compared to the GB results. This difference in magnitude could be attributed to the domain-aggregated satellite values compared to point-derived results of the GB stations, which may smooth out the AOD fields. In a recent study by Logothetis et al. (2021), a sensitivity analysis was performed between coarse (1^∘) and fine (0.1^∘) spatial resolution of spaceborne AOD retrievals which revealed AOD trends of lower magnitude for the coarse spatial resolution data. Moreover, the daily satellite value, corresponding to satellite overpass, contains less information compared to the continuous monitoring during daytime of cloudless days of the ground-based photometers.

Hence, GB stations are representing the close area around their location and the representativeness of AERONET stations to characterize aerosol load/trends in megacities should be considered individually in each case, according to local conditions.

Figure 18 (a) Satellite-derived (red line) and GB (blue line) statistically significant (P<0.1) AOD linear trends for the urban domain and the corresponding spatial CV (red bars) for every city; the numbers on top of the bars declare the pixel availability for the domain. (b) The same for the surrounding domain (brown line for the satellite-derived values).

4 Summary and conclusions

Motivated by the environmental challenges caused by increasing urbanization and to maximize the use of spaceborne aerosol products, in this study, we investigate the aerosol regime over the megacities of the world using satellite aerosol retrievals. We are taking advantage of the global coverage, the high sampling frequency (daily values) and the relatively fine spatial resolution ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ grid) of the 18 years MODIS-based AOD and DOD at 550 nm products, in order to examine the spatiotemporal variability of aerosol loads for the largest 81 cities of the world.

For all European and American cities, mean AOD ranges mainly from 0.08 to 0.20. For all African and Asian cities but three (Bangalore (0.2), Dar es Salaam (0.17) and Johannesburg (0.09)), mean AOD ranges from 0.25 up to almost 0.90. There are cities which lie in the proximity of deserts or in the path of transported mineral dust particles, which were found to have considerable dust contribution (up to 70 %) and the associated mean DOD values (up to 0.4) can be considered as a constant background aerosol source for those megacities. Results of this contribution can be used by policymakers for defining the legislations on air quality urban thresholds.

The majority of cities (∼60 %) have greater mean AOD values over their urban agglomeration domain. Mostly Chinese and Indian megacities tend to have higher AOD in the surrounding areas of the city center. Shanghai is the city with the largest difference (13 % greater mean AOD values for the surrounding domain), but in general for the cities grouped in this category, the declinations between the two domains are lower compared to the first category, because of the high mean AOD values of the urban domain. Finding inequities in the exposure at urban scales using satellite remote sensing, may be a useful tool for air pollution assessment and finally taking diverse reduction measures at community level. As AOD differences are also observed for different months/seasons during the year, AOD intraannual variability at megacities' scale have been quantified. Low intraannual variability (temporal CV 10 %–40 %) was found for Chinese, eastern Indian and eastern Mediterranean megacities, while high values (>50 %) of temporal CV was recorded for eastern American megacities.

Although Chinese cities were found with the highest mean AOD values (up to ∼0.90), they also exhibit the highest AOD decrease in absolute values up to ∼0.3 per decade (or 30 % per decade), in response to the rigorous emission control measures implemented in the country, especially after 2010. The effectiveness of those measures also reflects the fact that, for Chinese cities, the AOD decrease was found to be higher for the urban domain (up to ∼20 %). Decreasing AOD values were also derived for US/Canadian and European megacities (up to ∼0.07 and ∼0.03 per decade) due to a series of air pollution control policies in the last decades. The maximum worldwide negative AOD trends in relative terms, ∼30 %–55 %, per decade were derived for most US cities, which simultaneously have low mean AOD values (less than ∼0.15). The highest AOD increase worldwide in absolute (up to ∼0.16 per decade for Kolkata) and relative terms ($\sim +\mathrm{70}\mathit{\%}$ per decade for Bangalore) was found for Indian megacities. Statistically significant positive AOD trends were found for all Indian cities, reflecting the increasing urbanization and industrialization of the country. The AOD increase for Indian megacities was found to be greater (up to ∼10 %) for the urban domain. Statistically significant AOD increase was also found for Middle East, South African and some Latin America megacities (up to ∼0.1, 0.06 and 0.03, respectively).

For cities where long-term ground-based AERONET measurements of AOD were available, the extent at which those measurements can capture the spatial and temporal AOD variability was investigated with respect to the spatial and temporal variability that were derived from satellite data. For GB stations within the urban agglomerations, a good agreement of the long-term mean AOD was found (R∼0.93) and with coincident signs on AOD trends in 75 % of the selected stations. The resulted discrepancies are attributed, apart from the satellite retrievals-related limitations (one overpass per day, high elevated pixels, low data availability, etc.) and the GB retrieval limitations (large temporal gaps due to instrument issues/calibration, shorter operating time periods, etc.), to the point (GB) versus area-averaged (satellite) comparison. It was found that for areas with non-homogeneous aerosol fields (great spatial CV, e.g., LOSA in our analysis), great differences (opposite signs) were recorded between satellite and GB trends. Those findings highlight the importance of the GB station location selection for future planning of aerosol measuring sites, to achieve representative AOD measurements for a specific city.

We acknowledge that the total column optical property of AOD, analyzed in this study to describe aerosol load variability, is not always proportional to the surface particulate matter (PM) concentrations, which is a parameter better describing the cities' air quality, is monitored and regulated by the cities' authorities and is related directly to health effects. However, PM concentrations are derived by in situ measurements and a dense network of those measurements is needed to describe the PM distribution of a city. Although these ground air quality monitoring stations are valuable, there are relatively few and unevenly distributed networks within cities around the world, especially in developing countries. One way to provide consistent PM data worldwide is by combining available PM measurements with satellite AOD observations and chemical transport models (e.g., HEI, 2020). Therefore, the use of satellite AOD is a very fitting source of information in order to have global coverage and high spatiotemporal resolution of aerosol loads over urban areas, keeping in mind that the agreement between satellite columnar AOD and ground-based PM concentrations is determined at a large degree by the vertical structure of aerosol layers (e.g., Gkikas et al., 2016).

According to our findings, long-term high-resolution spaceborne AOD retrievals can be utilized for detecting spatial and temporal aerosol variability at an urban scale, helping towards the current and future assessments of aerosol-related impacts in megacities. The high resolution of MODIS aerosol product provides access to cities' aerosol monitoring information which can serve as the basis for health-related and other prediction services. Future work linking the AOD changes with population and emission trends in these cities will reveal the linkage and will enhance the performance of the air quality projections for the next decades.

Appendix A

Table A1 Analytical table of 81 cities with the highest population up to 2018 (adopted by UN, 2018a, 2019a). The abbreviations of the statistical concept column stand for city proper (CP), urban agglomeration (UA) and metropolitan area (MA). Major metropolitan areas (MMAs) are shown in bold, and the cities located in the Guangdong–Hong Kong–Macau Greater Bay Area are shown in italics.

Appendix B

Table B1 AERONET stations that have been used in the analysis and the corresponding urban agglomerations.

Appendix C

Figure C1 Panels (a)–(g) correspond to different geographical domains. For every panel, the long-term mean AOD for the six different sectors normalized with the mean value of sector 1 (S1) is given with different colors, and symbols correspond to different cities. Solid/dashed lines denote inland/coastal cities, respectively. The bar plot shows the long-term mean AOD of S1 for every city in the geographical domain, denoted with the same color used in the upper plot. The dashed line corresponds to the mean from all cities, while the dotted line gives the mean AOD for the cities of the corresponding geographical domain.

Appendix D

Figure D1 Long-term monthly mean AOD values (intraannual variability) normalized against the mean value of month June. The different colors and symbols correspond to different cities. Solid/dashed lines denote inland/coastal cities, respectively.