GIST-PM-Asia v1: development of a numerical system to improve particulate matter forecasts in South Korea using geostationary satellite-retrieved aerosol optical data over Northeast Asia

Abstract

To improve short-term particulate matter (PM) forecasts in South Korea, the initial distribution of PM composition, particularly over the upwind regions, is primarily important. To prepare the initial PM composition, the aerosol optical depth (AOD) data retrieved from a geostationary equatorial orbit (GEO) satellite sensor, GOCI (Geostationary Ocean Color Imager) which covers a part of Northeast Asia (113–146${}^{\circ }$ E; 25–47${}^{\circ }$ N), were used. Although GOCI can provide a higher number of AOD data in a semicontinuous manner than low Earth orbit (LEO) satellite sensors, it still has a serious limitation in that the AOD data are not available at cloud pixels and over high-reflectance areas, such as desert and snow-covered regions. To overcome this limitation, a spatiotemporal-kriging (STK) method was used to better prepare the initial AOD distributions that were converted into the PM composition over Northeast Asia. One of the largest advantages in using the STK method in this study is that more observed AOD data can be used to prepare the best initial AOD fields compared with other methods that use single frame of observation data around the time of initialization. It is demonstrated in this study that the short-term PM forecast system developed with the application of the STK method can greatly improve PM${}_{10}$ predictions in the Seoul metropolitan area (SMA) when evaluated with ground-based observations. For example, errors and biases of PM${}_{10}$ predictions decreased by $\sim$ 60 and $\sim$ 70%, respectively, during the first 6 h of short-term PM forecasting, compared with those without the initial PM composition. In addition, the influences of several factors on the performances of the short-term PM forecast were explored in this study. The influences of the choices of the control variables on the PM chemical composition were also investigated with the composition data measured via PILS-IC (particle-into-liquid sampler coupled with ion chromatography) and low air-volume sample instruments at a site near Seoul. To improve the overall performances of the short-term PM forecast system, several future research directions were also discussed and suggested.

Introduction

It has been reported that there is a strong relationship between exposure to atmospheric particulate matter (PM) and human health (Brook et al., 2010; Brunekreef and Holgate, 2002; Pope and Dockery, 2006). PM has become a primary concern around the world, particularly in East Asia, where high PM pollution episodes have occurred frequently, mainly due to the large amounts of pollutant emissions from energetic economic activities. In an effort to understand the behaviors and characteristics of PM in East Asia, chemistry-transport models (CTMs) have played an important role in overcoming the spatial and temporal limitations of observations, and also enable policy makers to establish scientific implementation plans via atmospheric regulations and policies. To improve the performance of the PM simulations, integrated air quality modeling systems that consist of CTMs, meteorological models, emissions, and data assimilation using ground- and satellite-borne measurements have been introduced (Al-Saadi et al., 2005; Park et al., 2011; Song et al., 2008). However, accurate simulations of PM distributions with CTMs have been challenging, because of many uncertainties from emission fluxes, meteorological fields, and chemical and physical parameterizations in the CTMs. For example, the Korean Ministry of Environment (MoE) has recently started to implement air quality forecasts for PM${}_{10}$, PM${}_{2.5}$ and ozone over the Seoul metropolitan area (SMA), the largest metropolitan area in South Korea. However, the forecasting accuracy for high PM${}_{10}$ alert (81–120 $\mathrm{\mu }$g m${}^{-\mathrm{3}})$ in the current system has been low ($<$ $\sim$ 60 %) since 2013. Thus, urgent improvements in the PM${}_{10}$ predictions are necessary.

General structure of (a) numerical weather prediction (NWP), (b) conventional chemical weather forecast CWF), and (c) advanced chemical weather forecast system.

[Figure omitted. See PDF]

In this context, an improved short-term PM forecast system was developed and introduced, based on an analogy to the system of numerical weather prediction (NWP). Figure 1a presents a flow diagram of an NWP in which regional meteorological modeling is conducted using two important inputs: (i) boundary conditions (BCs) from global meteorological models and (ii) initial conditions (ICs) prepared via data assimilation using ground-measured data and balloon-, ship-, aircraft-, and/or satellite-borne measurements. In contrast, a conventional chemical weather forecast (CWF) (e.g., forecasts for ozone and PM) has been carried out only using meteorological fields and pollutant emissions (Fig. 1b). In the short-term PM forecast system proposed here (Fig. 1c), one more input is added to the conventional CWF system: the initial distribution of PM composition. To prepare the initial PM composition, a scheme that uses geostationary satellite-derived aerosol optical depths (AODs), is developed in this study. Similarly, the BCs for the CTM runs are obtained from global CTM simulations.

In the improved CWF system, AOD data retrieved from low Earth orbit (LEO) satellite sensors, such as Moderate Resolution Imaging Spectroradiometer (MODIS) and Multi-angle Imaging SpectroRadiometer (MISR) can be used to set up the ICs for the short-term PM forecast (Benedetti et al., 2009; Liu et al., 2011; Saide et al., 2013). While these AOD data have an advantage in spatial coverage compared with those obtained from point stations, the use of the LEO satellite-derived AODs has another limitation in acquiring continuous observations over a certain area due to the capabilities of the LEO sensors in their orbital periods and viewing swath widths.

Such limitations in using LEO satellite observations can be overcome with the help of geostationary equatorial orbit (GEO) satellite sensors providing semicontinuous observations over a specific part of the Earth during the day (Fishman et al., 2012; Lahoz et al., 2011; Zoogman et al., 2014). Recently, aerosol optical properties (AOPs) from the Geostationary Ocean Color Imager (GOCI) have become available. GOCI is the first multi-spectral ocean color sensor onboard the Communication, Ocean, and Meteorological Satellite (COMS), launched over a part of Northeast Asia (113–146${}^{\circ }$ E; 25–47${}^{\circ }$ N) in 2010, providing semicontinuous AOD, single scattering albedo (SSA), and fine-mode fraction (FMF) over a domain of Northeast Asia (Lee et al., 2010). With GOCI AOD data, a novel approach was developed to investigate transboundary PM pollution over Northeast Asia (M. E. Park et al., 2014).

In this study, we carried out hindcast studies (forecast studies with past data) to find the “best” method to improve the performance of the short-term PM forecasting using the GOCI AODs. To do this, we developed a model, Geostatistical Interpolation of Spatio-Temporal data for PM forecasting over Northeast Asia (GIST-PM-Asia) v1 that includes (i) a spatiotemporal-kriging (STK) method to spatiotemporally combine the GOCI-derived AODs, (ii) “observation operators” to convert the CTM-simulated PM composition into AODs and vice versa, and (iii) selection of “control variables” (CVs) through which the distribution of AODs can be converted back into the distributions of the PM composition to be used as the ICs. The uses of the STK method, observation operators, and CVs are illustrated in Fig. 1. The main advantages of using the STK method are discussed in detail in the main text. Several sensitivity studies were also conducted to improve the understanding of forecasting errors and biases in the short-term PM forecasting system developed.

With these research objectives and methodology, this paper is organized as follows: the hindcast framework is first described in detail in Sect. 2. In Sect. 3, the hindcast results with various configurations are evaluated with ground-based observations during the high PM episodes in SMA to find the “best” configuration for future short-term PM forecasts. Finally, a summary and conclusions are provided in Sect. 4.

Methodology

The initial aerosol composition was prepared using AOD data from both the GOCI sensor and CTM model simulations. For the CTM simulations, the Community Multi-scale Air Quality (CMAQ; v5.0.1) model (Byun and Ching, 1999; Byun and Schere, 2006) and the Weather Research and Forecast model (WRF; v3.5.1) (Skamarock and Klemp, 2008) were used. The STK method and 12 different combinations of observation operators and CVs were also used for preparing the distributions of the 3-D PM composition over the GOCI-covered domain. The CMAQ model simulations with the 12 different configurations were carried out and the performances were then tested against ground-measured AOD, PM${}_{10}$, and PM${}_{2.5}$ composition. The details of these components are described in the following sections.

WRF and CMAQ model configurations.

WRF (v3.5.1)		CMAQ (v5.0.1)
Microphysics scheme	WRF single-moment 3 class	Chemical mechanism	SAPRC-99
Long- and short-waveradiation	Rapid Radiation Transfer Modelfor GCMs (RRTMG)	Aerosol module	AERO-6
Planetary boundarylayer	Yonsei University scheme	Chemistry solver	Euler backward iterative (EBI) solver
Land-surface model	Noah-MP	Photolysis module	In-line photolysis calculations

Domains of CMAQ model simulations (black), GOCI sensor coverage (blue), and the Seoul metropolitan area (red). Also shown are seven AERONET level-2 sites (circles), 58 NAMIS PM${}_{10}$ sites (crosses), and a PM composition observation site (triangle) in the greater Seoul area, respectively.

[Figure omitted. See PDF]

Meteorological and chemistry-transport modeling

The WRF model provided meteorological data with 15 km $\times$ 15 km horizontal grid spacing and 26 vertical layers extending up to 50 hPa. To obtain highly resolved terrestrial input data, the topography height from the NASA Shuttle Radar Topography Mission (SRTM) 3 arcsec database (http://dds.cr.usgs.gov/srtm/version2_1/SRTM3) and the land use information provided by Environmental Geographic Information Service (EGIS; http://egis.me.go.kr) were used. Initial and boundary meteorological conditions for the WRF simulation were provided by the National Centers for Environmental Protection (NCEP) final operational global tropospheric analyses (http://rda.ucar.edu/datasets/ds083.2). To improve 3-D temperature, winds and water vapor mixing, objective analysis was employed by incorporating the NCEP ADP Global surface and upper air observation data. The meteorological fields were provided with 1 h temporal resolution and were then converted into the input fields for the CMAQ model simulations by the Meteorology–Chemistry Interface Processor (MCIP; v4.1) (Otte and Pleim, 2010).

The CMAQ model is a chemistry-transport model that simulates the chemical fates and transport of gaseous and particulate pollutants. In this study, the CMAQ modeling covered Northeast Asia, from 92 to 149${}^{\circ }$ E and 17 to 48${}^{\circ }$ N, using 15 km $\times$ 15 km horizontal grid spacing (Fig. 2) with 14 terrain following $\mathit{\sigma }$ coordinates, from 1000 to 94 hPa. The configurations of the WRF model and CMAQ simulation used in this study are described in Table 1.

Anthropogenic emission inputs were processed by Sparse Matrix Operator Kernel Emissions in Asia (SMOKE-Asia; v1.2.1), which has been developed for processing anthropogenic emissions for Asia. Details of SMOKE-Asia were described in Woo et al. (2012). Biogenic emissions were prepared using the Model of Emission of Gases and Aerosol from Nature (MEGAN; v2.0.4) (Guenther et al., 2006) with the MODIS-derived leaf area index (Myneni et al., 2002), MODIS land-cover data sets (Friedl et al., 2002), and the meteorological input data described above. For the consideration of biomass burning emissions, daily fire estimates provided by Fire Inventory from NCAR (FINN) were used (Wiedinmyer et al., 2011). Asian mineral dust emissions were not considered in this study. Thus, the periods for model evaluation were selected during periods when mineral dust events did not take place.

To take full advantage of the AOD data sets intensively measured during the Distributed Regional Aerosol Gridded Observation Network in Asia (DRAGON-Asia) campaign, modeling episodes were chosen for the campaign period from 1 March to 31 May 2012. First, background CMAQ model simulations were conducted for the 3-month DRAGON period with 10-day spin-up modeling. After this, initial conditions were prepared using the STK method, observation operators and CVs via the combination of GOCI AODs with the background modeling AOD. Analysis was carried out for 12 h from 12:00 in local time (LT) on 10 selected high PM pollution days. The hindcast hours are referred to as H $+$ 0 to H $+$ 12. In this study we paid more attention to the performance of the first 12 h PM${}_{10}$ hindcast results; the analysis of the hindcast results after 13 h is also discussed briefly in Sect. 3.3.

In the hindcast analysis, different hindcast runs with 12 combinations of different observation operators and CVs were conducted, as discussed in Sect. 2.4 and 2.5. We selected one episode from March (28 March), five episodes from April (8, 9, 14, 17, and 23 April), and four episodes from May (6, 13, 15, and 16 May), all in 2012, for the analysis associated with three criteria of (i) on the selected days the average PM${}_{10}$ from 12:00 to 18:00 LT was above 70 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ over SMA; (ii) on the selected days, the daily coverage of the GOCI AOD data was at least 20 % over the GOCI domain; and (iii) on the selected days, dust events were not recorded over South Korea according to the Korea Meteorological Administration (KMA). Additional hindcast runs were also conducted from 7 March 12:00 to 19 March 11:00 LT for evaluating the performances of the hindcast runs for less-polluted episodes. In this study, we focused on SMA because we were particularly interested in this area. However, the system introduced here can be applied to other areas inside the GOCI domain where surface PM observation data are available.

Observation data

GOCI AOD

As mentioned previously, GEO satellite sensors have important advantages compared with LEO satellite sensors, such as semicontinuously (with 1 h intervals) producing AOP data over a specific domain of interest. Despite this temporal advantage, it has been difficult for most GEO satellite sensors to produce accurate AOPs, because they have only one or two visible channels. In contrast, the GOCI instrument has six visible and two near-infrared channels, and can produce multi-spectral images eight times per day with a spatial resolution of approximately 500 m $\times$ 500 m with coverage of 2500 km $\times$ 2500 km, including part of Northeast China, the Korean peninsula, and Japan (Fig. 2). Using the 1 h resolved multi-spectral radiance data from GOCI, the uncertainties of AOP retrievals can be dramatically reduced (M. E. Park et al., 2014). The GOCI AOPs were retrieved with multi-channel algorithms that can provide hourly AOP data including AOD, FMF, and SSA at 550 nm (Choi et al., 2015). Compared with the algorithms from two previous studies (Lee et al., 2010, 2012), the GloA2 algorithm uses an improved lookup table for retrieving the AOPs, using extensive observations from the Aerosol Robotic Network (AERONET) and monthly surface reflectance observed from GOCI, and provides 1 h resolved AOP data at eight fixed times per day (from 09:30 to 16:30 LT) with 6 km $\times$ 6 km spatial resolution. In this study, the AOD data from the GOCI AOPs were used (because the SSA and FMF data need further improvements) and also compared with collection 5.1 10 km MODIS aerosol products from the Aqua and Terra satellites (Levy et al., 2007; Remer et al., 2005) and collection 6 3 km MODIS aerosol products from the Aqua and Terra satellites (Munchak et al., 2013) to present the relative performances of GOCI AOD. The AERONET AOD data were also used for assessing the relative accuracy of the GOCI AODs. Figure 3a, b, and c show the scatter plot analyses of three satellite-retrieved 10 km MODIS AODs, 3 km MODIS AODs, and GOCI AODs vs. AERONET level 2 AODs over the GOCI domain during the DRAGON-Asia campaign. All the satellite data were sampled within spatial and temporal differences of 3 km and 10 min from the AERONET observations. It should also be noted that the GOCI and MODIS data were compared with the AERONET data without the application of the kriging method. First, it was found that GOCI provided more frequent AOD data ($N$ $=$ 2276) than 3 km MODIS ($N$ $=$ 629) and that GOCI AODs data show comparable regression coefficient ($R$ $=$ 0.85), root mean square error (RMSE $=$ 0.25), and mean bias (MB $=$ $-$0.19), compared with 3km MODIS data ($R$ $=$ 0.89; RMSE $=$ 0.16; MB $=$ 0.06). This indicates that the GOCI AOD data not only have comparable quality to the MODIS AOD data but also provide a higher number of data over the GOCI domain. In Fig. 3d, the daily spatial AOD percent coverages of the Aqua/Terra MODIS and GOCI sensors are compared. It was found that there are a large number of daily missing pixels in the observations of both satellite sensors (the average percent coverages of Aqua MODIS, Terra MODIS and GOCI AODs during the period were about 9, 10, and 29 %, respectively).

Scatter plots of (a) 10 km Aqua/Terra MODIS AODs vs. AERONET level-2 AODs, (b) 3 km Aqua/Terra MODIS AODs vs. AERONET level-2 AODs and (c) GOCI AODs vs. AERONET level-2 AODs at 550 nm during the DRAGON campaign over the GOCI domain. $N$, $R$, RMSE, and MB represent the number of observations, the regression coefficient, root mean square error, and mean bias, respectively. Hourly resolved Aqua/Terra MODIS and GOCI spatial coverages (%) are also shown in (d) from 1 March to 31 May 2012.

[Figure omitted. See PDF]

Ground-based observations

AERONET is a global ground-based sunphotometer network managed by the NASA Goddard Space Flight Center, providing spectral AOPs including AOD, SSA, and particle size distributions, available at http://aeronet.gsfc.nasa.gov (Holben et al., 1998). To match the wavelength of GOCI AOD with AERONET AOD, the AOD data at 550 nm were calculated via interpolation, using AODs and Ångström exponent data between 440 and 870 nm from the DRAGON-Asia level 2.0 data. AOD data from 29 AERONET sites inside the GOCI domain were used for validating GOCI and STK AOD products, and those from six AERONET sites in SMA were selected for evaluating the performance of hindcast AODs.

To analyze hindcast surface aerosol concentrations, the PM${}_{10}$ observations provided by the National Ambient Air Monitoring System (NAMIS) network in South Korea were used. The NAMIS network, operated by the MoE has collected air pollutant concentrations of PM${}_{10}$ measured by an automatic $\mathit{\beta }$-ray absorption method with a detection limit of 2 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ at 5 min intervals. We selected 58 NAMIS sites in SMA, the locations of which are shown in Fig. 2, and used 1 h averaged data for the analysis during the selected episodes.

Ion concentrations of PM${}_{2.5}$ were also measured using a particle-into-liquid sampler coupled with ion chromatography (PILS-IC) and a low air-volume sampler with a Teflon filter in Yongin City, located downwind of Seoul (Fig. 2). Details on the measurement methods are described in Lee et al. (2015) and are not repeated here. The 1 h averaged sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-})$, nitrate (NO${}_{\mathrm{3}}^{-})$, and ammonium (NH${}_{\mathrm{4}}^{+})$ concentrations, measured by the PILS-IC, and 24 h averaged SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, organic carbon (OC), and elementary carbon (EC), measured by the low air-volume sampler, were used for further comparison during the selected episodes (Sect. 3.4). The observed OC concentrations were multiplied by a factor of 1.5 to estimate organic aerosol (OA) concentrations (He et al., 2011; Huang et al., 2010).

Spatiotemporal kriging

Kriging is a geostatistical interpolation method to estimate unmeasured variables and their uncertainties, using correlation structure of measured variables. An atmospheric application study of the kriging method to estimating PM${}_{10}$ exceedance days over Europe reported that ST kriging showed comparable performances to those of the EnKF (ensemble Kalman filter) approach (Denby et al., 2008).

In this study, the STK method was used to fill out the missing pixels (Fig. 3d) with the spatial and temporal GOCI AOD data. The AOD fields produced by ST kriging can be prepared with a horizontal resolution of 15 km $\times$ 15 km from 10:00 to 16:00 LT over the GOCI domain. In this study, the AOD data at 12:00 LT (H $+$ 0) during the selected episode days were used for preparing the initial conditions. The details and general application of the STK method are presented in Appendix A. One advantage of using ST kriging in this study framework is to use large numbers of observational data (GOCI AODs), compared with other methods. In fact, the GOCI AOD data are densely available temporally (with 1 h intervals) and spatially (compared with MODIS AODs; see Fig. 3a and b). This was the primary reason for using the STK method in this study. For example, when initial AOD fields were prepared at a certain time (e.g., at noon, 12:00 LT: H $+$ 0), the STK method used not only GOCI AOD data at 11:30 LT or 12:30 LT but also GOCI AOD data at 09:30, 10:30, and 13:30 LT, unlike other methods. In the case of 4 April 2012 (a high PM pollution episode during the DRAGON-Asia campaign), other interpolation methods (e.g., Cressman, bilinear, and nearest-neighbor methods) could use only the GOCI AOD data of $\sim$ 88 000 for the preparation of the initial AOD field at 12:00 LT, whereas the STK method used the GOCI AOD data of $\sim$ 280 000 (3 times more AOD data). Sequential data assimilation (DA) methods such as OI (optimal interpolation) and 3DVAR (three-dimensional variational data assimilation) can use the same number of observations as the STK method. However, they required four data assimilation steps (i.e., 4 h time window for DA) (Tang et al., 2015) to include observations from 09:30 to 13:30, thus greatly increasing the computational cost for daily assimilation.

If the observation data are densely available and the differences between the observations and model-simulated data are large (i.e., the model simulations include relatively large errors and biases), there is less “practical need” to use the CTM-simulated data in the process of data assimilation. That is, it would be more desirable if the values of the unobserved (missing) pixels could be filled in based on “more reliable” observation data (here, GOCI AODs). This would be particularly true, when the CTM-predicted AODs are systematically underestimated compared with GOCI or AERONET AODs (as will be shown in Fig. 5a). Additionally, computation costs of the STK method are so low that the STK AOD can be calculated rapidly. For example, the 1-day process for preparing the AOD fields over the GOCI domain takes only $\sim$ 20 min with two 3.47 GHz Xeon X5690 six-core processors and 32 GB of memory in the current application of the STK method. Thus, it can be applied directly to the daily CWF due to the relatively cheap computation cost. Again, computation time (rapid calculation) is a central issue in daily (short-term) chemical weather forecasts. The calculation of daily three-dimensional semivariograms takes most of the computation time (regarding the details of calculation of the daily three-dimensional semivariogram, refer to Appendix A and Fig. A1).

Connected with these discussions, in the application of the STK method to the GOCI AODs, the “optimal number” of observation data is necessary to balance the accuracy of the data and the computational speed. From many sensitivity tests (not shown here), the optimal number of observations for most missing (white) pixels is approximately 100. That is, the use of more observation data above this optimum number does not meaningfully enhance the accuracy of AODs of the missing pixels but simply takes more computation time. This number of observation data is usually available for most of the missing (white) pixels of the GOCI scenes from nearby grids both/either at the concurrent scene spatially within $\sim$ 100 km and/or at the temporally close snapshots within 3 h. Based on these reasons, the STK method was chosen for this study.

Observation operator

An observation operator (or forward operator) describes the relation between observation data and model parameters. For example, the observation operator in this study converts the aerosol composition into AODs (and vice versa). Based on the aerosol composition and the relative humidity (RH) from the model simulations, simulated AODs at a wavelength of 550 nm (${\mathit{\tau }}_{\text{CMAQ}})$ were calculated with the following observation operator: $${\mathit{\tau }}_{\text{CMAQ}}=\sum _{s=\mathrm{1}}^N\sum _{l=\mathrm{1}}^M{\mathit{\alpha }}_{s,\text{dry}}{f}_s\left({\text{RH}}_l\right)[C{]}_{s,l}{H}_l,$$ where $N$ and $M$ denote the number of aerosol species ($s)$ and model layer ($l)$, respectively, ${\mathit{\alpha }}_{s,\text{dry}}$ the mass extinction efficiency (MEE) of the species ($s)$ at 550 nm under the dry condition, $f_s$(RH${}_l)$ the hygroscopic enhancement factor for the species ($s)$ as a function of RH at the layer of $l$, [$C$]${}_{s,l}$ the mass concentration of the species ($s)$ at the layer of $l$, and $H_l$ the height of layer $l$. Here, [$C$]${}_{s,l}$ is selected as the control variable (refer to Sect. 2.5).

Values used in observation operators for estimating aerosol optical properties (AOPs).

Method for estimating	Aerosol	Hygroscopic	${\mathit{\alpha }}_{\text{OC}}^{\text{a}}$	${\mathit{\alpha }}_{\text{BC}}^{\text{b}}$	${\mathit{\alpha }}_{\text{SSAM}}^{\text{c}}$	${\mathit{\alpha }}_{\text{SSCM}}^{\text{d}}$
aerosol optical properties	speciation	aerosols
Chin et al. (2002)	(NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, OC, BC, dust (7 size bins), sea salt (2 modes)	(NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, OC, BC, sea salt	2.67	9.28	1.15	0.13
Martin et al. (2003)	(NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, OC, BC, dust (7 size bins), sea salt (2 modes)	(NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, OC, BC, sea salt	2.82	8.05	2.37	0.94
Malm and Hand (2007)	NH${}_{\mathrm{4}}$NO${}_{\mathrm{3}}$, (NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, organic matter, soil, coarse mass, sea salt	NH${}_{\mathrm{4}}$NO${}_{\mathrm{3}}$, (NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, sea salt	4.00	10.00	1.37	1.37

Dry mass extinction efficiencies (m${}^{\mathrm{2}}$ g${}^{-\mathrm{1}})$ at 550 nm of ${}^{\text{a}}$ OC, ${}^{\text{b}}$ BC, ${}^{\text{c}}$ sea salt in accumulation mode and ${}^{\text{d}}$ sea salt in coarse mode. Note: in cases of Chin et al. (2002) and Martin et al. (2003), the AOPs for sulfate were used for calculating AOPs for NH${}_{\mathrm{4}}$NO${}_{\mathrm{3}}$ and (NH${}_{\mathrm{4}}{)}_{\mathrm{2}}$SO${}_{\mathrm{4}}$.

Mass extinction efficiencies (MEEs) calculated for (a) SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$, (b) OAs, (c) BC and (d) sea salt at a wavelength of 550 nm as a function of RH (%) from three observation operators. For the GOCART and GEOS-Chem operators, 50 % of OAs and 20 % of BC are assumed to be hydrophilic. In sea-salt MEEs, the accumulate and coarse modes are represented as solid and dashed lines, respectively.

[Figure omitted. See PDF]

In this study, three observation operators were used for calculating AODs and updating initial PM composition for the hindcast studies. The differences in the observation operators are caused mainly by the differences in ${\mathit{\alpha }}_{s,\text{dry}}$ and $f_s$(RH${}_l)$ of Eq. (1). The first observation operator was selected from the Goddard Chemistry Aerosol Radiation and Transport (GOCART) model (Chin et al., 2002; hereafter GOCART operator). Hygroscopic growth rates for SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, BC, and sea-salt aerosols were considered separately in this operator. The second observation operator was from the GEOS-Chem model (the GEOS-Chem operator). The detailed aerosol speciation and MEE values were described in Martin et al. (2003). The final observation operator is based on the study of Malm and Hand (2007) (the IMPROVE operator). This observation operator was based on the reconstruction method with the MEEs and hygroscopic enhancement factors at 550 nm for different types of aerosol species. Table 2 summarizes the characteristics of the three observation operators chosen in this study. To consistently consider the characteristics of the three observation operators, aerosol types ($s$ in Eq. 1) were classified into seven groups: SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, OAs, BC, sea salt, and others, which mainly consist of PM${}_{2.5}$ trace elements (Reff et al., 2009). In the classification, internal mixing states of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ were assumed. It should also be noted that the consideration of NO${}_{\mathrm{3}}^{-}$ is important to correctly estimate AOD and aerosol mass loading in East Asia (R. S. Park et al., 2011, 2014; Song et al., 2008). Figure 4 shows the wet MEE values (${\mathit{\alpha }}_{s,\text{wet}}$; product of ${\mathit{\alpha }}_{s,\text{dry}}$ and $f_s$(RH${}_l)$ in Eq. 1) calculated for SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$, OAs, BC and sea salt at a wavelength of 550 nm as a function of RH, indicating that the three different operators can create large differences in the wet MEE values.

Definition of model configurations.

Configuration	Observation operator	Control variable
A1	Chin et al. (2002)	Total aerosol mass concentration
A2		SO${}_{\mathrm{4}}^{\mathrm{2}-}$ mass concentration
A3		SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs mass concentration
A4		SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$ and OAs mass concentration
B1	Martin et al. (2003)	Total aerosol mass concentration
B2		SO${}_{\mathrm{4}}^{\mathrm{2}-}$ mass concentration
B3		SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs mass concentration
B4		SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$ and OAs mass concentration
C1	Malm and Hand (2007)	Total aerosol mass concentration
C2		SO${}_{\mathrm{4}}^{\mathrm{2}-}$ mass concentration
C3		SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs mass concentration
C4		SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$ and OAs mass concentration

Selection of control variables

To prepare the distributions of the aerosol composition, the STK AOD fields should be converted into the 3-D aerosol composition. To do this, the differences between the STK AODs and background AODs (often called “observational increments”: $\mathrm{\Delta }$AOD${}_k$ $=$ AOD${}_{\text{STK},k}$ $-$ AOD${}_{\text{bg},k}$; $k$: grid cell) should be added to the background model-derived aerosol composition at each grid cell in connection with the observation operators (Eq. 1). Which aerosol species is/are selected for allocating $\mathrm{\Delta }$AOD${}_k$? We selected four types of control variables (CVs) of particulate species. First, all particulate species were selected as CVs. In this case, $\mathrm{\Delta }$AOD${}_k$ was distributed to all particulate species, with the particulate fractions calculated from the background CMAQ model simulations. The second CV was the selection of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentration. Despite the large contribution of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ to both AOD and PM concentration in East Asia, model-estimated SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations have shown large systematic underestimations, compared with observed SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (R. S. Park et al., 2011, 2014). This can be related to either (or both) the uncertainty in SO${}_{\mathrm{2}}$ emissions in East Asia or (and) the uncertainty in the parameterizations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production in the CTM models (Kim et al., 2013; Lu et al., 2010; Smith et al., 2011; R. S. Park et al., 2014). In addition, there is also large uncertainty in the levels of hydroxyl radicals (OH) due to uncertain daytime HONO chemistry, OH reactivation, in-plume processes and others (Archibald et al., 2010; Han et al., 2015; Karamchandani et al., 2000; Kim et al., 2009; Kubistin et al., 2010; Lelieveld et al., 2008; Song et al., 2003, 2010; Sörgel et al., 2011; Stemmler et al., 2006; Zhou et al., 2011). Obviously, these uncertainties can influence the levels of H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ and thus particulate sulfate concentrations in the atmosphere. In this case, aerosol mass concentrations (except for SO${}_{\mathrm{4}}^{\mathrm{2}-})$ were the same as those of the background aerosol concentrations. Third, SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs were chosen to be changed. Although OAs are one of the major particulate species, it is well known that OA concentrations are also systematically underestimated due to two reasons: (i) the uncertainty in the parameterizations of the secondary OA formation (Donahue et al., 2006, 2011; Dzepina et al., 2009; Hodzic et al., 2010; Matsui et al., 2014; Slowik et al., 2010), and (ii) the uncertainty in emission inventories for anthropogenic and biogenic OA precursors (Guenther et al., 1999; Han et al., 2013; Sakulyanontvittaya et al., 2008; Tsimpidi et al., 2010; Wyat Appel et al., 2008). In this case, the mass concentration of surface OAs is assumed to be equal to the mass concentration of surface SO${}_{\mathrm{4}}^{\mathrm{2}-}$, based on the ground-based measurement studies over East Asia (Lee et al., 2009; Zhang et al., 2007, 2012). Thus, $\mathrm{\Delta }$AOD${}_k$ accounted for the increments of concentrations from OAs and SO${}_{\mathrm{4}}^{\mathrm{2}-}$, which are changed independently from the background concentrations. Finally, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and OAs were selected to be changed. In this case, $\mathrm{\Delta }$AOD${}_k$ was distributed to the four species selected, with the fractions of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ calculated from background simulations. The method to change the OA concentration in the fourth selection of CVs was the same as the method in the third selection of CVs. The fourth selection of CVs was also made to consider the thermodynamic balance among SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ concentrations (Bassett and Seinfeld, 1983; Saxena et al., 1986; Seinfeld and Pandis, 2012; Song and Carmichael, 1999; Stelson et al., 1984). It should be noted that background modeling-derived vertical profiles and the size distributions of aerosol species were used for converting 2-D AOD to 3-D PM composition in all the STK cases. With the combinations of the three different observation operators and four choices of CVs (Table 3), 12 hindcast runs were made for high PM episodes during the DRAGON-Asia campaign.

Results and discussion

In Sect. 3, the performances of the STK method are evaluated via comparisons with the AERONET AOD in the GOCI domain (Sect. 3.1). Sensitivity analyses were then conducted to examine the impacts of the observation operators and CVs on the accuracy of the hindcast runs (Sect. 3.2). After that, the overall performances of the hindcasts were evaluated with ground-based observations during the high PM${}_{10}$ episodes over SMA (Sect. 3.3). A comparative analysis of the PM composition between hindcast results and observations was also conducted to further investigate/analyze the performance of the hindcast system (Sect. 3.4). In addition, hindcast results for the periods of less-polluted episodes are also shown with the best configuration (Sect. 3.5).

Evaluation of STK AODs

Figure 5a–c show scatter plot analyses of background CMAQ-simulated AODs, spatial kriging AODs (i.e., kriging only with the GOCI AODs from one scene) and STK AODs vs. AERONET level 2 AODs over the GOCI domain during the DRAGON-Asia campaign. First, it can be found that the CMAQ-predicted AODs are underestimated significantly compared with the AERONET AODs. As discussed in Sect. 2.3, this was the main reason that we used the STK method in this study. More weight should be given to observations because the CTM modeling produces significant biases. Second, STK AODs show improved correlations, compared with the AODs estimated via the spatial kriging method. Also, the STK AOD data show equivalent levels of errors and biases, compared with GOCI AOD data. If one compares Fig. 3b with Fig. 5c, it can be seen that the ST kriging can effectively produce the AOD fields (also note the increase in $N$).

Scatter plots of (a) background CMAQ model AODs, (b) spatial kriging AODs, and (c) STK AODs vs. AERONET level-2 AODs at 550 nm. Plots of ST kriging with kriging variances (KVs) less than or equal to 0.04 (d) and larger than 0.04 (e) are also shown. The color scale shown in (e) and presents the KVs of STK AODs. The number of spatial kriging AODs in (b) is smaller than those of (a) and (c) due to the missing hourly AOD fields by the anomaly in GOCI.

[Figure omitted. See PDF]

Figure 5d and e show the scatter plot analysis of the STK AOD products versus the AERONET AOD data with kriging variances (KVs). It is found that the STK AOD data with KV $\le$ 0.04 show a similar scattering pattern and accuracy to those of GOCI AOD. In contrast, some overestimated outliers from the STK AOD data in Fig. 5e (e.g., 1.0–2.0 in the $x$ axis and 2.0–4.0 in the $y$ axis) show different patterns than those from the GOCI AOD data. This may be explained by the relatively large KVs ($>$ 0.04) of such overestimated outliers. The KV generally increases when the observations near a certain prediction point are not available or when nearby observations have relatively large errors. Thus, when the GOCI observations are contaminated by optically thin clouds and they are not removed perfectly, this can increase the local variances due to their high cloud optical depth (COD). These factors can affect the quality of the STK AOD products. In this study, only the STK AOD products having small KVs (less than 0.04) were used for preparing the initial condition of each data processing step. Therefore, the initial PM concentrations did not changed where the STK AOD having large KVs (larger than 0.04; i.e. the right-bottom corner area in Fig. A2f). Collectively, it appears that the STK method is a reasonable tool for obtaining realistic AOD values at locations where the GOCI observations are not available.

Sensitivity of observation operators and control variables to AOD and PM${}_{10}$ predictions

To investigate the best combination of the observation operators and CVs, the AOD and PM${}_{10}$ hindcast runs and sensitivity analyses with the 12 different configurations (Table 3) were performed. For this, the hindcast AOD and PM${}_{10}$ from 13:00 to 19:00 LT (H $+$ 1 to H $+$ 6) on 10 selected episode days were compared with the ground-measured AOD and surface PM${}_{10}$. The observations from the six AERONET sites and nearest NAMIS PM${}_{10}$ stations within 10 km from the AERONET locations were selected for this comparison study (Fig. 2). The AOD values for the background CMAQ model simulations without the application of the STK method (noSTK) were also calculated with the GEOS-Chem observation operator.

Performance metrics for AOD and PM${}_{10}$ hindcasts on the 10 selected episodes at six AERONET sites and nearby NAMIS PM${}_{10}$ stations in SMA.

Configuration	AOD ($N^{\text{a}}$ $=$ 277)									PM${}_{10}$ ($N$ $=$ 340)
	IOA${}^{\text{b}}$	MFE${}^{\text{c}}$	MFB${}^{\text{d}}$	$R^{\text{e}}$	RMSE${}^{\text{f}}$	MB${}^{\text{g}}$	MNE${}^{\text{h}}$	MNB${}^{\text{i}}$		IOA	MFE	MFB	R	RMSE${}^{\text{j}}$	MB${}^{\text{j}}$	MNE	MNB
noSTK	0.48	113.2	$-$113.2	0.61	0.60	$-$0.53	70.0	$-$70.0		0.47	89.0	$-$88.5	0.54	55.15	$-$48.40	58.9	$-$58.4
A1	0.62	37.4	$-$22.1	0.46	0.36	$-$0.16	32.5	$-$13.7		0.60	35.4	$-$22.7	0.44	36.07	$-$15.80	31.1	$-$14.7
A2	0.60	39.8	$-$20.9	0.41	0.37	$-$0.15	35.8	$-$11.3		0.58	39.4	$-$34.7	0.50	37.13	$-$24.65	31.6	$-$25.8
A3	0.63	38.7	$-$22.5	0.46	0.36	$-$0.16	34.0	$-$13.5		0.64	33.0	$-$23.1	0.52	33.15	$-$17.07	28.4	$-$16.2
A4	0.63	37.4	$-$22.0	0.47	0.35	$-$0.16	32.6	$-$13.5		0.64	36.2	$-$28.3	0.53	34.58	$-$19.79	30.3	$-$20.4
B1	0.54	43.1	$-$27.1	0.33	0.40	$-$0.18	36.4	$-$16.2		0.53	41.1	$-$30.0	0.31	40.01	$-$20.90	33.9	$-$20.0
B2	0.51	44.7	$-$25.2	0.27	0.41	$-$0.17	39.5	$-$13.3		0.53	43.8	$-$39.2	0.37	40.94	$-$27.50	34.1	$-$28.5
B3	0.56	42.3	$-$25.8	0.35	0.39	$-$0.18	36.6	$-$15.2		0.56	38.0	$-$29.6	0.39	37.43	$-$21.65	31.2	$-$21.0
B4	0.55	41.9	$-$24.7	0.34	0.39	$-$0.17	36.3	$-$14.2		0.56	40.7	$-$33.9	0.42	38.30	$-$23.84	32.8	$-$24.4
C1	0.50	44.4	$-$37.5	0.28	0.43	$-$0.26	34.3	$-$26.3		0.55	35.8	$-$14.5	0.32	38.41	$-$9.82	33.4	$-$5.4
C2	0.47	45.7	$-$34.2	0.20	0.43	$-$0.24	36.7	$-$22.7		0.55	36.3	$-$26.0	0.37	36.86	$-$19.43	30.3	$-$17.6
C3	0.53	41.7	$-$30.5	0.34	0.40	$-$0.22	34.1	$-$20.5		0.60	32.9	$-$19.6	0.44	34.07	$-$14.92	29.1	$-$12.5
C4	0.53	41.7	$-$32.4	0.35	0.41	$-$0.23	32.4	$-$22.6		0.61	34.8	$-$21.8	0.44	34.78	$-$15.68	30.4	$-$14.0

${}^{\text{a}}$ The number of paired data, ${}^{\text{b}}$ index of agreement, ${}^{\text{c}}$ mean fractional error, ${}^{\text{d}}$ mean fractional bias, ${}^{\text{e}}$ Pearson product–moment correlation coefficient, ${}^{\text{f}}$ root mean square error, ${}^{\text{g}}$ mean bias, ${}^{\text{h}}$ mean normalized error, and ${}^{\text{i}}$ mean normalized bias. The units of all of metrics are dimensionless except ${}^{\text{j}}$ in microgram per cubic meter.

Figure 6 shows the soccer plot analysis of the 13 hindcast AODs (left panel) and PM${}_{10}$ (right panel) during the first 6 h of the short-term PM hindcasting on the 10 selected episode days. In the soccer plot, mean fractional bias (MFB) and mean fractional error (MFE) (described in Appendix B) are plotted on the $x$ and $y$ axes, respectively. Using this plot, the relative discrepancy can be presented by the distances from the origin of the plot and particular characteristics, such as systematic bias, can also be shown as a group of scatter points. Detailed statistical metric values are shown in Table 4. All the AODs and PM${}_{10}$ with the application of the STK method are much better than those from the noSTK simulation, with reduced errors and biases. Percentage decreases in MFE with the STK hindcasts were found to be 60–67 % for AOD and are 50–63 % for PM${}_{10}$. The MFB also decreased by 67–82 % for AOD and by 56–84 % for PM${}_{10}$. The noSTK case showed a strong negative bias (i.e., underprediction) and the 12 STK cases also showed less, yet still negative, biases. These negative biases are considered to be systematic, because of the negative bias of the GOCI AOD data (Fig. 6). Additionally, the negative biases are due to underestimation of CMAQ-simulated SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OA concentrations (Carlton et al., 2008, 2010; R. S. Park et al., 2011, 2014). This issue has been discussed in Sect. 2.5 and is investigated further in Sect. 3.4.

Soccer plot analysis for AOD (left panel) and PM${}_{10}$ (right panel) data from the first 6 h observations and the modeled data at six selected sites. BL (denoted by black diamond) represents the case of the bilinear interpolation method discussed in Sect. 3.2.

[Figure omitted. See PDF]

On the other hand, there are relatively small differences in errors and biases among the 12 STK cases (Fig. 6). Several differences among the 12 sensitivity cases were investigated further. First, the error and bias patterns for the AOD values were different from those for the PM${}_{10}$ predictions, being associated with the different observation operators. For example, the STK cases with the IMPROVE observation operator (cases C1, C2, C3, and C4) exhibited a relatively small bias for PM${}_{10}$ predictions, although they did not in the AOD predictions. This was likely caused by small wet MEE values of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ in the IMPROVE observation operator (represented by the green line in Fig. 4). In Eq. (1), the concentrations of converted aerosol species are inversely proportional to the MEEs of aerosol species. In the CV cases, the selections of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs (i.e., A3, B3, and C3) and SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and OAs (i.e., A4, B4, and C4) showed better performances for both the AOD and PM${}_{10}$ predictions.

Time series of hourly PM${}_{10}$ for the six sites over SMA for 9 April (a), 6 May (b), and 16 May (c) in 2012. Observed concentrations are shown as black circles and the model outputs as the colored lines with their own markers explained in the legend.

[Figure omitted. See PDF]

To show the degree of enhanced performances via the STK GOCI data, we also carried out some hindcast simulations, using the initial conditions prepared with single-frame GOCI data at 11:30 LT. The grids that did not have AOD observations were not filled out in this runs. In Fig. 6, the MFBs and MFEs of the bilinear interpolation method (denoted as BL) were $-$45.05 and 59.52 for AOD and $-$46.13 and 53.30 for PM${}_{10}$, respectively. It is shown in Fig. 6 that the use of the single-frame GOCI data without filling any gap cannot sufficiently improve the performance, compared with the cases of the STK simulations.

Figure 7 shows the performances of the short-term hindcast system with the 13 different configurations via comparisons between the hourly averaged PM${}_{10}$ observations and model PM${}_{10}$ predictions at the six NAMIS sites, on 9 April and 6 and 16 May 2012, respectively. Only 3-day and six-site results were selected and presented here, and more comprehensive performance evaluations are presented in Sect. 3.3. While noSTK failed to reproduce the high PM pollution, all the STK cases showed significant improvements in the surface PM${}_{10}$ predictions. However, there was a tendency in that the hourly peaks of PM${}_{10}$ were not well captured by the STK cases.

Performance metrics for PM${}_{10}$ hindcasting on the 10 selected episodes at 58 NAMIS PM${}_{10}$ stations in SMA. Abbreviations are the same as those in Table 3.

Configuration	H $+$ 1 to H $+$ 6 ($N$ $=$ 4823)									H $+$ 7 to H $+$ 12 ($N$ $=$ 4921)
	IOA	MFE	MFB	R	RMSE	MB	MNE	MNB		IOA	MFE	MFB	R	RMSE	MB	MNE	MNB
noSTK	0.45	99.6	$-$98.7	0.44	62.98	$-$54.59	63.9	$-$62.6		0.55	64.7	$-$36.9	0.30	56.76	$-$17.77	56.5	$-$12.9
A1	0.62	42.2	$-$30.9	0.47	40.64	$-$21.41	35.6	$-$19.7		0.62	43.9	1.5	0.37	49.17	5.27	51.2	19.1
A2	0.57	49.1	$-$43.4	0.48	43.81	$-$30.49	38.5	$-$30.4		0.60	45.1	$-$4.0	0.34	49.81	0.60	49.9	13.1
A3	0.64	40.5	$-$30.4	0.50	39.46	$-$21.83	34.2	$-$20.1		0.63	43.5	5.8	0.39	50.51	9.17	52.7	23.9
A4	0.63	44.6	$-$36.3	0.52	40.70	$-$24.99	36.3	$-$24.7		0.62	43.6	1.2	0.38	49.44	5.10	50.6	18.5
B1	0.54	48.8	$-$39.6	0.35	45.12	$-$27.64	38.8	$-$26.1		0.58	46.0	$-$3.6	0.31	49.18	0.71	50.8	14.1
B2	0.51	53.9	$-$48.9	0.36	47.76	$-$34.14	41.0	$-$33.9		0.59	46.3	$-$7.4	0.33	48.73	$-$2.37	49.1	9.3
B3	0.56	45.9	$-$37.9	0.40	43.36	$-$27.48	36.9	$-$25.8		0.61	44.6	0.7	0.35	48.81	4.07	51.0	17.9
A4	0.56	49.7	$-$43.0	0.43	44.46	$-$30.05	38.9	$-$29.6		0.60	45.2	$-$3.1	0.34	48.87	1.07	50.0	14.0
C1	0.60	40.4	$-$22.7	0.39	40.82	$-$15.98	35.9	$-$11.7		0.58	45.9	6.6	0.32	51.68	9.84	56.0	26.5
C2	0.56	43.7	$-$34.3	0.40	42.22	$-$25.43	35.6	$-$22.9		0.58	45.9	2.0	0.31	51.37	5.63	53.6	20.7
C3	0.63	39.0	$-$27.3	0.47	39.00	$-$20.13	33.3	$-$17.5		0.61	44.1	7.3	0.37	51.46	10.64	54.3	26.3
C4	0.63	41.2	$-$29.9	0.48	39.45	$-$21.30	34.5	$-$19.4		0.61	44.2	4.7	0.36	50.69	8.29	53.2	23.2

Consequently, it can be concluded that the combination of the GOCART observation operator and CVs of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs (represented by A3) leads to the best results in the current hindcast system (Table 4). The use of the GOCART observation operator and CVs of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and OAs (represented by A4) could also provide a comparable performance to A3. However, it appears that the differences among the 12 STK cases were relatively small.

Averaged PM${}_{10}$ of the noSTK case from H $+$ 1 to H $+$ 6 (a) and from H$+$7 to H$+$13 (b), and the averaged concentrations of case A3 at the same time series (c and d) for the selected 10 days. Averaged NAMIS PM${}_{10}$ observations are shown with colored circles.

[Figure omitted. See PDF]

Overall performance evaluation of PM${}_{10}$ hindcast over SMA

In this section, PM${}_{10}$ observations from the hindcast experiments were compared with the PM${}_{10}$ observations from 58 NAMIS sites to evaluate the overall performance of the current hindcast system in SMA. Table 5 provides the statistical metrics that were calculated separately from the first and the second 6 h hindcast results. The main characteristics of the statistical analysis in Table 5 are similar to those at the six sites discussed in the previous section. First, both errors and biases of PM${}_{10}$ distributions were significantly reduced after the application of the STK method. The MFEs and MFBs in the 12 h STK simulations decreased by $\sim$ 40 and $\sim$ 80 %, respectively.

A distinctive difference was also found in the model performances for the first and the second 6 h runs. During the first 6 h, all the hindcast results showed negative biases, with the MFB of $\sim$ $-$100 % for the noSTK cases and $\sim$ $-$40 % for the STK cases. The performances of the A3 and A4 cases are somewhat better than those of the other STK cases (Table 5). Collectively, the MFEs and MFBs of the STK cases are a factor of 2–4 smaller than those of the noSTK cases during the first 6 h.

Figure 8 shows a comparison between the noSTK case and the A3 case, in terms of the PM${}_{10}$ predictions, during the first and the next 6 h in SMA with the 6 h averaged NAMIS PM${}_{10}$ observations. As shown, the A3 case produced better PM${}_{10}$ predictions during the first and the next 6 h. In addition, the A4 case (not shown) also provided similar results to the A3 case, as discussed in Sect. 3.2. It can be confirmed again that the A3 and A4 cases are able to produce better PM${}_{10}$ predictions against the PM${}_{10}$ observations in SMA.

Hindcast performances from H $+$ 13 to H $+$ 24 were also evaluated with the ground-measured NAMIS PM${}_{10}$ data. In short, the differences between all the STK and noSTK cases became smaller than those during the first 12 h (an approximate difference of 10 % was found at H $+$ 24, i.e., 24 h after the hindcast actually began). Based on this, it appears that the effects of using the initial PM composition on the hindcast performances may effectively last during the first 12 h. After 12 h, the effects started to diminish. This is due to several facts: (i) the regions for applying the initial PM composition in this study were limited only within the GOCI domain (relatively small region); (ii) although the initial PM composition was used, its effects can be offset by uncertainties and errors in emissions as time progressed; and (iii) the large uncertainties associated with the formation of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs in the CTMs can also limit the effects of the initial PM composition. The latter two are the reasons that there is strong necessity for both emissions and CTMs to be improved continuously, even though the initial PM composition is applied in the short-term forecast activities.

Evaluation of hindcast performance with observed PM composition

In the previous section, PM${}_{10}$ mass concentrations were simply predicted by the short-term hindcast system with 12 different combinations of observation operators and CVs. Although the purpose of this study is to develop a better PM forecast system for accurately predicting PM${}_{10}$ mass concentrations, it is still necessary to more carefully scrutinize the changes in the PM composition in accordance with the different selections of the CVs.

During the DRAGON-Asia campaign, the PM${}_{2.5}$ composition was measured for SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ with 30 min intervals and for SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, OC and BC with 24 h intervals using the PILS-IC instrument (semicontinuous measurements) and low-air-volume sampler with a Teflon filter (offline measurements), respectively, in Yongin City near SMA (Fig. 2). Thus, in this section, the selection of the CVs is further discussed with the observed PM${}_{2.5}$ composition.

Figure 9 shows the comparison between 1 h averaged SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ concentrations measured via the PILS-IC instrument and model-predicted concentrations during the selected days at the Yongin observation site. Only the STK cases with the GOCART observation operator (i.e., A1, A2, A3, and A4) were selected here. The STK cases showed significant changes in the PM composition with the selection of CVs. For example, the A2 and A3 cases tended to overestimate the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations but underestimated the NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ concentrations, whereas the A1 and A4 cases tended to capture relatively well the trend of the concentrations of the three particulate species. This phenomenon was driven by intraparticulate thermodynamics. That is, if larger amounts of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ are allocated into particles (like the cases of A2 and A3), then NO${}_{\mathrm{3}}^{-}$ tends to be evaporated, because SO${}_{\mathrm{4}}^{\mathrm{2}-}$ is more strongly associated with NH${}_{\mathrm{4}}^{+}$ (Song and Carmichael, 1999). As shown in Fig. 9a and b, when the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations increase (as in case A2), the NO${}_{\mathrm{3}}^{-}$ concentrations decrease accordingly, because NO${}_{\mathrm{3}}^{-}$ is evaporated out of the particulate phase as a form of HNO${}_{\mathrm{3}}$ (Song and Carmichael, 1999, 2001). Collectively, the “best” results were produced from the case A4, as shown in Fig. 9a–c.

Time-series comparison of 1 h averaged (a) SO${}_{\mathrm{4}}^{\mathrm{2}-}$, (b) NO${}_{\mathrm{3}}^{-}$, and (c) NH${}_{\mathrm{4}}^{+}$ concentrations measured with the PILS-IC instrument and model-predicted concentrations. In (d), 24 h averaged aerosol concentrations in PM${}_{2.5}$ from observations (PILS-IC instrument and low-air-volume sampler with Teflon filter) are compared with hindcast concentrations at the Yongin City site for 10 selected episodes.

[Figure omitted. See PDF]

Time series of hourly PM${}_{10}$ at six sites in SMA for 8 March (a), 10 March (b), and 11 March (c) in 2012. Observed concentrations are denoted as black circles and the modeled concentrations as colored lines.

[Figure omitted. See PDF]

The 24 h averaged PM${}_{2.5}$ compositions measured from the PILS-IC instrument and the low-air-volume sampler with a Teflon filter during the campaign period are also compared in Fig. 9d. Again, the observations of the SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ concentrations were obtained from both the PILS-IC instrument and the low-volume sampler, whereas the concentrations of OAs ($\cong$ [OC] $\times$ 1.5) and EC were only measured via the low-air-volume sampler. As shown in Fig. 9d, the SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ concentrations from both samplers showed good agreements (see circles and crosses in Fig. 9d). The A4 case (the red bars in Fig. 9d) again showed the best results in the comparison between the observed and predicted particulate composition, particularly in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs. In the previous discussion (see Sect. 3.2 and 3.3), the A3 and A4 cases showed the best performances for predicting PM${}_{10}$ mass concentrations over SMA. This is somewhat consistent with our analysis in this section. However, in the case of A3, it can capture the PM mass behaviors (Sect. 3.3) but does not capture the changes in the PM composition well (this section). Based on this, it is concluded that the A4 case would be the best configuration for accurately predicting the PM composition as well as the PM mass. However, this PM composition analysis was conducted with observations from only one site (Yongin City) in this study. Thus, to reach a firmer conclusion, more intensive analyses with observations from multiple sites are required in future.

Evaluation of short-term hindcast performances

To further evaluate the performance of the short-term hindcast runs, 48 h hindcast simulations with the configuration of A4 were carried out from 7 to 19 March. The observations from the six AERONET sites and the nearest NAMIS stations were analyzed in this study.

The time series of the first and the second 24 h averaged PM${}_{10}$ at the six sites on 8, 10, and 11 March 2012 are presented in Fig. 10. Again, reduced errors and biases were shown in the A4 STK simulations, compared with the noSTK simulation for polluted episodes (panels a and b in Fig. 10) and for the less-polluted episode (panel c in Fig. 10). Percent decreases with MFEs of the first 24 h A4 STK hindcast were $\sim$ 40 % for AOD and $\sim$ 10 % for PM${}_{10}$, and those with MFBs were $\sim$ 40 % for AOD and $\sim$ 100 % for PM${}_{10}$. In addition, slight improvements in the horizontal distributions of AOD and PM${}_{10}$ were also found. This was indicated by the increases of correlation coefficients (refer to Table S1 in the Supplement). The second 24 h STK hindcasts also reduced the errors and biases for AOD and PM${}_{10}$, although the improvements in the spatial distributions were not shown clearly. More detailed statistical metrics are presented in the Supplement (Table S1).

Summary and conclusions

For the purpose of improving the performance of short-term PM forecast in South Korea, an integrated air quality modeling system was developed with the application of the STK method using the geostationary satellite-derived AOD data over Northeast Asia. The errors and biases of the STK AOD showed relatively good agreement, compared with the AERONET observations. With the combinations of the STK method along with various observation operators and CVs, the errors and biases of AOD and PM${}_{10}$ predictions can be reduced significantly. It was shown that the selection of the observation operators greatly influence the performances of the STK hindcast systems. On the other hand, the choice of CVs tends to affect PM composition. The combination of the GOCART observation operator and the selection of CVs of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OAs (case A3) was found to be the best one for the PM${}_{10}$ mass prediction. All the hindcast runs with the application of the STK method, however, generally showed negative biases (i.e., underpredictions). This was primarily due to the underestimation of the GOCI AOD.

Reducing errors and biases in the current system is important for further development of the PM forecast system. One of the potential methods for reducing the errors and biases is to introduce the MODIS AOD data into the STK stage, together with the GOCI data. It is expected that doing this will further reduce the systematic biases, due to the relatively smaller biases of MODIS AOD (as shown Fig. 3). In addition, the combination of the GOCART observation operator and the selection of CVs of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and OAs (Case A4) was found to give the “best” results for the prediction of particulate composition at one observation site. However, more intensive measurements of the PM composition are needed for reaching a more solid conclusion.

The STK AODs used in the current study are expected to be used in other data assimilation methods. For example, in the 3DVAR method, the observation error covariance matrix, which presents the degree of errors of the observations, has been usually assumed by linear equations or a single constant value (Liu et al., 2011; Schwartz et al., 2012; Shi et al., 2011). However, as discussed with KVs in Sect. 3.1, the error covariance of the AOD observations can be improved and the use of the improved observation error covariance matrix can help to prepare more accurate AOD fields, for example, via a 3DVAR method. This study is now underway.

In future, planned GEO satellite sensors will give other opportunities to use semicontinuous AOD observations at high spatial and temporal resolutions. Upcoming GEO satellite sensors scheduled for launch between 2018 and 2020 include NASA's Tropospheric Emissions: Monitoring of Pollution (TEMPO) over North America, ESA's Sentinel-4 over Europe, and Korea Aerospace Research Institute (KARI)'s Geostationary Environment Monitoring Spectrometer (GEMS) over Asia. In the case of the GEMS instrument, it is being designed to provide backscattered UV/Vis radiances between 300 and 500 nm with a spatial resolution of 5 km $\times$ 5 km over a large part of Asia. Using advanced observations from the GEMS sensor, it is anticipated that the system developed here will be able to make significant contributions to further improvements in the performances of the PM forecasting system in Asia. This improved PM predictions and modeling framework can also be a core part of the air quality forecasting system, a more comprehensive health impact assessments, and radiative forcing estimation over (East) Asia in future.

Spatiotemporal-kriging method

The STK methods assume that measured variables in space and time ($\mathit{\tau }(s,t))$ can be regarded as a random function, consisting of a trend component ($m)$ and residual component ($\mathit{ϵ})$ of which the mean is zero: $$\mathit{\tau }(s,t)=m(s,t)+\mathit{\epsilon }(s,t).$$ The unobserved value ${\mathit{\tau }}^{\ast }(s,t)$ can be averaged with weight using measured values from the surroundings: $${\mathit{\tau }}^{\ast }(s,t)=\sum _{i=\mathrm{1}}^n\mathit{\tau }(s_i,t_i)w_i(s,t),$$ where $n$ is the number of observations in the local neighborhood and $w_i(s,t)$ is the kriging weight assigned to $\mathit{\tau }(s_i,t_i)$. The kriging weight is determined by a theoretical semivariogram.

In case of spatial kriging ($\mathit{\tau }\left(s\right))$, the semivariogram ($\mathit{\gamma })$ is the best fit to the semivariance (${\mathit{\gamma }}^{\ast })$ as a function of spatial lag ($h)$. Assuming the trend component $m\left(s\right)$ in $\mathit{\tau }\left(s\right)$ is constant over the local domain (i.e., the ordinary kriging method), the semivariance is defined as $$\begin{array}{cc}{\mathit{\gamma }}^{\ast }\left(h\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}N\left(h\right)}}\sum _{i=\mathrm{1}}^{N\left(h\right)}\left[\mathit{\tau }\right(s_i)-\mathit{\tau }(s_i+h){]}^{\mathrm{2}}\\ & & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}N\left(h\right)}}\sum _{i=\mathrm{1}}^{N\left(h\right)}\left[\mathit{\epsilon }\right(s_i)-\mathit{\epsilon }(s_i+h){]}^{\mathrm{2}},\end{array}$$ where $N\left(h\right)$ is the number of paired observations at a spatial distance of $h$, and ${\mathit{\tau }}_i(s_i+h)$ is the $i$th observation (in this study, AOD) separated by $h$ from the observation located at $s_i$. The semivariogram is then depicted by a theoretical model which is the best-fitting curve to the semivariance by minimizing the least squares error. For example, a spherical semivariogram ($\mathit{\gamma })$, which is commonly used in the theoretical models of the atmospheric studies, is estimated by finding three optimal parameters: (i) nugget  ($c_n)$, (ii) range ($a)$, and (iii) partial sill (${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}})$: $$\begin{array}{cc}& \mathit{\gamma }\left(h\right)=c_n+{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}\left[{\displaystyle \frac{\mathrm{3}h}{\mathrm{2}a}},-,{\displaystyle \frac{h^{\mathrm{3}}}{\mathrm{2}{a}^{\mathrm{3}}}}\right](\text{for}h\le a),\\ & & \mathit{\gamma }\left(h\right)=c_n+{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}(\text{for}h>a).\end{array}$$ The range parameter indicates the maximum lag in which the variation of semivariogram is meaningful (Cressie, 1992).

Two daily three-dimensional semi-variograms from the GOCI AOD data on 8 April 2012: (a) fitted by the spherical model (Eq. A4), and (b) estimated by the sampled observations.

[Figure omitted. See PDF]

To combine the spatial and temporal data for preparing the spatiotemporal semivariograms, the temporal information can be converted into the spatial information (Gräler et al., 2012). First, the spatial and temporal semivariograms are estimated independently using the spherical model from the daily GOCI AOD data. Second, the ratio of the spatial range parameter ($a_{\text{s}})$ of the spatial semivariogram to the temporal range parameter ($a_{\text{t}})$ of the temporal semivariogram (i.e., spatiotemporal-scale factor, km h${}^{-\mathrm{1}})$ is used to convert the unit of temporal lag into the unit of spatial distance. Consequentially, the 3-D spatiotemporal AOD data are converted into the 2-D spatial AOD fields. After that, the spatiotemporal semivariogram is provided to predict the AOD fields with 15 km $\times$ 15 km spatial resolution from 10:00 to 16:00 LT over the GOCI domain. For the STK method, the “gstat” (Pebesma, 2004) and the “spacetime” (Pebesma, 2012) software packages in the R environment for statistical computing were used (R Development Core Team, 2011). Figure A1 presents an example of the 3-D semivariograms from the fitted model (left) and sample from the GOCI data on 8 April. The mean nugget ($c_{\text{n}})$, range ($a)$, and partial sill (${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}})$ of the spatiotemporal model semivariogram were 0.025, 583 km, and 0.227, respectively, during the entire DRAGON-Asia campaign. The average spatiotemporal-scale factor of $\sim$ 34 km h${}^{-\mathrm{1}}$ was calculated indicating that the AODs observed before or after 1 h at certain locations show a similar correlation pattern to those measured simultaneously at $\sim$ 34 km apart in the STK model. Figure A2 shows an example of spatial distributions of GOCI AOD from 10:30 to 13:30 LT and STK AOD at 12:00 LT with a criteria of KVs of less than 0.04.

Spatial distributions of GOCI AOD from 10:30 to 13:30 LT (a–d) and STK AOD at 12:00 LT (e) on 7 April 2012. The STK AOD at 12:00 LT with a criteria of KVs of less than 0.04 is also shown in (f).

[Figure omitted. See PDF]

Statistical metrics

In this study, eight statistical metrics were used for validating the hindcast results (Chai and Draxler, 2014; Savage et al., 2013; Willmott, 1981; Willmott et al., 2009; Willmott and Matsuura, 2005). $$\begin{array}{cc}& \text{Index of agreement (IOA)}=\\ & & \mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left(O_i,-,M_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left(|,O_i,-,\stackrel{\mathrm{‾}}{O_i},|,+,\left|,M_i,\stackrel{\mathrm{‾}}{M_i},\right|\right)}^{\mathrm{2}}}}& \text{Mean fractional error (MFE)}=\\ & & {\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\displaystyle \frac{|M_i-O_i|}{\left(\frac{M_i+O_i}{\mathrm{2}}\right)}}\times 100& \text{Mean fractional bias (MFB)}=\\ & & {\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\displaystyle \frac{\left(M_i,-,O_i\right)}{\left(\frac{M_i+O_i}{\mathrm{2}}\right)}}\times 100& \text{Regression coefficient}\left(R\right)=\\ & & {\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(O_i,-,\stackrel{\mathrm{‾}}{O_i}\right)(M_i-\stackrel{\mathrm{‾}}{M_i})}{\sqrt{\sum _{i=\mathrm{1}}^N(O_i-\stackrel{\mathrm{‾}}{O_i}{)}^{\mathrm{2}}}\sqrt{\sum _{i=\mathrm{1}}^N(M_i-\stackrel{\mathrm{‾}}{M_i}{)}^{\mathrm{2}}}}}\end{array}$$

$$\begin{array}{cc}& \text{Root mean square error (RMSE)}=\\ & & \sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(M_i,-,O_i\right)}^{\mathrm{2}}}& \text{Mean normalized error (MNE)}=\\ & & {\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left({\displaystyle \frac{|M_i-O_i|}{O_i}}\right)\times 100& & \text{Mean bias (MB)}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left(M_i,-,O_i\right)& \text{Mean Normalized bias (MNB)}=\\ & & {\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left({\displaystyle \frac{M_i-O_i}{O_i}}\right)\times 100\end{array}$$ In Eqs. (B1)–(B8), $N$ is the number of data and $M_i$ and $O_i$ are the model value and observation, respectively. The variables with an overbar are the arithmetic mean of the data.

Code availability

WRF and CMAQ source codes and R and NCL computer languages are available to the public. The source codes and computer languages may be downloaded by following instructions found at http://www2.mmm.ucar.edu/wrf/users/downloads.html for WRF, https://www.cmascenter.org/cmaq for CMAQ, http://cran.r-project.org for R, and https://www.ncl.ucar.edu/Download for NCL.

The STK module code used in this study was based on the instruction of Pebesma (2012) available at http://www.jstatsoft.org/v51/i07 and can be obtained by contacting S. Lee ([email protected]).

The Supplement related to this article is available online at doi:10.5194/gmd-9-17-2016-supplement.

Acknowledgements

This research was supported by the GEMS program of the Ministry of Environment, South Korea, as part of the Eco Innovation Program of KEITI (2012000160004). This work was also funded by the Korea Meteorological Administration Research and Development Program under Grant KMIPA 2015–5010 and was partly supported by the National Institute of Environmental Research (NIER). We thank all PI investigators and their staff for establishing and maintaining the AERONET sites of the DRAGON NE Asia 2012 campaign used in this study. We also thank the MODIS science team for providing valuable data for this research. NCL (2014) was used to draw the figures. The third author was supported by the research and development project on the development of global numerical weather prediction systems of the Korea Institute of Atmospheric Prediction Systems (KIAPS) funded by the Korea Meteorological Administration (KMA). Edited by: O. Boucher