Development of the Global to Mesoscale Air Quality Forecast and Analysis System (GMAF) and Its Application to PM 2.5 Forecast in Korea

: This paper presents the development of the global to mesoscale air quality forecast and analysis system (GMAF) and its application to particulate matter under 2.5 µ m (PM 2.5 ) forecast in Korea. The GMAF combined a mesoscale model with a global data assimilation system by the grid nudging based four-dimensional data assimilation (FDDA). The grid nudging based FDDA developed for weather forecast and analysis was extended to air quality forecast and analysis for the ﬁrst time as an alternative to data assimilation of surface monitoring data. The below cloud scavenging module and the secondary organic formation module of the community multiscale air quality model (CMAQ) were modiﬁed and subsequently veriﬁed by comparing with the PM speciation observation from the PM supersite. The observation data collected from the criteria air pollutant monitoring networks in Korea were used to evaluate forecast performance of GMAF for the year of 2016. The GMAF showed good performance in forecasting the daily mean PM 2.5 concentrations at Seoul; the correlation coefﬁcient between the observed and forecasted PM 2.5 concentrations was 0.78; the normalized mean error was 25%; the probability of detection for the events exceeding the national PM 2.5 standard was 0.81 whereas the false alarm rate was only 0.38. Both the hybrid bias correction technique and the Kalman ﬁlter bias adjustment technique were implemented into the GMAF as postprocessors. For the continuous and the categorical performance metrics examined, the Kalman ﬁlter bias adjustment technique performed better than the hybrid bias correction technique.


Introduction
Chemical transport models (CTMs) were first developed in the 1980s to investigate the transport [1][2][3] and the transformation [4,5] of air pollutants. The CTMs have been further developed by sophisticating atmospheric processes such as gas and aqueous chemical reactions [6], dry and wet deposition [7,8], and aerosol formation [9] to address mesoscale air pollution problems including acid deposition, tropospheric ozone (O 3 ), fine particulate matter (PM) pollution. The CTMs addressing air quality problems are named as air quality models. Most of mesoscale air quality models do not simulate atmospheric dynamics internally, requiring a mesoscale atmospheric model to provide the meteorological variables.
On the other hand, global atmospheric chemistry models have progressed rapidly to forecast and analyze the air quality despite of its coarse grid resolution and simplification of physical and chemical processes [10]. One of advantages of using global atmospheric chemistry models is that it is easy to incorporate polar orbiting satellite observations into data assimilation. The spatial and temporal resolution of polar orbiting satellite is too coarse to be used for mesoscale data assimilation. The remedy for this problem is to constrain the mesoscale model outputs by the global model outputs by using a grid nudging based fourdimensional data assimilation (FDDA) technique [11]. The grid nudging based FDDA is As shown in Figure 2, the GMAF imported CMAQ [23] for air quality forecast and analysis. The CMAQ version 5.3.1 was used with the carbon bond version 6 mechanism [24] as a gas-phase chemical mechanism, which includes 77 chemical species and 218 chemical reactions. The aerosol dynamics were modeled by AERO7 module, which includes 16 inorganic aerosol species, 33 organic aerosol species, as described in the github of CMAQ website [25]. The CMAQ divides sea salt into individual inorganic aerosol species and includes them except the coarse mode cations. Similarly, the CMAQ divides dust into individual aerosol species except the unknown composition fraction. The CMAQ treats the coarse mode cations of sea salt and the unknown composition fraction of dust as lumped aerosol species. The below-cloud scavenging module and the secondary organic aerosol formation mechanism were modified, and the evaporation loss of nitrate was added as detailed in Section 2.4. The WRF-ARW model version 3.6 [26] was used to provide the values of meteorological variables to the CMAQ.
Global assimilated datasets were used for generation of initial and boundary conditions of the outer domain and for providing gridded data to the grid nudging based FDDA. The WRF used the 0.25° global data assimilation system (GDAS) dataset and the CMAQ used the 0.4° CAMS dataset. The CAMS dataset was produced by the ECMWF's IFS using four-dimensional variational data assimilation of the polar orbiting satellite data. A hybrid, sigma-pressure vertical coordinate was used with 137 levels and the top level was extended to 0.01 hPa. The CAMS dataset provided the concentration of 14 gaseous species and seven aerosol species in 3 h intervals [27]. The gaseous species were nitric oxide (NO), nitrogen dioxide (NO2), sulfur dioxide (SO2), formaldehyde (HCHO), ozone (O3), peroxyacetyl nitrates (PAN), hydrogen peroxide (H2O2), nitric acid (HNO3), hydroxyl radical (OH), carbon monoxide (CO), methane, ethane, propane, and isoprene. The aerosol species were dust, sea salt, hydrophobic black carbon, hydrophilic black carbon, hydrophobic organic matter, hydrophilic organic matter, and sulfate. The dust and the sea As shown in Figure 2, the GMAF imported CMAQ [23] for air quality forecast and analysis. The CMAQ version 5.3.1 was used with the carbon bond version 6 mechanism [24] as a gas-phase chemical mechanism, which includes 77 chemical species and 218 chemical reactions. The aerosol dynamics were modeled by AERO7 module, which includes 16 inorganic aerosol species, 33 organic aerosol species, as described in the github of CMAQ website [25]. The CMAQ divides sea salt into individual inorganic aerosol species and includes them except the coarse mode cations. Similarly, the CMAQ divides dust into individual aerosol species except the unknown composition fraction. The CMAQ treats the coarse mode cations of sea salt and the unknown composition fraction of dust as lumped aerosol species. The below-cloud scavenging module and the secondary organic aerosol formation mechanism were modified, and the evaporation loss of nitrate was added as detailed in Section 2.4. The WRF-ARW model version 3.6 [26] was used to provide the values of meteorological variables to the CMAQ.
Global assimilated datasets were used for generation of initial and boundary conditions of the outer domain and for providing gridded data to the grid nudging based FDDA. The WRF used the 0.25 • global data assimilation system (GDAS) dataset and the CMAQ used the 0.4 • CAMS dataset. The CAMS dataset was produced by the ECMWF's IFS using four-dimensional variational data assimilation of the polar orbiting satellite data. A hybrid, sigma-pressure vertical coordinate was used with 137 levels and the top level was extended to 0.01 hPa. The CAMS dataset provided the concentration of 14 gaseous species and seven aerosol species in 3 h intervals [27]. The gaseous species were nitric oxide (NO), nitrogen dioxide (NO 2 ), sulfur dioxide (SO 2 ), formaldehyde (HCHO), ozone (O 3 ), peroxyacetyl nitrates (PAN), hydrogen peroxide (H 2 O 2 ), nitric acid (HNO 3 ), hydroxyl radical (OH), carbon monoxide (CO), methane, ethane, propane, and isoprene. The aerosol species were dust, sea salt, hydrophobic black carbon, hydrophilic black carbon, hydrophobic organic matter, hydrophilic organic matter, and sulfate. The dust and the sea salt were distributed into the three size bins: 0.030-0.55, 0.55-0.9 and 0.9-20 µm for the dust and 0.03-0.5, 0.5-5 and 5-20 µm for the sea salt.
The clean air policy support system (CAPSS) emission inventory [28] was used to generate the gridded anthropogenic emission of South Korea. The Korea-United States air quality study (KORUS-AQ) anthropogenic emission inventory [29] was used for the rest of the region. The biogenic emissions were generated by MEGAN version 2.1 [30]. The CAPSS

Forecasting Period and Area
The first national air quality standard for PM2.5 was promulgated in Korea in 2015 to begin installation of a PM2.5 monitoring instrument at the criteria air pollutant monitoring station. Most of the criteria air pollutant monitoring stations in the urban area were

Forecasting Period and Area
The first national air quality standard for PM 2.5 was promulgated in Korea in 2015 to begin installation of a PM 2.5 monitoring instrument at the criteria air pollutant monitoring station. Most of the criteria air pollutant monitoring stations in the urban area were equipped with a PM 2.5 monitoring instrument by 2016. Furthermore, KORUS-AQ field campaign was conducted in the spring of 2016 to provide the measurements from aircraft, ground sites and ships. The KORUS-AQ emissions were also developed to support the campaign [32]. Taking into account availability of monitoring data and emission inventory, this study chose the period from December of 2015 to November of 2016 for verification equipped with a PM2.5 monitoring instrument by 2016. Furthermore, KORUS-AQ field campaign was conducted in the spring of 2016 to provide the measurements from aircraft, ground sites and ships. The KORUS-AQ emissions were also developed to support the campaign [32]. Taking into account availability of monitoring data and emission inventory, this study chose the period from December of 2015 to November of 2016 for verification and evaluation of the GMAF. Four major cities in Korea were selected as target areas. Those are Seoul, Busan, Gwangju, and Daegu as displayed in Figure 3. The PM2.5 pollution in Seoul and Busan expressed in Figure 4, the first and second largest city in Korea, was the worst in 2016 in the period from 2015 to 2019. As shown in Table 1, the annual mean PM2.5 concentration was the highest and the number of days exceeding the PM2.5 standard of 35 μg/m 3 was the largest in 2016. Among meteorological parameters, precipitation, wind velocity and wind direction affect the year-to-year variations of the PM2.5 concentrations the most. The precipitation from June to September of 2016 was 52% less than the 30-year average in Seoul. This low precipitation reduced the wet scavenging to increase the PM2.5 concentrations. The PM2.5 concentrations from June to September of 2016 were higher than the other year averages of the same period by 33%. The number of days of exceeding the PM2.5 standard of 35 μg/m 3 from June to September of 2016 was larger than the other year averages of the same period by 258% in Seoul as shown in Figure 4. The PM 2.5 pollution in Seoul and Busan expressed in Figure 4, the first and second largest city in Korea, was the worst in 2016 in the period from 2015 to 2019. As shown in Table 1, the annual mean PM 2.5 concentration was the highest and the number of days exceeding the PM 2.5 standard of 35 µg/m 3 was the largest in 2016. Among meteorological parameters, precipitation, wind velocity and wind direction affect the year-to-year variations of the PM 2.5 concentrations the most. The precipitation from June to September of 2016 was 52% less than the 30-year average in Seoul. This low precipitation reduced the wet scavenging to increase the PM 2.5 concentrations. The PM 2.5 concentrations from June to September of 2016 were higher than the other year averages of the same period by 33%. The number of days of exceeding the PM 2.5 standard of 35 µg/m 3 from June to September of 2016 was larger than the other year averages of the same period by 258% in Seoul as shown in Figure 4.  Busan is located about 320 km southeast of Seoul. The observed precipitation from May to August in 2016 was comparable to the 30-year average in Busan, and therefore the precipitation is unlikely to contribute to the PM2.5 increase. Instead, the stagnant wind conditions that occurred in November 2016 increased the monthly mean concentrations and the number of days of exceedances as shown in Figure 5.   Busan is located about 320 km southeast of Seoul. The observed precipitation from May to August in 2016 was comparable to the 30-year average in Busan, and therefore the precipitation is unlikely to contribute to the PM 2.5 increase. Instead, the stagnant wind conditions that occurred in November 2016 increased the monthly mean concentrations and the number of days of exceedances as shown in Figure 5.

Grid Nudging Based FDDA Method
Data assimilation combines model with observation to improve performance of the forecast and the reanalysis [33]. The grid nudging based FDDA adds a relaxation term that is supposed to force the observations to the model [34] as where Y is the dependent variable, F is a discretized form of the governing equation, x is the independent spatial variable, and W is the nudging coefficient. The subscripts 'obs' and 'model' denote the observations and the model, respectively. Global assimilated datasets may be used instead of observations to combine the mesoscale model with the global model. The dependent variables, Y, are defined as = for the CMAQ, where is the meteorological variable of the WRF, is the chemical species concentration, pt is the air pressure at the top of the domain, ps is the air pressure at the surface and is the vertical Jacobian of the terrain-influenced coordinates. The variables nudged for the WRF were U, V winds, temperature and water vapor mixing ratio. The water vapor mixing ratio in the PBL and the surface temperature were not nudged. The nudging coefficients used for the WRF in this study were set as recommended by Baker et al. (2009) [35] and listed for completeness in Table 2.
The nudged variables of the CMAQ were selected to be consistent with CAMS dataset. The primary pollutants were selected to minimize the emission uncertainties in the outer domain. The O3 was added to the FDDA nudged variable list because its satellite observation data which the CAMS based on were reasonably accurate. In addition, the grid based FDDA on O3 enables incorporation of stratospheric intrusion into the mesoscale model domain. The grid spacing of the CAMS dataset is 0.4° (about 40 km in

Grid Nudging Based FDDA Method
Data assimilation combines model with observation to improve performance of the forecast and the reanalysis [33]. The grid nudging based FDDA adds a relaxation term that is supposed to force the observations to the model [34] as where Y is the dependent variable, F is a discretized form of the governing equation, x is the independent spatial variable, and W is the nudging coefficient. The subscripts 'obs' and 'model' denote the observations and the model, respectively. Global assimilated datasets may be used instead of observations to combine the mesoscale model with the global model. The dependent variables, Y, are defined as where α is the meteorological variable of the WRF, ϕ is the chemical species concentration, p t is the air pressure at the top of the domain, p s is the air pressure at the surface and J s is the vertical Jacobian of the terrain-influenced coordinates. The variables nudged for the WRF were U, V winds, temperature and water vapor mixing ratio. The water vapor mixing ratio in the PBL and the surface temperature were not nudged. The nudging coefficients used for the WRF in this study were set as recommended by Baker et al. (2009) [35] and listed for completeness in Table 2. Water vapor mixing ratio The nudged variables of the CMAQ were selected to be consistent with CAMS dataset. The primary pollutants were selected to minimize the emission uncertainties in the outer domain. The O 3 was added to the FDDA nudged variable list because its satellite observation data which the CAMS based on were reasonably accurate. In addition, the grid based FDDA on O 3 enables incorporation of stratospheric intrusion into the mesoscale model domain. The grid spacing of the CAMS dataset is 0.4 • (about 40 km in mid-latitudes) and the grid spacing of outer domain and inner nest of the GMAF are 27 and 9 km, respectively. Interpolation of the CAMS dataset to the inner nest may severely underestimate the peak values due to a large difference between the grid resolutions. Moreover, the emission uncertainties in the inner nest, Korea emission, were reasonably small. Therefore, the grid nudging based on FDDA was applied only to the outer domain. Table 2 presents the list of nudged chemical species with the corresponding nudging coefficients.

Modification of the CMAQ
The CMAQ was modified to enhance the model accuracy based on the recent findings deduced from model verification as described below.

Below-Cloud Scavenging
The CMAQ assumes that accumulation mode and coarse mode particles are completely absorbed into cloud and rainwater. This assumption is valid only for a cloud environment where water droplets are readily formed around particles due to high water vapor contents but not valid for a rain environment. It may lead to overestimation of below-cloud scavenging and subsequently to underestimation of particulate matter concentrations in the below-cloud region.
In order to remove the complete absorption assumption in a rain environment, the CMAQ was modified to calculate the below-cloud scavenging rate by where c(t) is the particulate species concentration at time t and the Λ is the scavenging coefficient [36]. The below-cloud scavenging coefficient was calculated by semiempirical formula derived by Slinn (1983) [37] with the empirically estimated mean raindrop diameter [38]. The calculated below-cloud scavenging coefficients ranged 10 −6~1 0 −2 s, consistent with the field measurements [39]. If the particles are completely absorbed into rainwater as the CMAQ assumes, then the below-cloud scavenging coefficient should be equal to 1, which is several orders of magnitude larger than the one estimated here. More importantly, it remains constant whereas it should increase with the precipitation intensity. As a result, the CMAQ seriously overestimates the below-cloud scavenging rate and subsequently underestimates PM 2.5 concentration especially when the precipitation intensity is low.
The WRF-CMAQ simulation was performed from 21 March to 27 April in 2016 to test the modified below-cloud scavenging module. Figure 6 compares the modeled and observed precipitations at Seoul. The modeled and observed surface PM 2.5 concentrations were also compared. The modeled precipitation agrees very well with the observation except 21 April. The dates of the modeled precipitation between 0.1 mm/day and 1.0 mm/day are 29, 30 March, and 3, 6,8,9,10,14,19,22 April. As also shown in Figure 6, the CMAQ without modification considerably underestimated the PM 2.5 concentrations on these dates of little precipitation. As noted earlier, this is because the CMAQ severely overestimated the below-cloud scavenging rates by assuming complete absorption of particles into rainwater. The CMAQ with the proposed modification reduced overestimation of the below-cloud scavenging rates, resulting that the modeled PM 2.5 concentration agreed well with the observation.
Atmosphere 2021, 12, x FOR PEER REVIEW 9 of 23 of the below-cloud scavenging rates, resulting that the modeled PM2.5 concentration agreed well with the observation.

Secondary Organic Aerosol Formation
In the present study, the AERO7 module including 12 primary organic aerosols (POA) and 21 secondary organic aerosols (SOAs) were used to model SOA formation. Compared to the previous versions of aerosol module, the AERO7 module expanded the SOA formation mechanism by adding pcSOA, organic nitrates from monoterpene oxidation and SOA from glyoxal and methylglyoxal. The pcSOA, standing for the potential SOA from combustion emissions, is a surrogate species designed to compensate missing and underestimated mass from combustion.
The pcSOA production rate of the AERO7 module was optimized for the southern and northern California sites and was known to overestimate the SOA in the continental US application [40]. Figure 7 displays the model calculated SOA compositions by season in Seoul. Here, the SOA was divided into four groups, pcSOA, SOA from biogenic VOCs (SOA_BVOC), SOA from anthropogenic VOCs (SOA_AVOC), and SOA from glyoxal and methylglyoxal (SOA_gly) for convenience. The pcSOA was the most dominant SOA in the winter due to low temperatures whereas the SOA_BVOC was the most dominant SOA in the summer due to high emission fluxes of biogenic VOCs and fast photochemical production. These seasonal characteristics of SOA compositions appeared similarly in most urban areas in Korea.

Secondary Organic Aerosol Formation
In the present study, the AERO7 module including 12 primary organic aerosols (POA) and 21 secondary organic aerosols (SOAs) were used to model SOA formation. Compared to the previous versions of aerosol module, the AERO7 module expanded the SOA formation mechanism by adding pcSOA, organic nitrates from monoterpene oxidation and SOA from glyoxal and methylglyoxal. The pcSOA, standing for the potential SOA from combustion emissions, is a surrogate species designed to compensate missing and underestimated mass from combustion.
The pcSOA production rate of the AERO7 module was optimized for the southern and northern California sites and was known to overestimate the SOA in the continental US application [40]. Figure 7 displays the model calculated SOA compositions by season in Seoul. Here, the SOA was divided into four groups, pcSOA, SOA from biogenic VOCs (SOA_BVOC), SOA from anthropogenic VOCs (SOA_AVOC), and SOA from glyoxal and methylglyoxal (SOA_gly) for convenience. The pcSOA was the most dominant SOA in the winter due to low temperatures whereas the SOA_BVOC was the most dominant SOA in the summer due to high emission fluxes of biogenic VOCs and fast photochemical production. These seasonal characteristics of SOA compositions appeared similarly in most urban areas in Korea. Atmosphere 2021, 12, x FOR PEER REVIEW 10 of 23 The CMAQ with the AERO7 module considerably overestimated the PM2.5 mass in the winter. As the pcSOA was a dominant species in the winter, the PM2.5 overestimation might be lessened by adjusting the pcSOA concentration by multiplying the pcSOA weighting factor (PWF). As shown in Figure 8, the AERO7 module with the PWF of 1.0 overestimated the observed monthly mean PM2.5 concentrations at Seoul by 12~49%. As the PWF decreased to 0.4 and to 0.0, the modeled PM2.5 concentrations were getting smaller. The modeled monthly mean PM2.5 concentrations best agreed with the observations when the PWF was set to 0.4. In addition, we conducted statistical performance evaluation to search an optimum value of the PWF. As listed in Table 3, the model performances in PM2.5 forecast at four major cities in Korea in 2016 were evaluated by correlation coefficient (R), normalized mean bias (NMB), normalized mean error (NME), probability of detection (POD), false alarm ratio (FAR), and Heidke skill score (HSS), the definition and the usage of which were presented in Section 2.6. As the PWF increased, the R increased indicating performance improvement. However, there was no clear relation between the PWF and the NME. The POD increased with the PWF indicating performance improvement whereas the FAR decreased with the PWF indicating performance decline. The CMAQ with the AERO7 module considerably overestimated the PM 2.5 mass in the winter. As the pcSOA was a dominant species in the winter, the PM 2.5 overestimation might be lessened by adjusting the pcSOA concentration by multiplying the pcSOA weighting factor (PWF). As shown in Figure 8, the AERO7 module with the PWF of 1.0 overestimated the observed monthly mean PM 2.5 concentrations at Seoul by 12~49%. As the PWF decreased to 0.4 and to 0.0, the modeled PM 2.5 concentrations were getting smaller. The modeled monthly mean PM 2.5 concentrations best agreed with the observations when the PWF was set to 0.4. The CMAQ with the AERO7 module considerably overestimated the PM2.5 mass in the winter. As the pcSOA was a dominant species in the winter, the PM2.5 overestimation might be lessened by adjusting the pcSOA concentration by multiplying the pcSOA weighting factor (PWF). As shown in Figure 8, the AERO7 module with the PWF of 1.0 overestimated the observed monthly mean PM2.5 concentrations at Seoul by 12~49%. As the PWF decreased to 0.4 and to 0.0, the modeled PM2.5 concentrations were getting smaller. The modeled monthly mean PM2.5 concentrations best agreed with the observations when the PWF was set to 0.4. In addition, we conducted statistical performance evaluation to search an optimum value of the PWF. As listed in Table 3, the model performances in PM2.5 forecast at four major cities in Korea in 2016 were evaluated by correlation coefficient (R), normalized mean bias (NMB), normalized mean error (NME), probability of detection (POD), false alarm ratio (FAR), and Heidke skill score (HSS), the definition and the usage of which were presented in Section 2.6. As the PWF increased, the R increased indicating performance improvement. However, there was no clear relation between the PWF and the NME. The POD increased with the PWF indicating performance improvement whereas the FAR decreased with the PWF indicating performance decline. In addition, we conducted statistical performance evaluation to search an optimum value of the PWF. As listed in Table 3, the model performances in PM 2.5 forecast at four major cities in Korea in 2016 were evaluated by correlation coefficient (R), normalized mean bias (NMB), normalized mean error (NME), probability of detection (POD), false alarm ratio (FAR), and Heidke skill score (HSS), the definition and the usage of which were presented in Section 2.6. As the PWF increased, the R increased indicating performance improvement. However, there was no clear relation between the PWF and the NME. The POD increased with the PWF indicating performance improvement whereas the FAR decreased with the PWF indicating performance decline. The R, the NMB, the NME, the POD, and the FAR are simple statistical metrics evaluating a partial performance of the model and therefore they may not agree with each other as happened here. This study used the HSS to decide an optimum value of the PWF, because the HSS evaluates the overall performance. The HSS averaged over the four cities was the highest when the PWF was equal to 0.4 or 0.7. The HSS at Seoul, the most populated city in Korea, was the highest when PWF was equal to 0.4. This value of the PWF is the same as the value previously derived using Figure 8. Therefore, the PWF was set equal to 0.4 in this study.
The PM supersite located at Bulkwang in Seoul collected PM 2.5 on a quartz filter every three days and analyzed them for organic carbon (OC). The observed OC at the Bulkwang PM supersite was compared with the modeled OC in Figure 9. The CMAQ calculated organic matter (OM) concentrations instead of OC concentrations. Therefore, the calculated OM concentrations were converted to the OC concentrations by taking account of molecular formula of individual organic aerosol species constituting the OM. As shown in Figure 9, the modeled OC agreed well with the observation except the summer. As the pcSOA dominated except the summer, the model must have accurately estimated the pcSOA concentration with the PWF of 0.4. Atmosphere 2021, 12, x FOR PEER REVIEW 12 of 23 Figure 9. Monthly observed and forecasted organic carbon concentration at the Bulkwang PM supersite.

Evaporation Loss of Nitrate
A national reference method for PM2.5 in many countries including Korea specifies Teflon membrane filter for collection of ambient particles. Teflon membrane filter may cause a significant evaporation loss of particulate nitrate (NO3 − ) during and after the sampling in the warm environment [41,42]. Assuming equilibrium, the evaporation loss (∆NO ) is expressed as where TR is the daily average temperature of the sampled air and Ki is the dissociation constant for ammonium nitrate evaluated at the ambient temperature and the relative humidity for hour i as described by Frank (2006) [43]. The GMAF applied Equation (5) to the temperature and humidity calculated by the WRF to estimate the nitrate evaporation loss. The estimated nitrate evaporation loss was subtracted from the nitrate concentration calculated by the CMAQ to estimate the nitrate concentration left on the filter after evaporation.

Implementation of Bias Adjustment Techniques
The accuracy of air quality forecast may be improved by applying bias adjustment techniques. Both the hybrid bias correction technique and the Kalman filter bias adjustment technique were implemented into the GMAF as postprocessors. The hybrid bias correction technique uses either the difference or the ratio between the forecast and the observation [44]. The former is an additive correction and expressed as The latter is a multiplicative correction and expressed as where C is the corrected forecast value, F is the raw forecast value and O is the observation. Nday denotes the number of previous days representing a training period. Typical values for Nday are the numbers between 1 and 5. Kalman filter bias adjustment technique predicts the future bias by correcting the previous estimate of the bias using the difference between the bias estimate and the forecast error at the previous time. The predicted future bias is used to adjust the forecast. The Kalman filter bias adjustment technique was implemented into the GMAF by following the work by Delle Monache et al. (2006) [22]. Kalman filter bias adjustment technique is One other thing to note is that the model significantly over-predicted the OC in the summer. As most of SOA in this period is the SOA from biogenic VOCs as shown in Figure 7, it is highly likely that the emissions of biogenic VOCs, isoprene and terpenes were overestimated or that the conversion rates from the biogenic VOCs to the SOA were overestimated. More studies are necessary to address this over-prediction problem.

Evaporation Loss of Nitrate
A national reference method for PM 2.5 in many countries including Korea specifies Teflon membrane filter for collection of ambient particles. Teflon membrane filter may cause a significant evaporation loss of particulate nitrate (NO 3 − ) during and after the sampling in the warm environment [41,42]. Assuming equilibrium, the evaporation loss where T R is the daily average temperature of the sampled air and K i is the dissociation constant for ammonium nitrate evaluated at the ambient temperature and the relative humidity for hour i as described by Frank (2006) [43]. The GMAF applied Equation (5) to the temperature and humidity calculated by the WRF to estimate the nitrate evaporation loss. The estimated nitrate evaporation loss was subtracted from the nitrate concentration calculated by the CMAQ to estimate the nitrate concentration left on the filter after evaporation.

Implementation of Bias Adjustment Techniques
The accuracy of air quality forecast may be improved by applying bias adjustment techniques. Both the hybrid bias correction technique and the Kalman filter bias adjustment technique were implemented into the GMAF as postprocessors. The hybrid bias correction technique uses either the difference or the ratio between the forecast and the observation [44]. The former is an additive correction and expressed as The latter is a multiplicative correction and expressed as where C is the corrected forecast value, F is the raw forecast value and O is the observation. Nday denotes the number of previous days representing a training period. Typical values for Nday are the numbers between 1 and 5. Kalman filter bias adjustment technique predicts the future bias by correcting the previous estimate of the bias using the difference between the bias estimate and the forecast error at the previous time. The predicted future bias is used to adjust the forecast. The Kalman filter bias adjustment technique was implemented into the GMAF by following the work by Delle Monache et al. (2006) [22]. Kalman filter bias adjustment technique is sensitive to the error ratio, which determines the relative weighting between the observation and the forecast and needs to be optimized.

Forecast Performance Evaluation Metrics
The methods for forecast performance evaluation are introduced in this section. Selected forecast period and variables for forecast evaluation are also discussed.

Forecast Variables
The GMAF can forecast most of the criteria air pollutants but this study selected the PM 2.5 for performance analysis. As of 2016, 25 PM 2.5 monitoring stations were located in Seoul, 19 stations in Busan, six stations in Gwangju and seven stations in Daegu. Korea uses the PM 2.5 concentration averaged over all the stations in the city of interests to set a policy goal and to calculate the air quality index (AQI). Taking averages over the city improves the spatial representativeness of monitoring data and reduces the associated errors.
Although the recent advances in analytical chemistry have made it possible to measure PM 2.5 mass and composition at time intervals equal to or less than 1 h, most of PM 2.5 toxicity and PM 2.5 management studies still use the daily mean or annual mean concentrations as a metric of PM 2.5 . This study used the daily mean PM 2.5 concentration for forecast performance evaluation.

Performance Evaluation Metrics of Continuous Forecasts
Continuous forecast aims at forecasting the values of continuous variables. Temperatures, wind velocities and air pollutant concentration are good examples of continuous variables. The model performance of continuous forecast was measured by three statistical metrics in this study; R, NMB and NME. The R, correlation coefficient, gives a measure of linear relationship between two variables and ranges between −1 and 1. The R between the forecast and observed values is defined as where O i and F i are the observation and the forecast, respectively, and N is the number of paired data sets. The NMB, normalized mean bias, is the difference between observed and forecasted mean value normalized by the observed mean value and is defined by The NMB is a useful metric for measuring the systematic errors. However, as the number of datasets increases, under-forecasted values tend to offset over-forecasted values to improve the NMB even if the accuracy remains unchanged. The root mean square (RMSE) uses the square of the difference between the observed and forecasted value to measure the errors as defined by The squaring operation makes the RMSE strongly dominated by the largest errors. The RMSE is not recommended for the air quality forecast and analysis because of its high sensitivity to the largest concentrations [10,45].
The normalized mean error (NME) uses the absolute value of the difference instead of the squares of them as defined by The NME gives a measure of the overall forecast error and is used in this study instead of RMSE. Emery et al. (2017) [45] analyzed the model performance evaluation studies published from 2006 to 2012 in the United States and recommended the criteria and goal of the R, the NMB and the NME as listed in Table 4, which is used in this study as a reference. Categorical forecast aims at forecasting the occurrence of a specific event. The event of interests in this study is that the PM 2.5 concentration exceeds the Korean PM 2.5 standard of 35 µg/m 3 . The contingency table used for analysis of the relationship between two categorical variables is presented in Table 5, where the "A" is the number of event forecasts with event occurring. The "B" is the number of event forecasts with event not occurring. The "C" is the number of no-event forecasts with event occurring. Finally, the "D" is the number of no-event forecasts with event not occurring. The A, B, C, D in the contingency table may formulate various performance metrics, including probability of detection (POD), probability of false detection (PODF), false alarm ratio (FAR), success ratio (SR), critical success index (CSI), true skill statistic (TSS), Gilbert skill score (GS), Heidke skill score (HSS), bias (B), and accuracy (A) [46].
The POD ranges from 0 to 1 with 1 being a perfect forecast and is defined by The POD considers missed events only, neglecting false alarms and therefore it can be artificially improved by over-forecasting. The FAR complements the POD in this aspect and is defined by The FAR ranges from 0 to 1 with 0 being a perfect forecast. The FAR considers false alarm events only, neglecting missed events and therefore it can be artificially improved by under-forecasting.
The POD and the FAR are examples of simple skill metrics. The problem of simple skill metrics is that a good score does not necessarily mean a good overall forecasting performance. One can overcome this problem by combining two or more simple skill metrics. This kind of skill metrics, called here as complex skill metrics, are the CSI, the TSS and the HSS. The CSI and the TSS may not work properly for rarely occurring events. As exceedance of the PM 2.5 standard, the event in this study, may occur rarely in certain seasons and regions, we selected the HSS as a complex skill metric.
The HSS is the fraction of accurate forecast with correction of random forecasts as defined by where E is the number of correct random forecasts as defined by The HSS ranges from −∞ to 1 with 1 being a perfect forecast. The HSS above zero means that the forecast is better than a random guess and the HSS above 0.5 is deemed a good score.
Studies on categorical forecast performance evaluation of PM 2.5 forecast were much more limited than studies on continuous forecast performance evaluation [47]. The POD for the event of exceeding the 8-h O 3 standard in the Northeastern United States ranged between 0.07 and 0.37 whereas the FAR, often larger than the POD, ranged between 0.64 and 0.76 [48]. The POD for the event of exceeding the 24-h PM 2.5 standard in the Eastern United States was 0.28 whereas the FAR was 0.35 [49]. The POD for the event of exceeding the 24-h PM 2.5 standard in Korea ranged between 0.1 and 0.39 whereas the FAR ranged between 0 and 0.64 [50].

Results and Discussion
The performance evaluation results of the GMAF are presented in this section. Test results of performance improvement attained by incorporating the nitrate evaporation and the bias adjustment techniques are also presented.

Performance Evaluation of the Base Forecast
The GMAF was used to forecast the PM 2.5 for the year 2016. All the modifications described in Section 2.4 were included except nitrate evaporation loss. Table 6 shows the seasonal and yearly values of performance metrics, where the seasons were defined as winter (December, January, February), spring (March, April, May), summer (June, July, August), and fall (September, October, November). The average PM 2.5 concentrations in the winter and the spring were higher than those in the summer and the fall in all the four selected cities, which are very typical of the PM 2.5 concentrations in Korea. The annual values of R were higher than 0.78 in all four cities, satisfying the goal set in Table 4. The winter values of R were the higher than the other seasons and the summer values of R were lower than the other seasons except Busan. The NMBs at Seoul were positive in all the seasons, indicating over-forecasting, whereas those at Busan were all negative, indicating under-forecasting. Out of total 16 cases (four seasons and four cities), half of the cases had the NMB value less than ±10%, achieving the goal set in Table 4 and the other half of cases had the NMB value less than ±30%, satisfying the criteria also set in Table 4. The seasonal NME ranged between 19% and 27%, achieving the goal set in Table 4. Unlikely to the R, the NME did not exhibit a distinctive seasonal variation.
The number of days exceeding the 24-h PM 2.5 standard, which is the frequency of occurrence of events, ranged between 17 and 31 days during the winter and the spring because of severe PM 2.5 pollution in Korea. The frequency of occurrence of events in this study is several times higher than that in the literature quoted in Section 2.6.3. It may be noted here that categorical forecasts perform better for the frequent events than for the rare events. During the winter, the skill scores were very good; the POD were 0.68~1.0 and exceeded the FAR by 0.41~0.59, and HSS were 0.58~0.68. Although the skill scores went down in the spring and fall, all three skill scores were still deemed fairly good.
The categorical forecast performed poorly during the summer compared to the other seasons because of two reasons. Firstly, as noted above, the R values and the NME values in the summer were generally lower than those in the other seasons, indicating the GMAF might not perform well in the summer. Secondly, the events occurred less frequently in the summer due to lessened PM 2.5 pollution in Korea, which might worsen the categorical skill scores measured by the POD, the FAR and the HSS.
As discussed above, both the GMAF performed the best during the winter and worst during the summer. In order to investigate the effects of seasons on the forecast performance, the modeled seasonal PM 2.5 chemical compositions are displayed in Figure 10. The nitrate was a dominant PM 2.5 species in the winter. The PM 2.5 nitrate in the urban area is rather locally produced than regionally transported [51]. The local emissions of nitrogen oxides (NO x ) were well known and the homogeneous and heterogeneous chemical pathways were well defined such that the GMAF was able to accurately predict the PM 2.5 nitrate concentrations and subsequently the PM 2.5 mass concentrations during the winter.

Modeling of Nitrate Evaporation Loss
The PM supersite located at the Bulkwang in Seoul measured the PM2.5 mass and PM2.5 nitrate using Teflon membrane filter. Teflon membrane filter causes a significant evaporation loss of particulate nitrate during and after the sampling in the warm environment. Here, the nitrate before evaporation loss is defined as as-is nitrate, whereas the nitrate left on the filter after evaporation loss is defined as retained nitrate.
The monthly mean concentrations of the retained nitrate at the Bulkwang site were calculated by applying Equation (5) to the temperature and the humidity forecasted by WRF and the nitrate concentration modeled by CMAQ. As displayed in Figure 11, both the as-is and retained PM2.5 nitrate concentration were the highest during the winter months and the lowest during the summer months because the ratios of nitrate aerosol to nitric acid gas (NO3 − /HNO3) decrease with temperature increase. In the summer months, the nitric acid mainly stayed in the gas phase due to the warm temperature. It may be noted that the nitric acid gas is readily removed by wet and dry deposition to have shorter lifespan than the nitrate aerosol.
The difference between as-is nitrate and retained nitrate is the evaporation loss of nitrate taking place during and after Teflon filter sampling. The evaporation loss is less than 10% of as-is nitrate during the winter months whereas it was more than 80% of as-is nitrate during the summer months. Therefore, the retained nitrate was negligibly small by being less than 0.5 μg/m 3 during the summer months. During the summer, the PM 2.5 mass was dominated by the OM, the most of which was the SOA from biogenic VOC as shown in Figure 7. There was large uncertainty in estimating biogenic VOC emission and in modeling aerosol formation processes from biogenic VOCs, which probably degraded the model performance.

Modeling of Nitrate Evaporation Loss
The PM supersite located at the Bulkwang in Seoul measured the PM 2.5 mass and PM 2.5 nitrate using Teflon membrane filter. Teflon membrane filter causes a significant evaporation loss of particulate nitrate during and after the sampling in the warm environment. Here, the nitrate before evaporation loss is defined as as-is nitrate, whereas the nitrate left on the filter after evaporation loss is defined as retained nitrate.
The monthly mean concentrations of the retained nitrate at the Bulkwang site were calculated by applying Equation (5) to the temperature and the humidity forecasted by WRF and the nitrate concentration modeled by CMAQ. As displayed in Figure 11, both the as-is and retained PM 2.5 nitrate concentration were the highest during the winter months and the lowest during the summer months because the ratios of nitrate aerosol to nitric acid gas (NO 3 − /HNO 3 ) decrease with temperature increase. In the summer months, the nitric acid mainly stayed in the gas phase due to the warm temperature. It may be noted that the nitric acid gas is readily removed by wet and dry deposition to have shorter lifespan than the nitrate aerosol.
The difference between as-is nitrate and retained nitrate is the evaporation loss of nitrate taking place during and after Teflon filter sampling. The evaporation loss is less than 10% of as-is nitrate during the winter months whereas it was more than 80% of as-is nitrate during the summer months. Therefore, the retained nitrate was negligibly small by being less than 0.5 µg/m 3 during the summer months.  Figure 12 compares the observed nitrate with the modeled nitrate between April and October during which Equation (5) predicted a significant evaporation loss of nitrate as shown in Figure 11. As the nitrate was measured every three days, the modeled nitrate concentrations were also selected every three days in the Figure 12. The modeled retained nitrate agreed well with the observed nitrate in April, August, September, and October, indicating that Equation (5) calculated the nitrate evaporation loss accurately. The retained nitrate concentrations were kept below 5 μg/m 3 to lessen the nitrate overestimation appeared in the as-is nitrate concentration. However, the observed nitrate was closer to the modeled as-is nitrate than the modeled retained nitrate about half of May, June and July. This underestimation of the retained nitrate is probably due to inaccurate estimation of the temperature to which the nitrate evaporation loss is sensitive.  Table 7 shows the effect of inclusion of nitrate evaporation loss on the forecast performance in the spring and the fall during which the nitrate evaporation loss had a discernible effect on the PM2.5 mass concentration. The CTMs including the CMAQ are exposed to high degrees of uncertainty in modeling physical and chemical processes and in preparing input data such as emissions, meteorological parameters. Enhancement of one module does not necessarily mean improvement of the forecast performance. As the nitrate evaporation loss was included, the NMB decreased as expected. The R decreased in  Figure 12 compares the observed nitrate with the modeled nitrate between April and October during which Equation (5) predicted a significant evaporation loss of nitrate as shown in Figure 11. As the nitrate was measured every three days, the modeled nitrate concentrations were also selected every three days in the Figure 12. The modeled retained nitrate agreed well with the observed nitrate in April, August, September, and October, indicating that Equation (5) calculated the nitrate evaporation loss accurately. The retained nitrate concentrations were kept below 5 µg/m 3 to lessen the nitrate overestimation appeared in the as-is nitrate concentration. However, the observed nitrate was closer to the modeled as-is nitrate than the modeled retained nitrate about half of May, June and July. This underestimation of the retained nitrate is probably due to inaccurate estimation of the temperature to which the nitrate evaporation loss is sensitive.  Figure 12 compares the observed nitrate with the modeled nitrate between April and October during which Equation (5) predicted a significant evaporation loss of nitrate as shown in Figure 11. As the nitrate was measured every three days, the modeled nitrate concentrations were also selected every three days in the Figure 12. The modeled retained nitrate agreed well with the observed nitrate in April, August, September, and October, indicating that Equation (5) calculated the nitrate evaporation loss accurately. The retained nitrate concentrations were kept below 5 μg/m 3 to lessen the nitrate overestimation appeared in the as-is nitrate concentration. However, the observed nitrate was closer to the modeled as-is nitrate than the modeled retained nitrate about half of May, June and July. This underestimation of the retained nitrate is probably due to inaccurate estimation of the temperature to which the nitrate evaporation loss is sensitive.  Table 7 shows the effect of inclusion of nitrate evaporation loss on the forecast performance in the spring and the fall during which the nitrate evaporation loss had a discernible effect on the PM2.5 mass concentration. The CTMs including the CMAQ are exposed to high degrees of uncertainty in modeling physical and chemical processes and in preparing input data such as emissions, meteorological parameters. Enhancement of one module does not necessarily mean improvement of the forecast performance. As the nitrate evaporation loss was included, the NMB decreased as expected. The R decreased in  Table 7 shows the effect of inclusion of nitrate evaporation loss on the forecast performance in the spring and the fall during which the nitrate evaporation loss had a discernible effect on the PM 2.5 mass concentration. The CTMs including the CMAQ are exposed to high degrees of uncertainty in modeling physical and chemical processes and in preparing input data such as emissions, meteorological parameters. Enhancement of one module does not necessarily mean improvement of the forecast performance. As the nitrate evaporation loss was included, the NMB decreased as expected. The R decreased in all the four cities for the spring whereas it increased in Gwangju and Daegu for the fall. The NME increased in most cases. Therefore, it may be concluded that the performance of the continuous forecast was not improved although the nitrate evaporation loss was included. Moreover, the categorical forecast measured by the POD, the FAR and HSS was declined in most cases considered here.

Application of Bias Adjustment Technique
The bias adjustment techniques described in Section 2.5 were applied to the PM 2.5 forecasts at the four cities for the year 2016. We tested different values of the number of previous days (Nday) for additive and multiplicative bias correction and different values of the error ratio for the Kalman filter bias adjustment, and chose 2 as the number of previous days and 0.5 for the error ratio.
The additive and multiplicative bias correction techniques were designed to remove the biases by algebraic manipulation. They successfully eliminated the biases at all the four cities as shown in Table 8. However, the improvements of biases did not propagate to the other performance evaluation metrics. They did not generally improve the R and the NME except for Seoul. Moreover, they generally degraded the POD and the HSS.
Although the Kalman filter bias adjustment technique was less efficient in reducing the NMB than the additive and multiplicative bias correction techniques, it improved all the other performance metrics significantly. The averaged values over the four cities were improved by 5% for the R, 13% for the NME, 11% for the POD, 8% for the FAR, and 10% for the HSS.

Conclusions
The GMAF was developed by utilizing the currently available multiscale air quality modeling system and the global forecast and analysis dataset. The GMAF implemented the grid nudging based FDDA onto the CMAQ for the first time to take advantage of data assimilation done in the global modeling. Although the current forecast verification was limited to the next day forecast and analysis, the forecasting period can be extended up to five days permitted by the global forecast and analysis model without a significant loss of accuracies.
The GMAF lessened the tendency of underestimating the PM 2.5 concentrations in the case of little precipitation by improving the below-cloud scavenging coefficient estimation. The CMAQ-5.3.1 with carbon bond version 6 mechanism and AERO7 module significantly enhanced the accuracy of nitrate formation rate calculation by updating homogenous reaction rate constants and heterogeneous reaction mechanism. The SOA formation mechanism was also expanded by adding pcSOA, organic nitrates from monoterpene oxidation and SOA from glyoxal and methylglyoxal. The expansion of the SOA formation mechanism amended the underestimation problem of SOA frequently found in the previous version. However, a comparison with the observations indicated that the pcSOA in the winter and the SOA from biogenic VOCs in the summer might be overestimated. A more refined estimation of the pcSOA emissions and more accurate modeling of biogenic VOCs are recommended as future studies.
A severe PM 2.5 episode often arises in the winter and the spring in Korea. The Korea government agency relies on the numerical forecast during these PM 2.5 season in issuing a PM 2.5 advisory or warning to the people. The GMAF performed very well especially in the winter suiting the needs of the government. In the winter, the R between the observed and forecasted PM 2.5 was 0.87, the POD for event exceeding the national PM 2.5 standard was 1.00, and the FAR was only 0.48 at Seoul. The GMAF also showed good performance at the other three cities, Busan, Gwangju and Daegu.