Long-Range Transport of SO 2 from Continental Asia to Northeast Asia and the Northwest Pacific Ocean : Flow Rate Estimation Using OMI Data , Surface in Situ Data , and the HYSPLIT Model

This present study suggests a method to calculate the SO2 flow rate from a source area to receptor areas on a regional scale using Ozone Monitoring Instrument (OMI) SO2 products, surface in situ SO2 data, and the hybrid single particle Lagrangian integrated trajectory (HYSPLIT) model. The method was implemented to calculate the SO2 flow rate from continental Asia to northeast Asia and the Northwest Pacific Ocean. For the high SO2 events when SO2 was transported from continental Asia to Japan via the Korean Peninsula on 22–24 December 2006, the long-range transported SO2 flow rates were 14.0 (21.0) Mg ̈h ́1 OMI ̈ gird ́1 at Gangneung (Seoul) in Korea and 4.2 (5.3) Mg ̈h ́1 OMI ̈gird ́1 at Hiroshima (Kumamoto) in Japan. For the long-range transport of SO2 from continental Asia to the Northwest Pacific Ocean on 6–7 October 2008 (9–11 October 2006), the flow rates were 16.1 (16.2) Mg ̈h ́1 OMI ̈gird ́1 at Hokkaido, Japan (Vladivostok, Russia) and 5.6 (16.7) Mg ̈h ́1 OMI ̈gird ́1 at the Aleutian Islands, Northwest Pacific Ocean (Bering Sea). The mean rates of decrease in the SO2 flow rate per 1000 km were also calculated between continental Asia and the receptor areas. Uncertainties in the flow rate estimates were also assessed and discussed.


Introduction
Sulfur dioxide (SO 2 ) in the atmosphere is emitted from anthropogenic and natural sources [1][2][3][4][5].More than 70% of the global total SO 2 emission is from anthropogenic sources and half of this is from fossil-fuel combustion [6].In last two decades, atmospheric SO 2 has shown a decreasing trend in regions including North America and Europe [7]; however, the global total SO 2 has increased due to increased fossil-fuel combustion in China [8,9].Although China has started showing a decreasing trend in SO 2 emission since 2008 [10], SO 2 emission in China is still reported to account for more than 60% of the total Asian SO 2 emission.Today, SO 2 emission in China is reported to account for about 25% of the global total SO 2 emission and more than 90% of the total East Asian SO 2 emission [10,11].
Chemical reactions transform SO 2 into sulfate and sulfuric acid [5].These sulfuric species may directly degrade the local air quality in receptor areas while they also cause acid rain and poor visibility at the receptor areas via long-range transport (LRT) [12][13][14].SO 2 and various sulfuric aerosols are also known to affect radiative forcing [15].Therefore, to establish effective air-quality control strategies at receptor areas, it is necessary to quantify the contribution of long-range transported air pollutants such as SO 2 to the atmospheric environment of that area.
Recent studies have examined long-range transported SO 2 (LRT-SO 2 ) emitted from continental Asia and investigated its contribution to downwind regions.Some of these studies utilized SO 2 data from satellite sensors that enable the investigation of SO 2 LRT as they provide regional and global coverage over short time intervals from one to several days.Chemical transport models and the data obtained from aircraft and satellite (Moderate Resolution Imaging Spectroradiometer (MODIS), Multi-angle Imaging Spectroradiometer (MISR)) measurements have been utilized to estimate the contribution of LRT-SO 2 from East Asia to North America [14].The LRT of SO 2 from continental Asia to the Korean Peninsula has been detected using data from the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) and validated through comparison with Multi Axis Differential Optical Absorption Spectroscopy (MAX-DOAS) and in situ measurements [16].The meteorological mechanisms of transport were reported based on aircraft observations over Northeast China [17].Lu et al. [8] reported correlations between sulfur variations in East Asia and SO 2 emissions in China using the Ozone Monitoring Instrument (OMI) and SCIAMACHY data.Jeong et al. [18] estimated the contribution of LRT-SO 2 from China to increased SO 2 concentration in Seoul using a Conditional Potential Source Contribution Function (CPSCF) with SO 2 data obtained from OMI and in situ measurements.Hsu et al. [19] reported high SO 2 events over the Pacific Ocean due to LRT from China, as inferred from OMI, Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) data, and the hybrid single particle Lagrangian integrated trajectory (HYSPLIT) model.Mallik et al. [20] reported an episode of high SO 2 concentration over Northern India as a result of LRT from Africa using OMI, CALIPSO, and backward trajectory analysis data.Li et al. [21] reported a case study of transport and evolution of a SO 2 plume from Northern China using OMI, MODIS, and backward trajectory analysis data.The above studies revealed that LRT of SO 2 and sulfuric species emitted in China and East Asia can affect not only neighboring areas but also remote regions such as the Pacific Ocean and North America.They either reported the detection of LRT-SO 2 or estimated the contribution of LRT-SO 2 to the surface SO 2 concentration at receptor areas.
For better understanding of the LRT-SO 2 contribution to changes in the SO 2 amount at a receptor area, it is necessary to identify the total amount of LRT-SO 2 at the receptor area, meaning, in turn, that we need to calculate LRT-SO 2 flow rate.The flow rate information can be directly used to quantify the total LRT-SO 2 amounts for a duration at a receptor area if it is obtained either hourly or daily by Low Earth Orbit (LEO) or geostationary orbit (GEO) satellite measurement, respectively.Here, SO 2 flow rate denotes the mass of SO 2 that passes per unit time.In the present study, we calculate the flow rate of LRT-SO 2 from continental Asia to the Korean Peninsula, Japan, and the Northwest Pacific Ocean using OMI SO 2 column data, in situ SO 2 data, and HYSPLIT backward trajectory simulations.In addition, we investigate the changes in the LRT-SO 2 flow rate as a function of distance from continental Asia.Uncertainties in the flow rate estimates are also discussed.

Data
The SO 2 column data used to calculate the flow rate of the LRT-SO 2 were obtained from the OMI onboard NASA's Earth Observing System (EOS)/Aura satellite [22] that was launched on 15 July 2004 in a sun-synchronous ascending polar orbit at 705 km altitude with 13:45 local equator crossing time.The details of the SO 2 retrieval method have been described previously [23,24].SO 2 column data in the OMSO2e level-3 products (i.e., OMSO2e at a resolution of 0.25 ˝) were acquired from NASA's Goddard Earth Sciences Data and Information Services Center (GES-DISC) [25].We used the SO 2 data produced with the updated algorithm in 2015 [26].We averaged six pixels of OMI PBL SO 2 to enhance the accuracy (NASA's Goddard Earth Sciences Data and Information Services Center (GES-DISC) [27]).We used the OMI SO 2 data for the domain 110 ˝E-180 ˝E and 20 ˝N-70 ˝N (Figure 1) for the four years of the study period (2005)(2006)(2007)(2008) when the OMI CCD detector condition was good since the amount of available SO 2 data has been reduced in recent years due to row anomaly problems in the OMI CCD detector (NASA's Goddard Earth Sciences Data and Information Services Center (GES-DISC) [28]).During the study period, the flow rates of the LRT-SO 2 were calculated only for high SO 2 events.LRT-SO 2 event dates in the present study denote the dates when the contribution of LRT-SO 2 from continental Asia is considered to significantly change the atmospheric SO 2 level at a receptor area.The method for determining the event dates is described in Section 2.3.
Atmosphere 2016, 7, 53 3 of 14 We used the OMI SO2 data for the domain 110°E-180°E and 20°N-70°N (Figure 1) for the four years of the study period (2005)(2006)(2007)(2008) when the OMI CCD detector condition was good since the amount of available SO2 data has been reduced in recent years due to row anomaly problems in the OMI CCD detector (NASA's Goddard Earth Sciences Data and Information Services Center (GES-DISC) [28]).During the study period, the flow rates of the LRT-SO2 were calculated only for high SO2 events.LRT-SO2 event dates in the present study denote the dates when the contribution of LRT-SO2 from continental Asia is considered to significantly change the atmospheric SO2 level at a receptor area.The method for determining the event dates is described in Section 2.3.Surface SO2 concentration data, which are continuously available, were used to determine the SO2 event date.We used SO2 concentration data obtained from in situ measurements at surface level on the Korean Peninsula.Figure 1 shows the locations of the in situ measurement sites at Seoul, Incheon, Gangneung, and Wonju.These four sites benefit from high availability of continuously monitored surface SO2 data and a large amount of in situ data synchronous with the OMI SO2 column data.SO2 concentrations were measured using the ultraviolet fluorescence method with a sulfur dioxide analyzer (model 4108, Dasibi Enviromental Corp., Glendale, CA; US Environmental Protection Agency reference method EQSA-1086-061).The analyzer provides SO2 concentration every 10 min.We used the hourly mean SO2 concentration.

HYSPLIT Backward Trajectory Simulation
Three-day backward trajectories from the receptor sites (Seoul, Incheon, Gangneung, and Wonju) were calculated every hour at intervals of 0.1 km above ground level using the HYSPLIT model version 4.9 developed by the National Oceanic and Atmospheric Administration, Air Resources Laboratory (NOAA/ARL) [29].Re-analysis meteorological data from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) were used for the trajectory calculation.The global reanalysis data are available on a 2.5° grid at 6 h intervals (synoptic times) on pressure surfaces.HYSPLIT backward trajectory simulations were conducted for the calculation of average air mass transportation speed between source and receptor areas.The calculation of average air mass transportation speed is described in Section 3.2.Uncertainties in the trajectory endpoints have been reported to be ~20% of the travel distance [30].A more recent study of biases in backward trajectories reported that the errors in single trajectories may, in fact, be greater than 20%, especially if the error is in the first few time steps [31].Surface SO 2 concentration data, which are continuously available, were used to determine the SO 2 event date.We used SO 2 concentration data obtained from in situ measurements at surface level on the Korean Peninsula.Figure 1 shows the locations of the in situ measurement sites at Seoul, Incheon, Gangneung, and Wonju.These four sites benefit from high availability of continuously monitored surface SO 2 data and a large amount of in situ data synchronous with the OMI SO 2 column data.SO 2 concentrations were measured using the ultraviolet fluorescence method with a sulfur dioxide analyzer (model 4108, Dasibi Enviromental Corp., Glendale, CA; US Environmental Protection Agency reference method EQSA-1086-061).The analyzer provides SO 2 concentration every 10 min.We used the hourly mean SO 2 concentration.

HYSPLIT Backward Trajectory Simulation
Three-day backward trajectories from the receptor sites (Seoul, Incheon, Gangneung, and Wonju) were calculated every hour at intervals of 0.1 km above ground level using the HYSPLIT model version 4.9 developed by the National Oceanic and Atmospheric Administration, Air Resources Laboratory (NOAA/ARL) [29].Re-analysis meteorological data from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) were used for the trajectory calculation.The global reanalysis data are available on a 2.5 ˝grid at 6 h intervals (synoptic times) on pressure surfaces.HYSPLIT backward trajectory simulations were conducted for the calculation of average air mass transportation speed between source and receptor areas.The calculation of average air mass transportation speed is described in Section 3.2.Uncertainties in the trajectory endpoints have been reported to be ~20% of the travel distance [30].A more recent study of biases in backward trajectories reported that the errors in single trajectories may, in fact, be greater than 20%, especially if the error is in the first few time steps [31].

Study Locations and Dates
The flow rates of the LRT-SO 2 were calculated for LRT-SO 2 event dates.We define LRT-SO 2 event dates as those when the SO 2 at receptor sites is greatly enhanced by SO 2 long-range transported from continental Asia.The flowchart in Figure 2 is used to determine LRT-SO 2 events on the Korean Peninsula and Japan.Continuously monitored SO 2 data are useful in quantifying the background SO 2 level and are necessary to determine the LRT-SO 2 event since the difference between the background SO 2 level and the increase in SO 2 level due to the LRT event at a receptor area is required.Therefore, we used in situ measurement data, which provide high continuity of surface SO 2 data, together with OMI data to determine LRT-SO 2 events as shown in Figure 2.
Atmosphere 2016, 7, 53 4 of 14 The flow rates of the LRT-SO2 were calculated for LRT-SO2 event dates.We define LRT-SO2 event dates as those when the SO2 at receptor sites is greatly enhanced by SO2 long-range transported from continental Asia.The flowchart in Figure 2 is used to determine LRT-SO2 events on the Korean Peninsula and Japan.Continuously monitored SO2 data are useful in quantifying the background SO2 level and are necessary to determine the LRT-SO2 event since the difference between the background SO2 level and the increase in SO2 level due to the LRT event at a receptor area is required.Therefore, we used in situ measurement data, which provide high continuity of surface SO2 data, together with OMI data to determine LRT-SO2 events as shown in Figure 2. The following conditions (Figure 2) are applied in sequence to determine the dates of LRT-SO2 events when the contribution of LRT-SO2 from continental Asia to a receptor area is considered significant.If SO2 concentration rapidly and significantly increases for a short time interval at the receptor sites, LRT-SO2 event is very likely to take place due to the almost negligible SO2 emission in Korean Peninsula in comparison with that at upwind regions [11].Thus, Conditions 1 and 2 are applied to The following conditions (Figure 2) are applied in sequence to determine the dates of LRT-SO 2 events when the contribution of LRT-SO 2 from continental Asia to a receptor area is considered significant.If SO 2 concentration rapidly and significantly increases for a short time interval at the receptor sites, LRT-SO 2 event is very likely to take place due to the almost negligible SO 2 emission in Korean Peninsula in comparison with that at upwind regions [11].Thus, Conditions 1 and 2 are applied to ensure the both rapid and significant increase in SO 2 concentration for a short time interval at the receptor sites.
A low threshold value in Condition 1 allows dates with high surface SO 2 concentration and also possibly high OMI SO 2 column.On those dates, there can be larger differences between LRT-SO 2 level and background SO 2 level, which is controlled by local SO 2 emission, than that on dates with low surface SO 2 concentration (See Section 3.1 for the detailed description of the background SO 2 level calculation).For example, on the dates that satisfy Condition 1 with the threshold of 0.05% instead of 0.01%, the average in situ mixing ratio and OMI SO 2 column are 8 ppbv and 0.011 g¨m ´2, respectively.The OMI SO 2 column of 0.011 g¨m ´2 is smaller than 0.016 g¨m ´2, the sum of background OMI SO 2 level (0.004 g¨m ´2) and its uncertainty.The uncertainty of background OMI SO 2 column in this present study is the sum of 0.001 g¨m ´2, the standard deviation of the background SO 2 and 0.011 g¨m ´2, 1σ error of OMI SO 2 column according to the literature [11].Thus, there can be a large uncertainty in quantifying LRT-SO 2 .In the case of Condition 1 with the threshold 0.03%, the OMI SO 2 column of 0.018 g¨m ´2 is slightly larger but still similar to 0.016 g¨m ´2, the sum of background OMI SO 2 level and its uncertainty.In the case of Condition 1 with the threshold 0.01%, the OMI SO 2 column of 0.033 g¨m ´2 is about two times larger than 0.016 g¨m ´2, the sum of background OMI SO 2 level and its uncertainty, which is likely to lead to a smaller uncertainty in quantifying LRT-SO 2 than that of using higher threshold values in Condition 1.
Table 1 lists the LRT-SO 2 event dates determined using the five conditions shown in Figure 2. E1 through E7 represent high SO 2 events when LRT-SO 2 reaches the Korean Peninsula, whereas E8 and E9 refer to the events when LRT-SO 2 reaches Japan via the Korean Peninsula.E10 and E11 did not meet the above five conditions, which require the surface SO 2 in situ data since air mass pathways from the source areas do not pass over the Korean Peninsula where the surface SO 2 in situ data are available.However, E10 and E11 were considered as LRT-SO 2 events in the Northwest Pacific Ocean according to Hsu et al.
[19], who reported events in which LRT-SO 2 from continental Asia reached the Northwest Pacific Ocean using both OMI data and a transport model (HYSPLIT).The cross-sectional area in Equation ( 1) can be defined as the total cross-sectional area of the atmosphere that the long-range transported air mass flows through at a receptor area, and is located perpendicular to the long-range transported air mass transport direction at the receptor area.In the present study, the cross-sectional area is obtained by multiplying the OMI grid length by the atmospheric thickness.Atmospheric thickness is the distance between the surface and the Top Of the Atmosphere (TOA).The term "atmospheric thickness" eventually has no effect on the flow rate estimation since it is canceled out in Equation (3).The average flux term can be expressed as follows: The cross-sectional area in Equation ( 1) can be defined as the total cross-sectional area of the atmosphere that the long-range transported air mass flows through at a receptor area, and is located perpendicular to the long-range transported air mass transport direction at the receptor area.

Flow Rate Calculation
In the present study, the cross-sectional area is obtained by multiplying the OMI grid length by the atmospheric thickness.Atmospheric thickness is the distance between the surface and the Top Of the Atmosphere (TOA).The term "atmospheric thickness" eventually has no effect on the flow rate estimation since it is canceled out in Equation ( 3

2)
Here, α denotes the LRT-SO 2 column.LRT-SO 2 column (α) means the SO 2 column at a receptor area due to LRT from continental Asia.β denotes the average transportation speed of air mass.The transportation speed of air mass (β) can be defined as the average local transportation speeds of the air masses, which are originated within the PBL in a source area and are also within the OMI grid at a receptor site and the air mass here is assumed to contain SO 2 that originates from the source region.Replacing the average flux term in Equation ( 1) with the terms on the right hand side of Equation ( 2 The OMI grid length is 25 km, since OMI Level 3 data were used in this study.The methods to obtain the α and β are described in Sections 3.1 and 3.2 respectively.

Calculation of α
α was calculated as follows: In Equation ( 4), SO 2 VCD local represents the SO 2 vertical column density (VCD) due to local SO 2 emission at a receptor area at an OMI overpass time for each event in Table 1.SO 2 VCD local was calculated by averaging OMI SO 2 column data that were available for the 15 days before and after each event date, excluding the data measured on the event date.SO 2 VCD event represents the OMI SO 2 column measured for each event in Table 1.Table 2 lists the α values calculated for each event and receptor site.We found the mean value (0.033 g¨m ´2) of α to be much larger than 0.008 g¨m ´2, which is the sum of the mean SO 2 VCD local (0.006 g¨m ´2) and the standard deviation of SO 2 VCD local (0.001 g¨m ´2).This large α value, compared with SO 2 VCD local plus its standard deviation, clearly distinguishes LRT-SO 2 from local SO 2 .We also calculated LRT-surface SO 2 mixing ratio using the Equation (4) but for surface mixing ratio instead of VCD.In Figure 4, LRT-SO 2 mixing ratios were calculated using the in situ data averaged for a one-hour interval between 13:30 and 14:30 (LT).The α values are plotted against LRT-SO 2 mixing ratios, showing certain agreement (R = 0.51) in Figure 4.The high R value between LRT-SO 2 column and LRT-SO 2 mixing ratio implies that the contribution of long-range transported SO 2 to a receptor area is vertically similar, whereas the low R value means that the contribution of long-range transported SO 2 to a receptor area is vertically different.The slope and intercept of the regression equation in Figure 4 do not mean much in comparison with the R value.However, there exists the y-interceptor (0.02) of the regression equation in Figure 4 due to no zero values present at xand y-axis because only the data for the SO 2 event periods can be used to calculate either LRT-SO 2 column or LRT-SO 2 mixing ratio in this present study.
Atmosphere 2016, 7, 53 7 of 14 Here, α denotes the LRT-SO2 column.LRT-SO2 column (α) means the SO2 column at a receptor area due to LRT from continental Asia.β denotes the average transportation speed of air mass.The transportation speed of air mass (β) can be defined as the average local transportation speeds of the air masses, which are originated within the PBL in a source area and are also within the OMI grid at a receptor site and the air mass here is assumed to contain SO2 that originates from the source region.Replacing the average flux term in Equation ( 1) with the terms on the right hand side of Equation ( 2), we have The OMI grid length is 25 km, since OMI Level 3 data were used in this study.The methods used to obtain the α and β are described in Sections 3.1 and 3.2, respectively.

Calculation of 𝛼
α was calculated as follows: In Equation ( 4), SO2 VCDlocal represents the SO2 vertical column density (VCD) due to local SO2 emission at a receptor area at an OMI overpass time for each event in Table 1.SO2 VCDlocal was calculated by averaging OMI SO2 column data that were available for the 15 days before and after each event date, excluding the data measured on the event date.SO2 VCDevent represents the OMI SO2 column measured for each event in Table 1.Table 2 lists the  values calculated for each event and receptor site.We found the mean value (0.033 g•m −2 ) of  to be much larger than 0.008 g•m −2 , which is the sum of the mean SO2 VCDlocal (0.006 g•m −2 ) and the standard deviation of SO2 VCDlocal (0.001 g•m −2 ).This large α value, compared with SO2 VCDlocal plus its standard deviation, clearly distinguishes LRT-SO2 from local SO2.We also calculated LRT-surface SO2 mixing ratio using the Equation (4) but for surface mixing ratio instead of VCD.In Figure 4, LRT-SO2 mixing ratios were calculated using the in situ data averaged for a one-hour interval between 13:30 and 14:30 (LT).The α values are plotted against LRT-SO2 mixing ratios, showing certain agreement (R = 0.51) in Figure 4.The high R value between LRT-SO2 column and LRT-SO2 mixing ratio implies that the contribution of long-range transported SO2 to a receptor area is vertically similar, whereas the low R value means that the contribution of long-range transported SO2 to a receptor area is vertically different.The slope and intercept of the regression equation in Figure 4 do not mean much in comparison with the R value.However, there exists the y-interceptor (0.02) of the regression equation in Figure 4 due to no zero values present at x-and y-axis because only the data for the SO2 event periods can be used to calculate either LRT-SO2 column or LRT-SO2 mixing ratio in this present study.

Calculation of β
Sulfur dioxide is likely to be concentrated within the PBL at the source area.Therefore, the air mass that starts traveling within the PBL over continental Asia during the event in Table 1 is assumed to contain SO 2 from the source area.We used the seasonal mean PBL height information reported by von Engeln and Teixeira [32] in order to reflect seasonal PBL variations at the source region on continental Asia.Air mass backward trajectories from each receptor site were computed every 0.1 km above ground level (agl) for each event, using the HYSPLIT model.The source areas, which are identified in this present study, are coincident with the areas with high SO 2 emission reported in the Asian emission inventory study [8,11,33] and also the receptor modeling study [18].β was calculated by averaging the local transportation speeds of the air masses, which are originated within the PBL in a source area and are also within the OMI grid at a receptor site. Figure 5 shows an example of the vertical distribution of the backward trajectories that started traveling between continental Asia and Incheon on 11 February 2006.Table 2 lists the β values for each event and receptor site as well as the LRT-SO 2 flow rates for each event and receptor site using Equation (3) with the α and β values in Table 2.

LRT-SO2 Flow Rates
Table 2 lists α and β, and LRT-SO2 flow rate for each event over the Korean Peninsula, Japan, and the Northwest Pacific Ocean.In events E1-E7, when LRT-SO2 reached the Korean Peninsula, the maximum flow rate is 33.4 mg•h −1 OMI•gird −1 in Wonju on 22 December while the minimum flow rate is 2.9 mg•h −1 OMI•gird −1 in Wonju on 8 November 2007.In events E8 and E9 when LRT-SO2 reaches Japan via the Korean Peninsula, the maximum flow rate is 24.2, the decreasing tendency in LRT-SO2 flow rate over transport distance is mostly due to a decreasing trend of α since β does not show any trend.The decreasing trend of α could be attributed to dispersion of SO2 over transport distance between a source area and a receptor area.In Table 2, average LRT-SO2 flow rates over the Korean Peninsula, Japan, and the Northwest Pacific Ocean are calculated to be 13.1, 8.5, and 11.2 mg•h −1 OMI•gird −1 , respectively.We have the average flow rate over the Northwest Pacific Ocean higher than that over Japan due to the contribution of the volcanic emissions, which are explained in detail later in this section.
We investigated the variations in LRT-SO2 flow rates as a function of distance from continental Asia for events E8-E11.Figure 6a shows the backward trajectories that correspond to events E8 in Table 1, while Figure 6b shows the variation of the LRT-SO2 flow rate for event E8 as a function of distance from continental Asia.In Figure 6a, A, B, and C represent the locations of the air mass backward trajectory for event E8 on 22-24 December 2006 at Shenyang, Gangneung, and Hiroshima, respectively.The corresponding SO2 flow rates, shown in Figure 6b, have values of 22.7, 14.0, and 4.2 mg•h −1 OMI•gird −1 at A, B, and C, which are 510, 1090, and 1890 km from the source area, respectively.Figure 6c,d is the same as Figure 6a,b but for E9.In Figure 6c, A, B, and C 2, the decreasing tendency in LRT-SO 2 flow rate over transport distance is mostly due to a decreasing trend of α since β does not show any trend.The decreasing trend of α could be attributed to dispersion of SO 2 over transport distance between a source area and a receptor area.In Table 2, average LRT-SO 2 flow rates over the Korean Peninsula, Japan, and the Northwest Pacific Ocean are calculated to be 13.1, 8.5, and 11.2 Mg¨h ´1 OMI¨gird ´1, respectively.We have the average flow rate over the Northwest Pacific Ocean higher than that over Japan due to the contribution of the volcanic emissions, which are explained in detail later in this section.

LRT-SO 2 Flow Rates
We investigated the variations in LRT-SO 2 flow rates as a function of distance from continental Asia for events E8-E11.Figure 6a shows the backward trajectories that correspond to events E8 in Table 1, while Figure 6b shows the variation of the LRT-SO 2 flow rate for event E8 as a function of distance from continental Asia.In Figure 6a, A, B, and C represent the locations of the air mass backward trajectory for event E8 on 22-24 December 2006 at Shenyang, Gangneung, and Hiroshima, respectively.The corresponding SO 2 flow rates, shown in Figure 6b, have values of 22.7, 14.0, and 4.2 Mg¨h ´1 OMI¨gird ´1 at A, B, and C, which are 510, 1090, and 1890 km from the source area, respectively.Figure 6c,d is the same as Figure 6a,b but for E9.In Figure 6c, A, B, and C   Figure 7a shows the backward trajectories that correspond to events E10 and E11 in Table 1, while Figure 7b shows the variation of the LRT-SO2 flow rate for events E10 and E11 as a function of distance from continental Asia.In Figure 7a, A  Figure 7a shows the backward trajectories that correspond to events E10 and E11 in Table 1, while Figure 7b shows the variation of the LRT-SO 2 flow rate for events E10 and E11 as a function of distance from continental Asia.In Figure 7a, A The values for event E11 are 16.2, 9.0, and 16.7 Mg¨h ´1 OMI¨gird ´1 at A, B, and C, which are 1230, 3820, and 5120 km from the source area, respectively.We found a general decrease in LRT-SO 2 flow rates with distance from the source region except for locations between B and C. The increase in LRT-SO 2 between B and C in Figure 7b could be attributed to the contribution of SO 2 from continuous volcanic emissions (e.g., Sheveluch, Klyuchevskov, Bezymianny, Karymsky, Ebeko, Severgin, and Berga) over the Kamchatka Peninsula, which is located on the air mass pathway between B and C in Figure 7a.The increase in LRT-SO 2 between B and C in Figure 7b is thought to be associated with contribution of the volcanic emissions since the volcanic gas emissions except for large eruption cases are generally known to affect the lower troposphere [34][35][36][37].The rate of decrease in the LRT-SO 2 flow rate between locations A' and B' is 4.3 Mg¨h ´1 OMI¨gird ´1 per 1000 km distance from the source area for event E10, whereas that between locations A and B is 4.5 Mg¨h ´1 OMI¨gird ´1 per 1000 km distance from the source area for event E11.
Atmosphere 2016, 7, 53 11 of 14 continuous volcanic emissions (e.g., Sheveluch, Klyuchevskov, Bezymianny, Karymsky, Ebeko, Severgin, and Berga) over the Kamchatka Peninsula, which is located on the air mass pathway between B and C in Figure 7a.The increase in LRT-SO2 between B and C in Figure 7b is thought to be associated with the contribution of the volcanic emissions since the volcanic gas emissions except for large eruption cases are generally known to affect the lower troposphere [34][35][36][37].The rate of decrease in the LRT-SO2 flow rate between locations A' and B' is 4.3 mg•h −1 OMI•gird −1 per 1000 km distance from the source area for event E10, whereas that between locations A and B is 4.5 mg•h −1 OMI•gird −1 per 1000 km distance from the source area for event E11.

Error Estimation
The uncertainties in LRT-SO2 flow rate were calculated following the error propagation analysis of Clifford [38]:

Error Estimation
The uncertainties in LRT-SO 2 flow rate were calculated following the error propagation analysis of Clifford [38]: Here, flow rate represents the LRT-SO 2 flow rate.σ represents the uncertainty in LRT-SO 2 flow rate, and σ α and σ β represent the uncertainties in α and β, respectively.In this present study, OMI SO 2 column uncertainty of 45% was used for σ α following Lee et al. [39].We also calculated σ β , which is composed of the covariance of the distance and time for each trajectory as well as the uncertainty in air mass transportation speed.We used 20% for the uncertainty in air mass transportation speed, following Stohl [30].The LRT-SO 2 flow rate uncertainties calculated using Equation ( 5) are shown as numbers in brackets in Table 2.The uncertainties are also indicated by the error bars in Figures 6b,d and 7b.For all events from E1 through E11, the uncertainties in LRT-SO 2 flow rate range from 51% to 61% (average, 56%).

Conclusions
We report the first use of the OMI SO 2 column and backward trajectory simulation to calculate the LRT-SO 2 flow rate from continental Asia to the Korean Peninsula, Japan, and the Northwest Pacific Ocean.In general, we found a decreasing trend in LRT-SO 2 flow rates with distance from the source region.However, the flow rate increased with distance from continental Asia in the case of event E11 when large volcanic areas in the Kamchatka Peninsula, another significant SO 2 source, are located on the air mass pathway between continental Asia and the Northwest Pacific Ocean (Bering Sea).For future studies, the method for estimation of LRT-SO 2 can be applied to other downwind regions with significant SO 2 emission sources.The method can also be utilized for other atmospheric species (e.g., CO and aerosol).As the information of a transported air pollutant amount can be available at a receptor area, this method is expected to be useful for the efficient local air quality management.In addition, the application of this method to geostationary measurements (e.g., Geostationary Environmental Monitoring Spectrometer (GEMS)) allows quantifying the hourly amount of gases or aerosols that are long-range transported from a source area to receptor areas.

Figure 1 .
Figure 1.Locations of SO2 source and receptor sites.

Figure 1 .
Figure 1.Locations of SO 2 source and receptor sites.

Figure 2 .
Figure 2. Flowchart for the determination of high SO2 event dates over the Korean Peninsula and Japan.

• Condition 1 . 2 . 3 . 4 . 5 .
Dates with the top 0.01% of surface SO2 concentrations during the study period at the four sites on the Korean Peninsula.• Condition Daily mean surface SO2 concentration is three times higher than that of the previous day.• Condition Backward trajectory calculations show that the air mass is transported from the major SO2 source areas such as Shenyang and Beijing (Figure 1) to the Korean Peninsula.• Condition Both OMI SO2 column and in situ SO2 measurement data are available.Dates that satisfy Conditions (1-4) are determined as event dates for the Korean Peninsula.• Condition Backward trajectory calculations show that the air mass is transported farther into Japan.Dates that satisfy Conditions (1-5) are determined as event dates in Japan.

Figure 2 .
Figure 2. Flowchart for the determination of high SO 2 event dates over the Korean Peninsula and Japan.

‚ Condition 1 .‚ Condition 2 .‚ Condition 3 .‚ Condition 4 .‚ Condition 5 .
Dates with the top 0.01% of surface SO 2 concentrations during the study period at the four sites on the Korean Peninsula.Daily mean surface SO 2 concentration is three times higher than that of the previous day.Backward trajectory calculations show that the air mass is transported from the major SO 2 source areas such as Shenyang and Beijing (Figure1) to the Korean Peninsula.Both OMI SO 2 column and in situ SO 2 measurement data are available.Dates that satisfy Conditions (1-4) are determined as event dates for the Korean Peninsula.Backward trajectory calculations show that the air mass is transported farther into Japan.Dates that satisfy Conditions (1-5) are determined as event dates in Japan.

Figure 3 Figure 3 1 )Figure 3 .
Figure 3 summarizes the flow rate calculation for the LRT-SO 2 event dates at each receptor location of interest.The LRT-SO 2 flow rate in the present study represents the total mass of SO 2 transported from continental Asia into the atmosphere over a receptor area of size equal to an OMI grid length per unit time.The LRT-SO 2 flow rate can be expressed in terms of the average flux and cross-sectional area: Flow rate ´Mg¨h ´1 OMI¨gird ´1¯" average flux ´Mg¨km ´2¨h ´1¯ĉ ross-sectional area ´km 2 ¯(1)

Figure 3 .
Figure 3. Flow rate of the calculation method for the LRT-SO 2 event dates at each receptor location of interest.

Figure 5 .
Figure 5. Example of the vertical distribution of three-day HYSPLIT backward trajectories for air masses arriving in Incheon on 11 February 2006 (Local Time).
5 mg•h −1 OMI•gird −1 over the North Yellow Sea on 21 December 2006 while the minimum flow rate is 4.2 mg•h −1 OMI•gird −1 at Hiroshima on 24 December 2006.In events E10 and E11, when LRT-SO2 reaches the Northwest Pacific Ocean, the maximum flow rate is 16.7 mg•h −1 OMI•gird −1 at the Northwest Pacific Ocean (Bering Sea) Hokkaido 11 October 2006 whereas the minimum flow rate is 5.6 mg•h −1 OMI•gird −1 at the Northwest Pacific Ocean (Aleutian islands) on 7 October 2008.A general decrease in LRT-SO2 flow rates with increasing distance from continental Asia can be associated with reduction in either α or β according to Equation (3).As shown in Table represent the locations of the air mass backward trajectories for event E9 on 21-23 December 2006 over the North Yellow Sea, Seoul, and Kagoshima, respectively.The corresponding SO2 flow rates for E9, shown in Figure 6d, have values of 24.5, 21.0, and 5.3 mg•h −1 OMI•gird −1 at A, B, and C, which are 610, 1280, and

Figure 5 .
Figure 5. Example of the vertical distribution of three-day HYSPLIT backward trajectories for air masses arriving in Incheon on 11 February 2006 (Local Time).

Table 2
lists α and β, and LRT-SO 2 flow rate for each event over the Korean Peninsula, Japan, and the Northwest Pacific Ocean.In events E1-E7, when LRT-SO 2 reached the Korean Peninsula, the maximum flow rate is 33.4 Mg¨h ´1 OMI¨gird ´1 in Wonju on 22 December 2006, while the minimum flow rate is 2.9 Mg¨h ´1 OMI¨gird ´1 in Wonju on 8 November 2007.In events E8 and E9 when LRT-SO 2 reaches Japan via the Korean Peninsula, the maximum flow rate is 24.5 Mg¨h ´1 OMI¨gird ´1 over the North Yellow Sea on 21 December 2006 while the minimum flow rate is 4.2 Mg¨h ´1 OMI¨gird ´1 at Hiroshima on 24 December 2006.In events E10 and E11, when LRT-SO 2 reaches the Northwest Pacific Ocean, the maximum flow rate is 16.7 Mg¨h ´1 OMI¨gird ´1 at the Northwest Pacific Ocean (Bering Sea) Hokkaido 11 October 2006 whereas the minimum flow rate is 5.6 Mg¨h ´1 OMI¨gird ´1 at the Northwest Pacific Ocean (Aleutian islands) on 7 October 2008.A general decrease in LRT-SO 2 flow rates with increasing distance from continental Asia can be associated with reduction in either α or β according to Equation (3).As shown in Table represent the locations of the air mass backward trajectories for event E9 on 21-23 December 2006 over the North Yellow Sea, Seoul, and Kagoshima, respectively.The corresponding SO 2 flow rates for E9, shown in Figure 6d, have values of 24.5, 21.0, and 5.3 Mg¨h ´1 OMI¨gird ´1 at A, B, and C, which are 610, 1280, and 2240 km from the source area, respectively.Although the distances between the receptor locations (B and C) and the source area for event E8 in Figure 6a,b are slightly different from those for event in Figure 6c,d, we found that the LRT-SO 2 flow rates at receptor sites B and C in Figure 6b are closely similar to those for B and C in Figure 6d.The similarity between the LRT-SO 2 flow rates in Figure 6b,d may be associated with similarity between the two event periods in December 2006 as well as between the air mass pathways for events E8 and E9.The rate of decrease in the LRT-SO 2 flow rate between locations A and C is 13.4 Mg¨h ´1 OMI¨gird ´1 per 1000 km distance from the source area in Figure 6b compared with 11.8 Mg¨h ´1 OMI¨gird ´1 per 1000 km distance from the source area in Figure 6d.Atmosphere 2016, 7, 53 10 of 14 2240 km from the source area, respectively.Although the distances between the receptor locations (B and C) and the source area for event E8 in Figure 6a,b are slightly different from those for event in Figure 6c,d, we found that the LRT-SO2 flow rates at receptor sites B and C in Figure 6b are closely similar to those for B and C in Figure 6d.The similarity between the LRT-SO2 flow rates in Figure 6b,d may be associated with the similarity between the two event periods in December 2006 as well as between the air mass pathways for events E8 and E9.The rate of decrease in the LRT-SO2 flow rate between locations A and C is 13.4 mg•h −1 OMI•gird −1 per 1000 km distance from the source area in Figure 6b compared with 11.8 mg•h −1 OMI•gird −1 per 1000 km distance from the source area in Figure 6d.

Figure 6 .
Figure 6.(a) air mass pathways for event E8; (b) LRT-SO2 flow rate as a function of distance from the source area.A, B, and C in (a) and (b) correspond to Shenyang, Gangneung, and Hiroshima, respectively; (c) air mass pathways for event E9; (d) LRT-SO2 flow rate as a function of distance from the source area.A, B, and C in (c) and (d) correspond to the North Yellow Sea, Seoul, and Kumamoto, respectively.
Figure7ashows the backward trajectories that correspond to events E10 and E11 in Table1, while Figure7bshows the variation of the LRT-SO2 flow rate for events E10 and E11 as a function of distance from continental Asia.In Figure7a, A' and B' represent the locations of the air mass backward trajectory for event E10 on 6-7 October 2008 at Hokkaido and the Northwest Pacific Ocean (Aleutian Islands), respectively, and A, B, and C represent those for event E11 on 9-11 October 2006 at Vladivostok, Magadan, and the Northwest Pacific Ocean (Bering Sea), respectively.The corresponding LRT-SO2 flow rates are shown in Figure 7b.The values for event E10 are calculated to be 16.1 and 5.6 mg•h −1 OMI•gird −1 at A' and B', which are 2290 and 4760 km from the source area, respectively.The values for event E11 are 16.2, 9.0, and 16.7 mg•h −1 OMI•gird −1 at A, B, and C, which are 1230, 3820, and 5120 km from the source area, respectively.We found a general decrease in LRT-SO2 flow rates with distance from the source region except for locations between B and C. The increase in LRT-SO2 between B and C in Figure 7b could be attributed to the contribution of SO2 from

Figure 6 .
Figure 6.(a) air mass pathways for event E8; (b) LRT-SO 2 flow rate as a function of distance from the source area.A, B, and C in (a) and (b) correspond to Shenyang, Gangneung, and Hiroshima, respectively; (c) air mass pathways for event E9; (d) LRT-SO 2 flow rate as a function of distance from the source area.A, B, and C in (c) and (d) correspond to the North Yellow Sea, Seoul, and Kumamoto, respectively.
Figure7ashows the backward trajectories that correspond to events E10 and E11 in Table1, while Figure7bshows the variation of the LRT-SO 2 flow rate for events E10 and E11 as a function of distance from continental Asia.In Figure7a, A' and B' represent the locations of the air mass backward trajectory for event E10 on 6-7 October 2008 at Hokkaido and the Northwest Pacific Ocean (Aleutian Islands), respectively, and A, B, and C represent those for event E11 on 9-11 October 2006 at Vladivostok, Magadan, and the Northwest Pacific Ocean (Bering Sea), respectively.The corresponding LRT-SO 2 flow rates are shown in Figure 7b.The values for event E10 are calculated to be 16.1 and 5.6 Mg¨h ´1 OMI¨gird ´1 at A' and B', which are 2290 and 4760 km from the source area, respectively.The values for event E11 are 16.2, 9.0, and 16.7 Mg¨h ´1 OMI¨gird ´1 at A, B, and C, which are 1230, 3820, and 5120 km from the source area, respectively.We found a general decrease in LRT-SO 2 flow rates

Figure 7 .
Figure 7. (a) air mass pathways for event E10 (dashed line) and E11 (solid line).Open triangles indicate Kamchatka volcanic hotspots.(b) SO2 flow rate as a function of distance from the source areas.A' and B' in (a) and (b) correspond to Hokkaido and the Northwest Pacific Ocean (Aleutian Islands), respectively.A, B, and C in (a) and (b) correspond to Vladivostok, Magadan, and the Northwest Pacific Ocean (Bering Sea), respectively.

Figure 7 .
Figure 7. (a) air mass pathways for event E10 (dashed line) and E11 (solid line).Open triangles indicate Kamchatka volcanic hotspots.(b) SO 2 flow rate as a function of distance from the source areas.A' and B' in (a) and (b) correspond to Hokkaido and the Northwest Pacific Ocean (Aleutian Islands), respectively.A, B, and C in (a) and (b) correspond to Vladivostok, Magadan, and the Northwest Pacific Ocean (Bering Sea), respectively.

Table 1 .
LRT-SO 2 events determined for each receptor area.

Table 2 .
Estimated LRT-SO 2 flow rate and the parameters (α and β) used to estimate the long-range transported flow rate for each event.