Spatiotemporal Analysis of MODIS Aerosol Optical Depth Data in the Philippines from 2010 to 2020

Satellite remote sensing for air quality assessment provides information over a large spatial coverage and time period that shows the trends and effects of anthropogenic activities. Using data collected from the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra satellite from the years 2010 to 2020, the spatiotemporal variations to aerosol optical depth (AOD) in Koronadal City and Quezon City were studied. Validation showed a strong relationship between the MODIS AOD values and the Aerosol Robotic Network (AERONET) AOD values (R2 = 0.83) and a low root mean square error (RMSE) of 0.26. Annual variation in the AOD of the two study areas showed a peak AOD value in 2015 due to an immense biomass burning in Indonesia and a low AOD value in 2020 due to the COVID-19 lockdown. Koronadal City experienced a high AOD value during the fall season due to aerosols from biomass burning in Indonesia that were carried by the southwest monsoon. Quezon City experienced a high AOD value during spring from increased local sources, as well as long-range transport pollutants from East Asia that were carried by northeasterly winds. Overall, this study provides an understanding of the spatiotemporal variations in aerosols in the Philippines, which could be used in environmental management, air quality regulations, and health assessment studies. This shows the urgency of monitoring and mitigating poor air quality in the Philippines.


Introduction
Air pollution, including atmospheric aerosols, is the effect of not only local emissions but also foreign emissions from surrounding areas. The amount of air pollution transported depends on the country's topography and weather. Because of its growing industry, the anthropogenic emissions in Asia are greater than those in Europe and America. Residents of Asia are nine times more exposed to air pollution than their European and American counterparts. The primary source of these pollutants is emissions from the fuel combustion of motor vehicles [1]. Air pollution is currently considered as the largest environmental risk factor [2].
Atmospheric aerosols refer to suspended particles with a diameter of less than 10 µm. Aerosols can be in the form of natural sources, such as dust, sea salts, and volcanic ash, or can be from anthropogenic sources, such as the combustion of fossil fuels and biomass burning. Aerosol particles can affect the estimation of global and regional climate change as well as cloud formation [3][4][5] through their interaction with incoming solar radiation [6]. Some aerosols may pose risks to health and increase the likelihood of developing certain diseases, including pulmonary diseases, neonatal diseases, diabetes, and cerebrovascular diseases [7,8]. It is important to study the characteristics of atmospheric aerosols, as well as their spatial and temporal variations, so as to evaluate their roles in atmospheric dynamics and as health risks. Aerosols can be evaluated using different parameters: the aerosol optical depth (AOD), single scatter albedo (SSA), Ångström exponent (AE), and the complex refractive index. The AOD value represents the geographical distribution of aerosols or the column-integrated mass distribution. The AE is an optical property of an aerosol that measures the wavelength dependence in aerosol absorption, as well as a qualitative indicator of the size distribution of the aerosol [5,9]. There are several networks providing ground-based measurements of aerosol properties, such as the Aerosol Robotic Network (AERONET) [10], SKYNET [11], and Sun-Sky Radiometer Observation Network (SONET) [12]. However, there are still regions and areas with few studies and little to no monitoring [13,14]. Furthermore, the data measured at each site can only account for a small area. Therefore, the high demand for information about aerosol properties on a wide scale cannot be met at present.
In the last few years, the process of remote aerosol retrieval from satellite sensors has progressed. One of these advancements is the derivation of aerosol products from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard NASA's (National Aeronautics and Space Administration) Terra and Aqua satellites. Studies of and improvements to MODIS show that retrievals fall within the acceptable accuracy values [15]. There have been studies about aerosol trends in the Philippines [16][17][18][19][20][21]. Studies by [19,22] show that the Philippines, a maritime country, has the highest AOD concentration during October. This is due to the strong land and sea breeze created by warmer land surface temperatures. Long-range transport of anthropogenic emissions from East Asia and Indonesia also contributes to the elevated concentrations of aerosols [23][24][25].
This study aims to examine the annual and seasonal spatiotemporal variations in AOD values in the Philippines from 2010 to 2020, based on MODIS data at a spatial resolution of 10-by-10 array of 1 km pixels, and the V6.1 algorithm developed by NASA. The MODIS AOD values were validated with the ground monitoring network AERONET AOD using the CIMEL CE-318 photometer at the Manila Observatory and Notre Dame of Marbel University.

Related Literature
The Moderate Resolution Imaging Spectrometer (MODIS) is an instrument onboard NASA's Earth Observing System (EOS) satellites Terra, launched in December 1999, and Aqua, launched in May 2002. The Terra and Aqua satellites detect a wide spectral range of electromagnetic energy. The MODIS can acquire high radiometric-sensitive images (12 bit) in 36 spectral bands between 0.62 and 14.385 µm with a scanning width of 2330 km. The spatial resolutions of the sensor are 250 m, 500 m, and 1 km in several bands and can sweep the entire Earth surface every 1 to 2 days [26]. While Terra has a morning orbit, with an equator crossing time of 10:30 a.m., Aqua's is in the afternoon, with an equator crossing time of 1:30 p.m. [27]. MODIS can evaluate the AOD values over land and ocean using the Dark Target (DT) and Deep Blue (DB) algorithms. The DT algorithm measures the AOD in dark vegetative surfaces [15,[28][29][30][31] and over the ocean [15,32,33], while the DB algorithm measures the AOD in desert and vegetated surfaces [34,35]. Recent developments to the DB algorithm allow it to include other land areas [36][37][38].
As for ground monitoring networks, the Aerosol Robotic Network (AERONET) is a network of ground-based sun photometers that measure atmospheric aerosol properties. The AERONET employs a CIMEL Electronique 318A spectral radiometer that measures the sun and sky radiances obtained at a number of fixed wavelengths within the visible and near-infrared (VNIR) spectrum [10]. AERONET datasets are available in three quality levels: Level 1.0 (unscreened), Level 1.5 (cloud screened and quality controlled), and Level 2.0 (quality assured) [39]. The AERONET AOD uncertainty is 0.01 for larger wavelengths and 0.02 for shorter wavelengths [9]. Currently, there are approximately 400 stations in 50 countries, with 3 in the Philippines, namely, Manila Observatory in Quezon City, El Nido Airport in El Nido, and Notre Dame of Marbel University in Koronadal City. However, El Nido Airport stopped reporting ground measurements after 2013.

Study Area
The study area included two locations in the Philippines: Quezon City (latitude: 14.6 • N; longitude 121.1 • E; elevation: 63.0 m) and Koronadal City (latitude: 6.5 • N; longitude: 124.8 • E; elevation: 70.0 m). El Nido was not included because the ground station's limited data are not sufficient for the decade-long study. Figure 1 shows the location of the study areas, along with the corresponding AERONET ground monitoring stations. The Philippines is a tropical country and belongs to the Western North Pacific (WNP) boreal summer monsoon region. The wet season is experienced from May to October, while the dry season is experienced the rest of the year. Quezon City (QC) is located in the northeast part of Metropolitan Manila, on Luzon Island. It has a total area of 151.06 km 2 [40] and a population of 2,960,048 [41]. The city is situated on the Guadalupe Plateau and is largely rolling, with alternating ridges and lowlands. The lower part of the city has low-grade terrain, while the northern half is undulating and culminates at the Novaliches Reservoir [42]. Koronadal City (KC), on the other hand, is located on Mindanao Island and has a total land area of 277.0 km 2 . KC is surrounded by the mountain ranges of Roxas and Quezon. In general, the city has gently sloping terrain and 50% of the total land area is predominantly flat [43]. The latest Philippine Statistics Authority (PSA) census puts the total population of the city at 195,398 [41]. QC has an urbanization level of 100%, while KC has an urbanization level of 75.8% [44].

Dataset
The MODIS AOD data that were downloaded from the Terra satellite are reported in Level 2 (L2) granule-based and Level 3 (L3) global-gridded products. Daily MODIS AOD Level 2 products (MOD04_L2) at 10 km swaths from 2010 to 2020 were used for the analysis of AOD trends. In this study, the Deep Blue algorithm was utilized. The DB algorithm retrieves AOD data at 550 nm over global land in cloud-free and snow/ice-free scenes. Quality flags are also produced with each retrieval to ensure quality of retrieved AOD product.
To validate MODIS AOD data, this study used AERONET AOD data at Level 2.0, which are quality assured, meaning that they are cloud-screened and quality-controlled [10]. The study used the data from the two AERONET stations in the Philippines. These are located within the study areas, namely, Manila Observatory in QC and Notre Dame of Marbel University in KC.

Pre-Processing
The MODIS Terra reports AOD values at 550 nm; therefore, the following sub-dataset was extracted and then projected into 10 km grids: Deep_Blue_Aerosol_Optical_Depth_550_Land.
The collected AERONET AOD data included the following: the acquisition time, AOD value at 500 nm, and Ångström exponent at 440 nm to 670 nm. Since the MODIS and AERONET AOD values are reported in different wavelengths, interpolation for AERONET at 550 nm was calculated using Equation (1): where λ 1 and λ 2 are the wavelengths of the MODIS and AERONET, respectively; τ 1 is the calculated AOD value for AERONET at 550 nm; τ 2 is the collected AERONET AOD value at 500 nm; and α is the Ångström exponent for 440 nm to 670 nm from AERONET [9]. AERONET AOD measurements within ± 30 min of MODIS overpass time were considered.

Statistical Analysis
To extract the MODIS AOD values, a 3-by-3 pixel array centered at the study area ground site was set up. A condition requiring a minimum of 6 valid pixels was imposed on the extracted MODIS AOD data to minimize possible gaps. The mean of the pixel values was then compared with the corresponding AERONET AOD.
The data were analyzed using the following criteria: the arithmetic mean (x), standard deviation (σ), coefficient of determination (R 2 ), and root mean square error (RMSE).
where x and y are the arithmetic means of the MODIS and AERONET AOD, and x t and y t are the AOD values at time t, with n number of samples. An error of estimation was also presented using the following equation [45]:

Air Mass Back-Trajectory
Analysis of the wind back-trajectories was undertaken using the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model [46]. Meteorological conditions were determined using the NCEP/NCAR Reanalysis meteorological datasets [47]. Trajectories were calculated based on three altitude endpoints: 100 m, 500 m, and 1000 m. One trajectory was simulated every 6 h. Backward trajectory times for 72 h windows were plotted. Each measured airmass was associated with its corresponding 72 h end observation location.

Satellite Data Validation
Validation of MODIS data was undertaken by comparing the MODIS AOD values to the AERONET AOD values. There were a total of 213 MODIS AOD datapoints with corresponding AERONET AOD values. For QC, there were no available data for the years 2016, 2017, 2019, and 2020, while for KC there were no available data for 2018. This is because either their acquisition time were not within the ± 30 min AERONET measurement time, their locations are not within the 3-by-3 pixel array, or they had fewer than six valid pixels. Table 1 shows the summary of the number of annual data from the MODIS AOD matched with the AERONET data used for the analysis for each study area. High cloud cover and particularly optically thin cirrus pose challenges in AOD retrieval and quality for both ground-based and satellite data measurements [48,49]. Coverage limitations due to the contamination of clouds and site limitations are the primary cause of low data volumes in Southeast Asia, including the Philippines. The scatterplot shown in Figure 2 depicts the correlation of MODIS AOD data and AERONET AOD data from 2010 to 2020. The results show that the R 2 value for the data was relatively high (0.83). Based on the equation of the scatterplot, the coefficient is less than 1 (0.73). This means that the MODIS AOD values were less than the AERONET AOD values. This is in contrast to previous studies showing that the MODIS AOD overestimates AERONET AOD values [50][51][52]. This could be due to the standard deviation of the MODIS AOD values (0.16). The results also show that the satellite data had an arithmetic mean of 0.33 and a low RMSE (0.26). Figure 3a,b show the boxplots for the combined annual and seasonal trends in satellite AOD data. The hollow circles in the graph represent outliers in the data. For the Philippines, the months June, July, and August (JJA); September, October, and November (SON); December, January, and February (DJF); and March, April, and May (MAM) are referred to as summer, fall, winter, and spring, respectively. As shown in Figure 3a, the year 2014 recorded the highest satellite AOD value, as evidenced by an outlier; however, overall, 2015 had the highest mean for the AOD values. There was a general upward trend of AOD values from 2010 to 2015, but it suddenly declined in 2016.   Figure 4 shows the spatial distribution of the mean annual AOD value for the Philippines. The figure suggests a high AOD concentration in the East Asian and Indonesian areas, whereas the Pacific has relatively low AOD values. Eastern China showed continuously high AOD readings within the study years. High mean AOD values in China were due to continuous economic growth and energy consumption [53,54]. Meanwhile, the high AOD values for Indonesia were due to biomass burning. In 2015, Indonesia experienced extensive biomass burning due to El Niño-induced droughts [55,56]. Another major biomass burning event occurred in 2019, resulting in high AOD levels in Indonesia and nearby countries [57,58]. These haze events occur typically from July to October, during the country's dry season [57,59]. These high AOD levels also affect the Philippines. Low AOD values were recorded for 2020 due to the effects of the COVID-19 pandemic and lockdown. This is because strict lockdown protocols in the Philippines restricted the economic production, consumption activities, and transportation of the primary sources of atmospheric aerosols [60][61][62]. Figure 5 shows the boxplot for the annual concentrations for 11 years at the two sites. The mean and standard deviation are presented in Table 2. KC had the highest AOD concentration in 2014 (x = 0.39, σ = 0.28), while QC had the highest concentration in 2010 (x = 0.39). Overall, QC showed the highest concentration within the study period, with a total x of 0.35, compared to KC value of x = 0.32. To test the significance of the difference in concentration values from the two sites, hypothesis testing was performed using an independent t-test. The calculated p-value was 0.38 (p-value ≤ 0.05 means significant). This means that the variation in concentration values between the two sites is not significant, which could be attributed to the sources and quantities of aerosols received by the two study areas. QC is a highly urbanized city, where pollution is primarily caused by anthropogenic activities, but it is also affected by long-range aerosols from East Asia [23]. KC, on the other hand, is a city with a large percentage of its total land area covered by vegetation. Agriculture occupies 55% of the total area, open grasslands occupy 3%, and forest occupies 28% of total land area [63]. While KC also experiences biomass burning from local sources, the bulk of air pollution originates from maritime continent countries in Southeast Asia [19,22,64]. The amount of data could also affect the calculated average AOD values of both sites, since some years are missing data values. Note that "$" corresponds to Quezon City and " " corresponds to Koronadal City.  Figure 5. Boxplots of annual concentrations of evaluated MODIS AOD data from 2010 to 2020 in Quezon City and Koronadal City; median is shown by the middle line of the box, inter-quantile range is shown by the box, and whiskers represent the ± 2.7 + inter-quantile range.

Seasonal Variation
The boxplot of the seasonal variation between the two sites is shown in Figure 6, and the seasonal spatial distribution of AOD values is shown in Figure 7. For the seasonal variation in Figure 6, QC recorded the lowest mean AOD value (0.27) in fall, when the monsoon transition occurs, and the highest (0.40) in spring, as seen in Table 3. A wind back-trajectory analysis (Figure 8) of the northern part of the Philippines during this season showed that air volume comes from the Pacific Ocean by easterly winds, where no known large pollutant emitters are present. During spring, large concentrations of local sources, such as biomass burning and long-range transport pollutants from East Asia, are carried by strong daytime monsoon winds [22]. For KC, the highest AOD values occurred during fall. A wind back-trajectory analysis (Figure 9) during this season showed that the movement of air came from southwest and west of the site. During this time, Indonesia experiences a dry season that encourages biomass burning events that are made severe by high temperatures [64]. This increase in AOD values during fall in KC was a result of forest fires from Indonesia and Indochina [19,65].

Conclusions
This study examined the MODIS AOD trends over Koronadal City and Quezon City from 2010 to 2020. For a decade-long study, the number of data studied is a limitation. This is associated with the limited amount of quality assured and cloud-free data available from the ground-based network and satellite. Nevertheless, the results still show a trend in annual and seasonal AOD values congruent with previous studies.
The results show that the x and R 2 are 0.33 and 0.83, respectively, for the whole study period. The MODIS generally underestimates AERONET AOD values in the Philippines but has a low RMSE value of 0.26. Annual variation in AOD values in the Philippines showed an increase from 2010 to 2015. The mean AOD value in 2015 was recorded as the highest due to extensive biomass burning in Indonesia. Spatial analysis showed that, due to COVID-19 lockdown, the mean AOD value in the Philippines remained relatively low in 2020. Seasonally, KC had a higher AOD concentration during the fall season due to long-range transport pollutants from Indonesia caused by biomass burnings. A wind backtrajectory showed that the high AOD concentration during fall is because aerosols were carried by the southwest monsoon from Indonesia. QC, on the other hand, experienced higher AOD concentrations during spring from increased local anthropogenic activities, as well as long-range transport pollutants from East Asia carried by daytime monsoon winds. These results suggest that subtropical and monsoon climates, as well as the long-range transportation of aerosols, have significant effects on the air quality in the Philippines.
The findings of this study describe and evaluate the atmospheric aerosol trends in the Philippines for the past decade. For an archipelagic country, where some places are remote and inaccessible, using satellite data is highly advantageous because it provides extensive spatial coverage as well as real-time sets of data. This makes it easier to study air quality trends and their effects on health, the environment, and the climate. The authors suggest using other retrieval MODIS AOD algorithms to assess which is the most appropriate for use in different areas within the Philippines. Future studies should also use other satellite sources, such as the Modern-Era Retrospective analysis for Research and Applications (MERRA-2), Visible Infrared Imaging Radiometer Suite (VIIRS), and Himawari-8 satellite, to assess variations in the AOD in different regions in the Philippines, especially in highly urbanized regions. More in-depth studies could be used by the government to enact and initiate interventions in proper urban planning and agricultural practices, as well as health risk assessments to improve air quality in the Philippines.