Identification of atmospheric transport and dispersion of Asian dust storms

Backward trajectories of individual Asian dust storm (ADS) events were calculated using the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) at four representative stations in Korea. A total of 743 ADS events and associated 2229 (endings of altitudes at 1000, 1500, and 2000 m per ADS event) backward trajectories from four stations were traced from January 2003 to August 2015. Regardless of the locations of the observed stations and the threshold time divide, a recent increase in the ADS occurrence rate was statistically significant with a 99.9 % confidence limit. Winter and spring were high-occurrence seasons for the ADS, while it rarely occurred in summer. Angular distributions of dust transport indicated a dominance of northwesterly wind, as more than two-thirds of ADS events are azimuthally confined from 290 to 340. In addition, there is a tendency for stronger PM10 dust air concentration to be from the northwest. We found a strong inverse correlation between the number of days with ADS events and cumulative PM10 dust air concentration, indicating that the total amount of cumulative PM10 discharge was rather constant over time. If so, relatively shorter transport distances and a more continental dust passage over the Shandong peninsular would yield less PM10 in a shorter transport path but with a stronger concentration.


Introduction
Asian dust storms (ADSs) originate mostly in the Gobi and Taklamakan deserts, where high-speed surface winds and intense dust storms took up dry, dense clouds of fine-grain surface material.Once elevated, ADSs spread in the strato-sphere, often reaching to the upper troposphere to form the dust devil (Yumimoto et al., 2009).These dust clouds are then carried eastwards by prevailing westerly winds and pass over China, Korea, Taiwan, and Japan to the North Pacific (Hsu et al., 2012).These dust clouds were eventually transported for more than one full circuit around the globe in about 2 weeks (Uno et al., 2009).Such long-range transport is possible as the dust clouds are transported in atmospheric dust layers (Tratt et al., 2001;Schepanski et al., 2009;Uno et al., 2009).It has been estimated that ADSs contribute about 4 % of the global dust emission (Mahowald et al., 2010).
ADSs contained surficial minerals of natural origin (e.g., weathered soils) as well as pollutants of anthropogenic origin such as black carbon, heavy metals, and sulfates (Guo et al., 2004;Hsu et al., 2004;Ramana et al., 2010).These pollutants have an impact on climate variations as black carbon absorbs solar radiation (Jacobson, 2001(Jacobson, , 2012)), while sulfates scatter solar radiation energy (Ramana et al., 2010).Deposition of anthropogenic pollutants also affects the human epidermal keratinocytes (Choi et al., 2011), dispersion of bacterial cells (Yamaguchi et al., 2012), and marine biogeochemical environment (Lin et al., 2007;Hsu et al., 2009).In addition, longrange transport of ADSs influences cardiovascular disease in Taiwan (Chen and Yang, 2005), cardiopulmonary emergence rate in Taiwan (Chan et al., 2008), asthma in Japan (Watanabe et al., 2011), daily mortality in Korea (Lee et al., 2013), and food toxicity in Japan (Kobayashi et al., 2016).Hence it is natural to raise public concern of the possible adverse effects of ADSs in eastern Asia.
ADSs influences air quality and the local ecological environment including vegetation and soil.It accelerated the process of desertification by increasing the rate of evaporation.
Most of all, it has been observed that the frequency of ADS occurrence in eastern Asia increased as a result of desertification, overgrazing and over farming, irrigation, and lack of precipitation in central Asia.Of course, it is also true that ADSs are dependent on synoptic climatic conditions such as cyclone activity, air temperature, and precipitation in the source area (Natsagdorj et al., 2003;Qian et al., 2002;Hsu et al., 2013).
Over the past few decades, variation of ADSs and their climate control has been explored (Natsagdorj et al., 2003;Qian et al., 2004;Kim et al., 2008;Lee et al., 2010;Zhao et al., 2010).For instance, Natsagdorj et al. (2003) compiled dust storms in Mongolia from 1937 to 1999, on the basis of observational data from 49 stations and found that the annual mean number of days with dust storms were 20-37 days in the Gobi desert and nearby arid area.Wang et al. (2005) defined three different types of dust storms and their characteristics in China, using data from 701 meteorological observation stations from 1954 to 2000.It should be highlighted that the number of ADS events increased by more than twice from 1960 to 2000 (Zhang et al., 2003;Wang et al., 2007Wang et al., , 2008)).
Previous studies in Korea focused on the characteristics of dust particles and the magnetic concentration of dust particles in ADSs (Chun et al., 2001(Chun et al., , 2008;;Chung et al., 2003;Lee et al., 2006Lee et al., , 2013;;Kim et al., 2008).In the present study, we trace the air parcel trajectories of ADSs using the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Draxler and Hess, 1998).Temporal and spatial variation of ADSs allows for a more comprehensive view of the dust generation and transport in eastern Asia.In particular, tracing the ADS trajectories in Korea is pivotal because Korea is the first country encountered out of the source region by ADSs on its westerly dominating dust transport.It is the purpose of this study to gain insight into how atmospheric transport, dispersion, and deposition of dust particulates occur in eastern Asia.

The HYSPLIT model
The HYSPLIT has evolved from the earliest model in 1982 (Draxler and Taylor, 1982) from modeling long-range air parcel trajectories into simulations of pollutant transportation, dispersion, and deposition over global scales.The HYSPLIT model relies on the determination of air concentration as a cumulative summation of dust flow per unit grid cell.Each dust flow is considered an independent particle flow puffed by advection and is represented by its trajectory.Backward trajectories are constructed on the basis of Stochastic Time-Inverted Lagrangian Transport (STILT) model.The STILT is a widely used model for tracing atmospheric mixing between the source and the receptor point in terms of 2-dimensional upstream surficial fluxes.
As an open-source model, HYSPLIT is available on the web through the ARL READY system (http://ready.arl.noaa.gov/HYSPLIT.php),operated by the National Oceanic and Atmospheric Administration (NOAA) Air Resource Laboratory (ARL).In the present study, we used HYSPLIT from the most recent model in September 2015 (Stein et al., 2015).
The HYSPLIT model requires the meteorological data and vertical movement of atmospheric circulation as input, and it displays the analysis of the simulation outputs (Stein et al., 2015).The HYSPLIT model describes transport and dispersion dynamics of aerosol, incorporating boundary stability determined by turbulent velocity, a wind-blown dust emission algorithm, convectional plume rise produced by buoyancy of heat, wind velocity, atmospheric friction velocity, and in-cloud wet scavenging.One great advantage of using the HYSPLIT model is that both forward and backward trajectories are available with local or global airflow patterns to interpret the transport of pollutants.The HYSPLIT model is continuously evolving to cope with turbulent mixing process and to incorporate higher temporal frequency data available from the meteorological data (Stein et al., 2015).
In general, PM 10 dust air concentration less than 25 µg m −3 is considered to be recommended by the World Health Organization.According to KMA, advisory warning is issued when the hourly averaged PM 10 dust air concentration is expected to exceed 150 µg m −3 per hour.When the hourly averaged PM 10 dust air concentration is expected to exceed 400 µg m −3 per hour, a more significant danger warning is issued.In the present study, an individual ADS event was counted as a day with an advisory warning.

Results
Backward trajectories of individual ADS events were calculated using the HYSPLIT at four representative stations at BR, DJ, KS, and UR (Fig. 1a).The data archive includes a collection of HYSPLIT model since 2003 based on NCEP/NCAR reanalysis 1 at the time of ADS events.Trajectories of dust transport at altitudes of 1000, 1500, and 2000 m were traced for 72 h (Fig. 1b).A total of 743 ADS events and associated 2229 (three different endings of altitude per ADS event) backward trajectories from four stations were traced from January 2003 to August 2015 (Table 1).Then, ADSs were classified into six transport paths (Fig. 1b) on the basis of wind azimuth during the first 24 hours of each ADS event as N-NE (northerly 0 • to northeasterly 45), N-NW (northerly 360 • to northwesterly 315 • ), W-NW (westerly 270 • to northwesterly 315 • ), W-SW (westerly 270 • to southwesterly 235 • ), S-SW (southerly 180 • to southwesterly 235 • ), and S-SE (southerly 180 • to southeasterly 135 • ) (Table 1).
Annual cycles of the number of days each year with ADS events (n ADS ) are displayed (Fig. 2).A modern data set with PM 10 observations from 2003 to the present was combined with an old data set without PM 10 observations (Fig. 2).For the past 13 years, the cumulative number of days with ADS events was 111 at BR, 220 at DJ, 257 at KS, and 155 at UR (Fig. 3a, Table 2).For each station, annual and monthly variations of (n ADS ) were similar (Fig. 3b, c).However, it is interesting that (n ADS ) of BR was systematically lower than that  for other stations (Fig. 3b).For the past 13 years, the lowest frequency of ADSs occurred in 2012 (Fig. 3b).Winter and spring were high-occurrence seasons for the ADS, while the ADS rarely occurred in summer (Fig. 3c).
Temporal variations of dust air concentration observed in each station were represented with time (Figs. 4 and 5).For each station, the distribution of individual PM 10 dust air concentration was plotted as a function of time in lower panel (Figs. 4 and 5).In the upper panels (Figs. 4 and 5), results were rearranged in box plots in which the central box represents the interquartile of annual mean PM 10 dust air concentration and whisker lines extend beyond the maximum and the minimum.For BR, a total of 111 ADS events from 2003 to 2015 showed mean PM 10 dust air concentration of 424.7 ± 341.1 µg m −3 (Fig. 4a, Table 2).On the other hand, values of mean PM 10 dust air concentration in other stations were smaller than those at BR (Table 2).For instance, DR, KS, and UR recorded mean PM 10 dust air concentrations of 189.5±94.2µg m −3 (Fig. 4b), 190.7±62.3µg m −3 (Fig. 4c),   and 188.3 ± 55.7 µg m −3 (Fig. 4d), respectively.Days on which danger warnings were issued with PM 10 dust air concentrations exceeding 400 µg m −3 were 37 for BR (33.3 %), 5 for DJ (2.3 %), 7 for KS (2.7 %), and 3 for UR (1.9 %).Values of PM 10 dust air concentration tended to be larger over spring (Fig. 5).On a log-log plot, PM 10 dust air concentration observed from BR plot in a line with a slope of −1.5, meaning that an order increase in PM 10 dust air concentration correlates with a 50-fold decrease in ADS occurrence (Fig. 6a).For UR, the slope of maximum asymptotic power fitting yielded −4.0, an order increase in PM 10 dust air concentration correlating with a four-order decrease in ADS occurrence (Fig. 6a).Results for DJ and KS were confined within the trends of BR and UR (Fig. 6a).They were definitely not linear, but are tailed towards stronger values of PM 10 dust air concentration for lower occurrence (Fig. 6a).We have constructed cumulative probability distribution functions for each station (Fig. 6b).Values of median PM 10 dust air concentration were 318.72 µg m −3 for BR, 162.13 µg m −3 for DJ, 164.49µg m −3 for KS, and 171.62 µg m −3 for UR (Fig. 6b).Tracing the dust transport using the HYSPLIT was represented with an angle histogram plot, which is a polar plot showing the distribution of prevailing wind in angle bins of 5 degrees (Fig. 7).For the past 13 years, the spatial distribution of ADS events is prominently northwesterly (Fig. 7).

Discussion
Because of its short operation history, only the modern data are available in BR (Fig. 2a).More extended temporal coverages over 50 years of ADS events were available in DJ (Fig. 2b), KS (Fig. 2c), and UR (Fig. 2d).It is notable that the annual occurrence rate of ADS events has increased recently (Fig. 2).For instance, the annual mean occurrence rate of ADSs in Daejeon was 4.1 ± 4.0 from 1960 to 2000 and 17.2 ± 6.5 from 2001 to 2015, which were significantly different (p = 7.611 ± 10 −7 ) to each other according to the criteria of Welch's t test with a significance of 99.9 % (Table 3).Regardless of the locations of the observed stations and the threshold time divide (either 1998 or 2001), the recent increase in ADS events was statistically significant with a 99.9 % confidence limit.It is also apparent that there may be also a decline in ADS occurrence with respect to 2007/2008.Results for Welch's t test are of the order of 10 −2 , slightly over the statistical threshold for the acceptance of the null hypothesis.Despite its failure of the null hypothesis, we cannot completely rule out the possibility as the temporal separation in 2007/2008 leaves only 8 data points for recent intervals (Table 3).
Year-to-year variations (Fig. 3b) and month-to-month variations (Fig. 3c) of ADS events were compared.It is apparent that ADSs are heavily concentrated in the dry seasons, from February to June (Fig. 3c).The highest dust air concentration that occurred was 2371 µg m −3 on 8 April 2006 at BR, 862 µg m −3 on 12 September 2010 at DJ, 1342 µg m −3 on 2 April 2007 at KS, and 440 µg m −3 on 2 April 2007 at UR (Figs. 4 and 5).Occurrence of more frequent ADSs and  their seasonality was also documented in neighboring countries including Japan and Taiwan (Yang, 2002;Watanabe et al., 2011;Chien et al., 2012;Kimura, 2012).Then, it can be understood that seasonal variation of ADSs is a regional, natural phenomenon in the eastern part of Asia.Several factors are known to be closely related to the occurrence and strength of dust storms including extinction of pasture, abandoning cropland without vegetation cover, overexploitation of forests and shrubs, enhancement of mining activities, unregulated wild roads desertification, increase in human or factory transport, and decline of annual precipitation (Xuan et al., 2004;Aoki et al., 2006;Bian et al., 2011).Increase in desertification results in high dust emission from the places with less vegetation cover with dry and loose sandy soils in the Gobi and Taklamakan deserts (Laurent et al., 2005(Laurent et al., , 2006;;Zhang et al., 2008).The recent increase in ADS events is contemporaneous with the increased desertification in the Gobi and Taklamakan deserts.For instance, the number of the dust days in the Gobi desert have been almost tripled in recent years when compared to those in the 1960s (Natsagdorj et al., 2003).On another occasion, it has been reported that sandy desertification in northern China increased rapidly with mean annual areal expansion of 2460 km 2 (Zhang et al., 2008).In addition, a decrease in the amount of soil moisture and increase in mean wind speed provides more frequent generations of dust storms (Natsag-dorj et al., 2003).It should be highlighted that the highest frequency of ADSs arose from the Gobi deserts, mostly in spring due to the development of lower air temperature in winter and high frequencies of cyclone activities (Qian et al., 2002).In addition to natural pedogenic enhancement, we cannot ignore the contribution of anthropogenic particulate matters supplied by fossil fuel combustion, coal burning, and industrial plants.Although anthropogenic particulate matters represent only 5-30 % of ADS volumetrically, they are harmful as they have a strong tendency to react with heavy metals.Considering ongoing demand for fossil fuel combustion, it is reasonable to suggest that pollutants of anthropogenic origin are also responsible for the increase in ADSs.
To evaluate the differences between the data sets observed from four different stations, a Kolmogorov-Smirnov test of nonparametric and distribution-free statistics was applied.When the value of P is insignificant (e.g., P <0.01), we can reject the null hypothesis of no difference between the two data sets.According to the Kolmogorov-Smirnov test, data sets KS and UR are similar (P = 0.7040) and those between DJ and KS are somewhat similar (P = 0.0280) with a 98 % confidence limit.Contrary to the initial impression based on the correlation between ADS occurrence and PM 10 dust air concentration (Fig. 6), however, all other pairs of data sets are statistically different (Table 4).Angular distribution of dust transport indicates a dominance of northwesterly wind  The maximum difference between the cumulative distribution (D) and corresponding probability (P ) according to the Kolmogorov-Smirnov test.
(Fig. 7).In fact, the prevailing ADS is from a narrow sector between 290 and 340 • (Fig. 7).There is a tendency for stronger dust densities to be from the northwest (Fig. 7).
The HYSPLIT backward trajectories at different altitudes of 1000, 1500, and 2000 m were counted as individual paths in the present study.They do not show meaningful differences statistically, implying that atmospheric turbulent mixing was minimal.Such directional consistency for different altitudes of the HYSPLIT model might result from the relatively low and flat geographic conditions.For instance, both eastern China and western Korea are low in elevation, bridged by the shallow Yellow Sea which extends 900 km in a north-south direction and 700 km in an east-west direction.
As the ADS is a randomly occurring natural hazard, the number of individual ADS events (n ADS ) should follow the exponential probability density distribution with the weaker dust air concentration occurring more frequently (Fig. 6a).As anticipated, ADS occurrence and PM 10 dust air concentration displayed an inverse power relation in all four stations (Fig. 6a).But why did each station record a different power relation?
Thermodynamic equilibrium of wind-blown dust requires competing balance among uprising buoyancy, gravitational settlement, and flow resistant drag force.Strong inverse correlation between the number of days with ADS events (n ADS ) and cumulative PM 10 dust air concentration implies that the total amount of cumulative PM 10 discharge is rather constant over time.As a result of the nearly constant PM 10 discharge, relatively shorter transport distance along with more continental dust passage through Shandong peninsular would yield less but stronger PM 10 dust air concentration in BR (Fig. 8a).In other words, longer settlement time intervals were required for DJ, KS, and UR (Fig. 8a) simply because they are located farther than BR from the dust sources (Fig. 1).
Mean PM 10 dust air concentration was estimated by dividing the total amount of cumulative PM 10 flux to the total number of ADS events (n ADS ).In fact, mean PM 10 dust air concentration (PM 10mean ) is equivalent to the arithmetic mean of PM 10 dust air concentration for each station.However, it is true that individual dust air concentration distribution is far from being Gaussian or lognormal (Fig. 6b).Instead, median PM 10 dust air concentration (PM 10median ) estimated from the cumulative dust air concentration distribution is useful to reflect individual PM 10 dust air concentration distribution.Regardless of the distribution of natural hazards, the median is inherently more stable than the mean with respect to the uncertainty and will change less over time (e.g., Atkinson and Goda, 2011).Nonetheless, as far as the present study is concerned, PM 10mean and PM 10median can both reflect the mean PM 10 dust air concentration as they are positively correlated (Fig. 8b).

Conclusions
The present study dealt with the Asian dust storm (ADS) outbreaks affecting Korea from January 2003 to August 2015.A total of 743 ADS air parcel backward trajectories reaching Korea were identified by means of a Lagrangian integrated trajectory (HYSPLIT) at three different ending altitudes of 1000, 1500, and 2000 m.In all four stations where the ADS was monitored, we found that ADS occurrence rate recently increased.This increase in ADS occurrence was statistically significant with a 99.9 % confidence level regardless of the threshold time divide of 1997/1998 or 2000/2001.Monthly variation of ADS occurrence was definitely nonuniform, as ADS was mostly concentrated in the colder seasons of winter and spring.Instead, ADS events rarely occurred from June to September.Majorities of ADS events are azimuthally confined in narrow intervals of 290-340 • on angle histograms, indicating that a northwesterly distribution of dust transport was prominent.Such angular dependence of ADS occurrence agrees well with the higher PM 10 dust air concentration from the northwest.We propose that the total amount of cumulative PM 10 discharge was rather constant over time in Korea, as there is an inverse correlation between ADS occurrence and PM 10 dust air concentration.Such constant PM 10 flux allows weaker PM 10 concentration for longer transport, and vice versa.
Data availability.HYSPLIT is available at http://ready.arl.noaa.gov/HYSPLIT.php.All input data supporting the findings of this study are listed in the Supplement in ASCII format.

Figure 2 .
Figure 2. Number of days each year with ADS events (n ADS ) at (a) BR, (b) DJ, (c) KS, (d) UR.Modern data set includes the PM 10 observations.

Figure 3 .
Figure 3.Comparison of ADS events for the past 13 years.(a) Number of days each year with ADS events (n ADS ) in a histogram, (b) annual variations of ADS events, (c) monthly variations of ADS events.

Figure 4 .
Figure 4. Annual variations of PM 10 dust air concentration at (a) BR, (b) DJ, (c) KS, (d) UR.The lower panel shows a distribution of individual dust air concentration.The upper panel displays box plots for which the central box represents the interquartile of annual mean dust air concentration and the whisker lines are extending beyond the maximum and the minimum.In the present study, an individual ADS event was counted as the day with an advisory warning, with PM 10 dust air concentration exceeding 150 µg m −3 per hour.For comparison another reference line was marked for a more significant danger warning with PM 10 dust air concentration exceeding 400 µg m −3 for an hour.

Figure 5 .
Figure 5. Monthly variations of PM 10 dust air concentration at (a) BR, (b) DJ, (c) KS, (d) UR.The lower panel shows a distribution of individual PM 10 dust air concentration.Data symbols and reference lines are as in Fig. 4.

R
Figure 6.(a) Power-law fit to the ADS events as a function of PM 10 dust air concentration.The two reference lines are the maximum and minimum asymptotic power fitting with slopes of −4.0 and −1.5, respectively.(b) Cumulative distribution of PM 10 dust air concentration for the ADS data.Median values of PM 10 dust air concentration (PM 10median ) were 318.72 µg m −3 for BR, 162.13 µg m −3 for DJ, 164.49µg m −3 for KS, and 171.62 µg m −3 for UR.

Figure 7 .
Figure 7. Spatial distribution of ADSs on a rose diagram for (a) BR, (b) DJ, (c) KS, (d) UR.The left panel shows an angular distribution of the number of days each year with ADS events (n ADS ).The right panel shows an angular distribution of PM 10 dust air concentration.

Figure 8 .
Figure 8.(a) An inverse correlation between ADS occurrence (n ADS ) and annual mean PM 10 dust air concentration (PM 10mean ).(b) A positive correlation between median PM 10 dust air concentration (PM 10median ) and annual mean PM 10 dust air concentration (PM 10mean ).

Table 1 .
Azimuthal dependence of ADS transport path.

Table
Summary of ADS discharge from January 2003 to August 2015.

Table 3 .
Statistical significance of the increased occurrence rate of ADS events.The Welch's t test is used to test the hypothesis that two populations have equal means.The Welch's t test is an extension of Student's t test and is more reliable when the two sample sets have unequal variances and unequal sample sizes.The value of P is the probability of obtaining significantly different sample means between two data sets.

Table 4 .
Kolmogorov-Smirnov comparison of two data sets.