Aerosol anomalies associated with occurrence of recent strong earthquakes (> M 8.0)

To promote understanding of lithospheric-atmospheric interactions during seismic activity, this study investigated the spatial and temporal variations of aerosol optical depth (AOD) associated with eight recent strong earthquakes (EQs; M w > 8.0) with focal depths of less than 30 km. All EQs occurred after 2009 and therefore have AQUA satellite moderate - resolution imaging spectroradiometer data for the EQ-af-fected areas. This study aims to extract aerosol anomalies before and after large EQs along the coastline. Spatially, AOD anomalies occurred around the EQ epicenters at a monthly scale with longer duration prior to the EQ compared to that after EQ, which may be related to the ground stress variation. Generally, AOD spatial variations increased, then decreased, and finally increased again before returning to normal levels over an 8-day time frame around each event. Temporally, at an annual scale, high anomalous AOD values (> 3.0) appeared before all eight strong EQs with high frequencies of anomalous AOD values occurring both prior to and after EQs. Compared with the control values, AOD values within the affected areas showed an abnormal uplift and downtrend before each EQ event. The AOD values over the control pixels were significantly lower than those over the selected AOD anomaly pixels during the same periods. To eliminate the potential bias from satellites, the result of AERONET daily average AOD data at 500 nm showed that a sudden increase of AOD appeared before and after EQs. AOD anomalies occurred before, during, and after each of the eight strong EQs in offshore regions, partially because of the release of surface latent heat and related atmospheric reactions that increased aerosol concentrations, excluding the effects of the wind field near the earth. The mechanisms behind these AOD anomalies still need to be identified.


INTRODUCTION
Atmospheric aerosol refers to a multiphase system composed of various solid and liquid particles as well as gas carriers suspended in the atmosphere with particle diameters mostly distributed between 0.01 and 100 μm. Natural dust, oceanic particles, and human activities are all sources of aerosol particles. Aerosol optical depth (AOD) is the most basic optical characteristic of atmospheric aerosol and is an important index used to characterize atmospheric turbidity (Wang et al. 2018). With aerosols, AOD is more likely associated with air pollution than with seismic activity. In fact, a correlation between aerosols and seismic activity was first reported back in the 1970s. Tributsch and other researchers identified that aerosol anomalies were related to the occurrence of earthquakes (EQs) (Tributsch 1978). Before an EQ (especially strong EQ), the red sun, red moon, and other light phenomena can be observed mainly because a large number of aerosol particles are emitted from the surface of the EQ region and its surrounding areas into the atmosphere (Ming 2010). Pre-EQ lithospheric seismicity always accompanies the aerosol activity response trace. According to field observations, gases are released, causing atmospheric turbidity prior to strong EQs such as the Haicheng EQ in 1975 and Tangshan EQ in 1976. A significant increase in aerosol concentrations was detected within a radius of 300 km during the Chi-Chi Taiwan EQ in 1999 (Hao et al. 2000). The abnormal behavior of aerosols before and after EQs has also been documented by satellite remote sensing. Ganguly (2016) reported that AOD was 40% higher than normal prior to the occurrence of the Nepal EQ on 25 April 2015. AOD variations obtained using moderate-resolution imaging spectroradiometer (MODIS) measurements acquired by the AQUA satellite were observed to peak 10 days before the strong (M w 7.8) EQ that occurred in Ecuador on 16 April 2016 (Akhoondzadeh et al. 2018). Abnormal AOD values during the Wenchuan EQ in 2008 were observed by Qin et al. (2014a). Results of individual case studies show that AOD has a certain correlation with seismic activity. However, most earlier studies focused on either individual cases or individual regions, without systematic consideration of the correlation between aerosols and EQs on a global scale.
When an EQ occurs, different phenomena are observed underground, on the surface, and from space. Days or weeks before an EQ, the injection of radon and other radioactive materials increases the ionization level in the atmosphere (Virk and Singh 1994;Heinicke et al. 1995;Yasuoka et al. 2006;Omori et al. 2007) and the number of charged aerosols before EQs can increase one or two orders of magnitude (Pulinets et al. 1997). Abnormal ionospheric disturbances have been widely reported to be associated with EQs (Shen et al. 2011(Shen et al. , 2017. One of the possible mechanisms is that air ionization by radon, followed by formation of large ion clusters of aerosol size prior to strong EQs coupled with the vertical turbulent transfer of these charged aerosols has an effect on the electric field, resulting in seismo-ionospheric effects (Liperovsky et al. 2005;Sorokin et al. 2005;Zhuo et al. 2015). Both theory and monitoring results indicate the importance of aerosols leading up to EQs. Currently, there is a lack of monitoring results to develop any seismogenic mechanism to explain this association. Therefore, this study selected strong EQs (M w > 8.0) to explore the correlation between aerosol and the seismic activities and explore its underlying mechanisms.
At present, the monitoring network of aerosols is limited and the distribution of monitoring sites is relatively sparse, although the changes of aerosols detected by satellite remote sensing can make up for the defects of ground monitoring. The MODIS aerosol product was developed more than 20 years ago with six versions of its algorithms being published. The product uses two algorithms: Dark Target (DT) and Deep Blue (DB). Both national and international scholars have carried out verification work for MODIS AOD products. MODIS C4 (fourth edition) AOD products performed well in ocean areas (R 2 = 0.92) but had large deviations over land areas. In response the MODIS fifth edition (C5) AOD products were developed. The establishment of the global Aerosol Robotic Network (AERONET) facilitated the validation of MODIS AOD products. Levy et al. (2010) verified that MODIS C5 AOD products had a high level of accuracy using AERONET data. In 2014, the MODIS team upgraded the DT and DB algorithms to form the sixth version (C6) AOD product, with a slightly smaller error and better coverage rate than that of C5 as verified by Levy et al. (2013). In this study, MODIS AOD data were selected with relatively reliable quality assurance for longterm monitoring of AOD fluctuations in seismogenic zones.
Based on the above research, the main purpose of this study is to explore the anomaly variation of atmospheric aerosols before and after global EQs (M w > 8.0), and to characterize aerosol anomalies at various spatial and temporal scales in seismic regions by satellite and in-situ data. The common characteristics of aerosol anomaly variations in seismogenic zones are summarized below.

Data
The twin MODIS sensors have been flying on the TERRA satellite since 2000 and the AQUA satellite since 2002, creating an incredible dataset of global Earth observations. This study selected the AQUA MODIS data because they avoid other interferences caused by various external factors such as orbit. AQUA was launched later than TERRA and its data are more reliable and accurate for the chosen period of our EQ study. Monthly and 8-day MO-DIS data acquired by the National Aeronautics and Space Administration's Earth Observing System AQUA satellite were used to retrieve AOD (Level 3 aerosol product) data from 2004 to 2017. This study used the AOD product from AQUA MODIS C6, including MYD08_E3 and MYD08_ M3 with "Combined Dark Target and Deep Blue AOD at 0.55 micron for land and ocean: Mean of Daily Mean". The data set, presented on a Gaussian grid of the earth's surface having 1° × 1° spatial resolution, was downloaded from: https://ladsweb.modaps.eosdis.nasa.gov. The AQUA satellite, launched in May 2002, is a polar-orbiting satellite that operates at an altitude of approximately 700 km, providing aerosol data every 1 -2 days. It crosses the equator at about 1:30 pm (ascending orbit) local sun time with a scanning swath of 2330 km (cross track). The expected error of the MODIS AOD measurements is ±0.05 (±0.15 × AOD over land) (Levy et al. 2013). The AOD trends calculated using AERONET instantaneous data at MODIS overpass times were consistent with those of AERONET daily data, which suggests that the AOD trends derived from the satellite measurements of 1 -2 overpasses are sufficient representatives of daily measurements (Remer et al. 2005;Levy et al. 2010). The EQ data were obtained from the US Geological Survey website (https://earthquake.usgs.gov).
In this study, the in-situ data are from the global aerosol automatic observation network AERONET (Aerosol Robotic Network) initiated and organized by NASA, which can realize automatic measurement, collection, storage, and transmission of data. The instrument used is the CE318 automatic scanning solar photometer produced by the French company CIMEL. There are 8 observation channels in the visible and near-infrared bands. At the end of 2010, more than 600 AERONET observation sites were distributed globally on land and sea. This study selected AERONET Level 2.0 products (Cloud Removal and Quality Inspection). The AERONET daily average data of AOD were obtained from the website (http://aeronet.gsfc.nasa.gov/). ERA-Interim is a global atmospheric reanalysis from 1979, which is continuously updated in real time. The data assimilation system that produces ERA-Interim uses the December 2006 version of the ECMWF Integrated Forecast Model (IFS Cy31r2), which originally covered dates from 1 January 1989, but later added with a past decade from 1 January 1979. The data assimilation is based on a 12-h fourdimensional variational analysis (4D-Var) with adaptive estimation of biases in satellite radiance data (VarBC). Values at four times a day, i.e., at 00:00, 06:00, 12:00, and 18:00 UTC are provided. To match with the MODIS satellite data, only local morning data of the 10 m U wind component and the 10 m V wind component were used, i.e., at 00:00 UTC. Spatial resolution is 0.75° in latitude and 0.75° in longitude over EQ preparation area.

Methodology
This study considered eight strong EQs (M w ≥ 8.0) with a focal depth of < 30 km, which occurred in various parts of the world during 2004 -2017 (Table 1). The appearance of EQ precursors usually occurs on large-scale areas around the impending faults. Dobrovolsky et al. (1979) proposed an empirical relationship (ρ = 10 0.43M ) between the supposed circular size of mechanical/thermal/electromagnetic precursors and magnitude of the eventual EQ where ρ is the radius in km of the interested area, called 'strain radius' and M is the magnitude of the impending EQ. Taking into account the empirical relationship and considering the magnitude 8 of the strong EQ, an area with a radius of about 2500 km was selected for data extraction.
The background values of AOD were based on mean AOD values, calculated using multiyear MODIS data by extracting the total column AOD values over the epicentral areas of these EQs. The background values of AOD were derived using Eq. (1): where G bac (x, y, t) is the background value, i.e., the multiyear mean AOD value for an area of latitudinal length x and longitudinal length y in period t (12 months or 8 days) during N years prior to a given EQ, excluding the EQ year. Both x and y were about 1°. The anomalies were derived using Eq. (2), which is used to define anomalies of thermal and gas fields (Tramutoli 1998(Tramutoli , 2007Tramutoli et al. 2001;Ouzounov et al. 2007Ouzounov et al. , 2018Cui et al. 2017). Index of ALICE AOD (x, y, t) is applied to detect anomalies, which represents the anomalous amplitude within a given region during a given time. Here ( , ) x y v is the standard deviation, as described by Eq. (3). Taking the 11 March 2011 earthquake in Japan as an example, time series of the AOD background field were selected from 2004 to 2010, given N = 7. Thus, G(x, y, t) represents the AOD value in 2011 for the latitudinal length x and longitudinal length y.

RESULTS
This study investigated the correlations between AODs and the occurrence of strong EQs in space and time. Spatial correlations were assessed at monthly and at 8-day time scales, while temporal ones were assessed over longer time scales. Figure 1 showed the spatial distributions of monthly AOD anomalies before, during, and after the occurrence of each EQ (cases 1 -8).

Characteristics of Spatial Patterns of Aerosol Anomalies-Monthly Scale
Monthly scale data were used to study the correlation between EQ activities and AOD values. The parameters of the strong EQs (M w > 8.0) and the general features of AOD anomalies are described in Table 1. Abnormal AOD phenomena were typically observed before the month in which the strong EQs occurred. From Table 1, AOD anomalies occurred before and after all EQs but with different durations. In general, the duration of the AOD anomaly in each EQ case was less than 6 months before the EQ, with longer AOD anomaly duration before EQs than after EQs. However, some cases did not have obvious AOD anomalies in the month in which the EQ occurred, which may be related to the locking status of the fault during EQ development (Cui et al. 2017). Moreover, the AOD monthly average results undercut AOD daily anomalies when EQs occurred. From a vertical perspective, abnormal AOD response areas are observed two months before the EQs in all cases, except for cases 6 and 8. Moreover, AOD anomaly enhancement after each event may be related to a series of aftershocks in the region. Figure 2 shows the abnormal AOD changes before, during, and after EQ with a 12-month interval around the occurrence of case 1 for an EQ with a magnitude of 9.1. On a spatial scale, there were obvious abnormal phenomena prior to and during this strong EQ. The abnormal AOD values during the month in which the EQ occurred were located in a narrow NE-SW zone around the epicenter with already high values in the northwest zone, two months before the event. One month after the event, the AOD anomaly was weaker but strengthened again in the second month after the event. At this point, the anomaly positions were more scattered, which is related to aftershock activity (Cui et al. 2013). The anomaly pattern may reflect aspects of regional geological structure and tectonic movement. However, it is noteworthy that the AOD values first increased and then decreased before case 1, becoming weaker after the EQ, which was similarly observed in other EQ cases. In conclusion, an obvious AOD anomaly zone occurred before strong EQs at a monthly scale, which is most likely related to the intensity of the EQ (M w > 8.0).

Characteristics of Spatial Patterns of Aerosol Anomalies-8-Day Scale
In addition to the monthly scale, this study explored the spatial correlation of aerosols and EQs over an 8-day scale within a 3-month interval, encompassing 2 months prior and  1. Spatial distributions of monthly AOD anomalies before, during, and after the occurrence of the studied earthquakes (1 -8). The epicenters, the AERONET sites, and the control points were indicated by the red stars, yellow circles, and the red '*', respectively.
1 month after each EQ event. Case 1 is again described as a typical EQ example, as shown in Fig. 3. An obvious AOD anomaly zone appeared near the epicenter during 2 to 9 February. Moreover, abnormal AOD values near the epicenter are still observed after the EQ (14 to 30 March), which are possibly related to aftershocks. Figures 4 and 5 showed the spatial distributions of 8-day AOD background field and AOD standard-deviation field before, during, and after the occurrence of 11 March 2011 Japan earthquake. The AOD background field in the study area gradually increased from January to April. However, standard-deviation calculation result showed that the standard deviation of AOD was small before March but was larger in the northern region of the epicenter in April, which might be related to climate changes. More importantly, the AOD anomaly regions in all cases were consistent before and after the EQ at this scale. Clearly, AOD spatial variations at this scale also increased, then decreased, and again increased before returning to normal levels. All eight cases had abnormal regions before the EQ event.

Characteristics of Temporal Patterns of Aerosol Anomalies-8 Scale
Spatially, the AOD anomaly in all cases is observed to be around the EQ epicenter. Herein, the correlation be-tween AOD and EQs over various time frames is explored. The temporal variations of 8-day AOD values are shown in Fig. 6. Every EQ event was subdivided into two parts: (1) The upper panel shows the original AOD values (black lines) and anomaly values (red vertical bars) and (2) The lower panel shows all EQs with M w > 5.0 that occurred within the study area. At an annual scale, high anomalous AOD values (> 3.0) appeared before all eight strong EQs with high frequencies of anomalous AOD values occurring prior to and after the EQs. Compared with the control values, AOD values within the affected areas showed an abnormal uplift and downtrend before each EQ event. In case 1, the initial AOD values obviously increased significantly before and after the strong EQ and the anomalous AOD values were large. Some abnormal AOD values also appeared at a later stage, which is possibly related to the strong aftershocks associated with case 1. Similar results are observed for all other EQs.
In this study, pixels without anomalous values located far from the epicenters and fault zones served as the control sites. The AOD values at the control sites were compared to those in the epicentral areas of the eight EQs (Fig. 7). The AOD values over the control pixels were lower than those over the selected AOD anomaly pixels for the same period. Using the method of anomaly, control values lower than 3 are obtained. Control AOD values also did not fluctuate strongly before or after EQs.

Characteristics of Temporal Patterns of Aerosol Anomalies-Daily Scale
To eliminate the potential bias from satellites, in-situ data of AOD, i.e., AERONET daily average data of AOD at 500 nm were chosen as the excellent data source for any of the 6 events (no AERONET sites around the remaining two EQs). The related AERONET information was listed in Table 1 and the spatial distribution of AERONET sites and the epicenters of EQs were shown in Fig. 1. Since there was no AERONET site around the epicenters of case 7 and case 8, temporal variations of AERONET AOD daily average data at 500 nm around the six EQ epicenters with the time series of three months prior to EQs and two months after EQs were shown in Fig. 8. Results of the six EQ cases showed obvious signs of AOD's sudden increase before or after the EQs. Taking the event occurred on 11 March 2011 as an example, the AOD continued to increase significantly a month before the EQ during 4 to 8 February 2011. The AOD fluctuations were relatively stable from 1 December 2010 to 3 February 2011 and were observed to be after the EQ, which may be mainly related to the aftershock activities. It can be seen that the aftershocks with the magnitude of greater than 5 were very obvious, which is consistent with the abnormal fluctuation of AOD after the EQ.

DISCUSSION
Both abnormal spatial and temporal variations of AOD are attributable to lithospheric crustal stress, geological processes, and geochemical processes before and after the occurrence of a strong EQ. From the perspective of crustal stress, seismic regions usually go through a stable quasilinear strain accumulation stage and then enter a non-linear strain stage. When strain accumulation has reached a certain threshold, the strain shows a slow-release effect, usually accompanied with a small increase in the pre-slip and EQ activity. Finally, a large rupture occurs, releasing the strain and resulting in a strong EQ. Before the strong EQ, stress extrusion and crushing of rocks could produce a thermal anomaly. Latent heat fluxes have been observed before EQs in unsaturated water vapor (Qin et al. 2014b;Wu et al. 2016). Traditionally, aerosols are considered to be the nuclei for water condensation in the atmosphere (Sorokin et al. 2001).
From the perspective of the geochemical processes, gaseous emissions from active tectonic areas are wellknown phenomena (Heinicke et al. 2006). Moreover, a large increase in the radon concentration and a marked thermal anomaly have been recorded before and after the occurrence of major EQs. For instance, the radon concentration reached a peak level of more than 10 times than that recorded before the Kobe EQ in 1995 (Igarashi et al. 1995). Radon Fig. 3. Spatial distributions of 8-day AOD anomalies before, during, and after the occurrence of 11 March 2011 Japan earthquake. Obvious AOD anomalies can be seen near the epicenters (black stars stand for the epicenter and the red '*' for the control site).     concentrations gradually increase over several months before a strong EQ. Radon decay could stimulate the nearsurface electric field, producing ionization of other gases (Yasuoka et al. 2006) and heat release. Ions created by ionization form centers of condensation that help release latent heat (Sorokin et al. 2001). Ion-induced nucleation process follows, increasing the number of cluster ions, which continue to grow in size. In this case, ion-induced nucleation can be considered as a catalytic exothermic reaction, in which ions play the role of the catalyzer. Moreover, these high ion concentrations lead to explosive nucleation processes, forming ion clusters of several microns in size (Pulinets and Ouzounov 2011). Experimental results of Aerosol Optical Thickness measurements from the AERONET network demonstrated that spikes in aerosol particle size distributions at around 1000 nm are observed before strong EQs . Therefore, aerosol anomalies typically occur before EQs. EQs occur because a deadlocked fault segment cannot undergo deformation due to the elastic deformation limit. In case 1 (Fig. 3), it could be seen that the AOD anomaly also obviously appeared in the southeast area of the epicenter from 25 January to 1 February and an obvious AOD anomaly zone then appeared near the epicenter during 2 to 9 February. It might be suspected that the AOD anomalies near the epicenter during 2 to 9 February are affected by the AOD high values of the southwest direction from 25 January to 1 February under the meteorological condition, not by seismic activity. To exclude the influence of meteorological factors on AOD migration, the trend of near-Earth 10-m wind field was extracted (see Fig. 9). From 25 January to 1 February, the wind direction was clearly observed to be basically northwest. During this period, the abnormal AOD of southwest direction affected by the northwest wind had a small effect on the AOD near the epicenter during 2 to 9 February. Therefore, it was further determined that the AOD anomalies near the epicenter during 2 to 9 February in Fig. 3 were most likely related to seismic activity. Meanwhile, it was reported that the stress orientation prior to EQ was NE-SW within the seismically active zone, consistent with the AOD anomaly phenomenon observed beforehand at both monthly and 8-day timescales. This suggests that there was a local stress anomaly related to the seismic activity (Lin et al. 2011). Typically, the anomaly durations prior to the EQs were longer than those after the EQs (Table 1), which may be related to EQ activity patterns. The release of heart energy and gaseous emissions during early seismic activity is slow. However, once the fracture zone forms during the EQ, a direct release of heat energy and gaseous emissions occurs following the EQ, triggering air ionization (ions start to form but have not grown to a low mobility level). As ion 'ageing' occurs, clusters form, which grow to aerosol size (from 1 to a few microns) . These clusters have extremely low mobility, which develops the NE-SW AOD anomaly along the active fault. Therefore, the aerosol anomaly preceding the EQs exists longer because it is related to the variation in ground stress during the early seismogenic process. This study considers all strong EQs that occurred on a coastline. Compared with aerosol retrievals over complex land surfaces, retrievals over the ocean are more accurate with less interference from non-seismic factors, which is conducive to capturing seismic aerosol anomalies. Despite Fig. 9. 10-m wind field near the earth. The black star represents the epicenter location of the case 1 and the yellow area is the land boundary. Latitude (North) and longitude (East) are in degrees. this, little research has explored the mechanism underlying seismic aerosol anomalies along the coast. In contrast to terrestrial aerosol anomalies, aerosol anomalies over the ocean, associated with offshore EQs, are likely related to the release of latent heat. The release of latent heat is caused by the condensation of water vapor ). The condensation of water vapor related to the release of latent heat will interact with salty ocean air to produce new aerosols, affecting aerosol concentrations along the coastline. In fact, a fog is observed very often over EQ-affected areas . The formation of such fogs cannot be separated from that of aerosols. The contrast in seismic anomaly index values for anomaly and control areas reflects the fact that aerosol concentrations are correlated with seismic activity.
Therefore, both seismic and geochemical processes likely lead to greater aerosol formation, which produces the AOD anomalies observed before and after strong EQ events, as well as the monthly and 8-day patterns related to strong EQs. Aside from seismic activities, aerosol anomalies could also be affected by climate change, human activities, AOD retrieval methods, and vegetation cover. In this study, the effects of background factors were eliminated using the method outlined in Eqs. (1) -(3). Moreover, large EQs along the coastline are less affected by human activities and vegetation cover, making the accuracy of data inversion higher than for other terrestrial settings. Therefore, the AOD anomalies described herein could be considered as solely EQ-related features.

CONCLUSION
AOD anomalies were found to occur around epicentral areas at a monthly scale in all eight strong EQs investigated in this study. The AOD anomaly enhancement after the EQ events was likely related to the series of aftershocks following strong EQs. The AOD spatial variations initially increased, then decreased, and again increased before returning to normal levels over an 8-day timeframe. On an annual scale, anomalous AOD values (> 3.0) appeared before all eight strong EQs with high frequencies of anomalous AOD values occurring prior to and after EQs. Compared with the control values, AOD values within affected areas showed an abnormal uplift and downtrend before each EQ event. The anomaly duration prior to the EQ was longer than that after the EQ, which may be related to the variation in ground stress. The contrast in seismic anomaly index values for anomaly and control areas reflects the correlation of aerosol concentrations with seismic activity. Through elimination of the effects of background factors and excluding the wind field of meteorological conditions, this study clarified that AOD anomalies observed both before and after strong EQs could be related solely to seismic features, based on the analysis of the lithosphere-atmosphere coupling and the variation in ground stress. Aerosol anomalies over ocean areas, associated with offshore EQs, are likely related to the release of latent heat. The condensation of water vapor caused by the release of latent heat interacts with salty ocean air to produce new aerosols, further affecting aerosol concentrations along the coastline. Continuous study of AOD values using satellite-based observations could be employed to investigate the lithospheric-atmospheric processes associated with strong EQs.