Statistical Guidance on Seasonal Forecast of Korean Dust Days over South Korea in the Springtime

Asian dust is a meteorological event in which dust particles are uplifted into the air by the strong surface wind in the source regions in the northern China and Mongolia, and are transported far away and deposited on the ground. Asian dust events that deposit dust in South Korea (called Korean dust events) may lead to severe economic and social damages. According to the National Institute of Meteorological Research (NIMR) in the Korea Meteorological Administration (KMA), the frequency and variation of Asian dust events are increasing gradually.

The source regions of Asian dust, that affect the Korean peninsula, cover most of northern China and Mongolia including the arid and semi-arid area of sand desert. (Zhang et al., 2008) divided the source regions of Asian which were the main contributors to long-lived dusts in northern China and Mongolia, into five regions. (Lim and Chun, 2006) divided these source regions into three regions (Gobi region, Inner Mongolia, and the northeast China) and showed that the source regions of Asian dust that affect the Korean peninsula are gradually extending eastward. (Chun et al., 2006) also showed that source regions of Asian dust have been expanding to Northeast China year by year because of the desertification in Northeast China due to drought, overgrazing and change of atmospheric circulation. (Kurosaki and Mikami, 2003) showed that the Asian dust events increased remarkably in the eastern part of the Asian Continent after 2000. (Chun and Lim, 2004) concluded that the intensity and occurrence of Asian dust have been increased gradually since 2000 and they have started to occur earlier in the year than before. (Tian et al., 2007) examined the relationship between Asian dust frequency observed in Japan and that in northern China with long-term observations and the National Center for Environmental Prediction in National Center for Atmospheric Research (NCEP/NCAR) reanalysis data.

The KMA began to announce the seasonal forecast of Asian dust events since 2001. Since 2008, objective models have been utilized for the seasonal forecast of Korean dust days over South Korea. The forecast modeling is changed every year in order to improve forecast accuracy. (Sohn et al., 2009a) applied a principal component regression model to the quantitative forecast and categorical forecast of Korean dust days (KDD). In this paper, KDD means the spatial average of monthly Korean dust days observed at 28 stations in South Korea. However, the skill scores were low and because the principle components were linear combinations of several kinds of data, they were difficult to explain. (Sohn et al., 2010) and (Park et al., 2011) considered the large-scale climate indices as additional predictors. Their accuracy of forecast is much improved; however, it was well-fit before the year 2000 but not after 2001.

This study aimed at improving the seasonal forecast model of KDD for each month in the springtime separately. In this study, we considered three kinds of data to examine the relationships among dust deposition, dust emission, and dust transport based on statistical analysis. The first one was meteorological data to estimate dust emission. Following the method of (Lim and Chun, 2006), we used meteorological factors which are considered to be important factors for the occurrence of Asian dust, in our study. As potential predictors, we used the NCEP/NCAR reanalysis data (hereafter, referred to as NNR) for ground temperature, rainfall amount, snowfall amount and ground wind speed. The second one was large-scale climate indices for the dust transport from source regions to the Korean peninsula. (Schwing et al., 2002) and (Kim et al., 2008) indicated that atmospheric and oceanic teleconnections govern many climate variations, the cause and effect relationships between these teleconnections are not well understood. In those cases, large-scale climate indices may be useful to explain the relationships. We examined whether large-scale climate indices can explain the dust transport from source regions to South Korea. The third one is the KDD.

In this study, we used a ternary forecast (below normal, normal and above normal) mode. It is well known that logistic regression model is useful to generate multinomial distribution when the target variable has a categorical value including ordinal scaling data (Myers et al., 2002). (Sohn et al., 2009b) applied a logistic regression model to the categorical forecast of heavy snowfall. It is easy to explain the effect of predictors via odds ratio using this type of model. In this study, the logistic regression model generated probability forecast in the form of a trinomial distribution (p_0, p_1, p₂) where p₀ is P (below normal), p₁ is P (normal) and p₂ is P (above normal). The trinomial distribution was converted into a ternary forecast using thresholds.

Forecast guidance consists of the estimated ordinal logistic regression model and the conversion algorithm to generate categorical forecast. The predictability of forecast guidance is evaluated by assessment measures called skill scores. Which can be of various types (von Storch and Zwiers, 1999; Sohn and Park, 2008). In this study, skill scores for the binary case and ternary case were used to assess the forecast guidance and find the optimal guidance.

In section 2, three kinds of data are presented and analyzed. Temporal variations in KDD and NNR are examined. In section 3, the potential predictors are selected based on correlation analysis between KDD and predictors. Skill scores are introduced and forecast strategy and guidance are explained. The forecast guidance is proposed for each month separately and the results are summarized in section 4.

For this study, we used KDD over South Korea for a period of 30 years (1981-2010). Table 1 shows the 30-year mean of KDD for each month, from which it is clear that KDD mainly occurred (about 85%) in the springtime. Our study focused on KDD in the springtime. Figure 1 illustrates time series plots of KDD for three month (March-May). It is evident that the change patterns of KDD were different from each other. And the three KDD series were not significantly correlated. It led us to develop the forecast guidance for each month separately.

The interdecadal variations in KDD series given in Table 2, showed that the mean levels (0.44, 0.99, and 3.74) of KDD increased in March. As shown in Fig. 1, KDD increased remarkably after the middle of 1990s, and its significant slope was 0.145 days per year. The mean levels (1.39, 2.94, and 2.72) of KDD increased further in April in the 1990s and slightly decreased in the 2000s. However, almost no change occurred in the mean level (1.05, 0.97, and 1.00) of KDD in May. The changes in standard deviation were similar to the corresponding mean levels.

Figure 1.  Time series plots of Korean dust days over South Korea for March (dot), April (square) and May (circle).

Figure 2.  Three source regions of Asian dust.

Figure 3.  30-year means of NNR in three regions (A, B and C) for each factor and month.

The target value in this study was a ternary forecast which was classified following the classification rules in Table 3, for each month in the springtime separately. KDD said to be below normal if KDD is less than (M-0.5×STD), above normal if greater than (M + 0.5×STD), and normal if otherwise, where M and STD are the 30-year mean and 30-year standard deviation respectively.

According to (Lim and Chun, 2006), source regions of Asian dust that affected the Korean peninsula in the springtime, are divided into three source regions (Fig. 2). The three regions are located in dry arid region (A, Gobi region: 35°-45°N, 100°-110°E), semi-arid region (B, Inner Mongolia region: 40°-45°N, 110°-120°E), and cultivated region (C, the northeast part of China: 40°-50°N, 120°-125°E).

The regional mean of NNR for monthly mean of ground temperature, monthly rainfall amount, monthly snowfall amount and monthly mean of ground wind speed over three source regions for a period of 30 years (1980-2009) were used as potential predictors. Figure 3 shows the 30-year means of NNR in three regions for each month. Ground temperatures in the three regions were similar to each other. Rainfall and snowfall amount in C were more than those in A and B, while the ground wind speed in B was higher than that in the others.

Table 4 shows the mean interdecadal variation in anomaly data of NNR. The ground temperature was found to have an increasing trend in all regions. Rainfall and snowfall amount showed a decreasing trend in all regions. Although the ground wind speed appeared to be decreasing in all regions, it did not show any significant trend. NNR indicated that all regions became drier and warmer, which might have increased floating dusts over the source regions; we considered these four factors as predictors for KDD.

We also used nine climate indices (AO, GML, MEI, NOI, ONI, PDO, PNA, SOI, and WPO; see in Table 5) as additional predictors. They were obtained from the Climate Diagnostic Center (CDC) and Climate Prediction Center (CPC) in National Oceanic and Atmospheric Administration (NOAA) and updated on a monthly basis. The AO is an index of the pressure gradient between the polar and subpolar regions of the Northern Hemisphere. The GML is the global mean land/ocean temperature index. The MEI is calculated as the first principal component of six observed variables fields (SLP, zonal and meridional components of the surface wind, SST, SAT, total cloudiness fraction of the sky) over the tropical Pacific. The NOI is an index of climate variability based on the difference in SLP anomalies at the North Pacific High and near Darwin. The ONI is a three-month running mean of SST anomalies in the Niño3.4 region (5°N-5°S, 120°-170°W) based on the 1971-2000 base period. The PDO is the first principal component of monthly SST anomalies in the North Pacific Ocean. The PNA is an index derived from the formula of (Wallace and Gutzler, 1981) using 500-hPa geopotential height values. The SOI is an index of climate variability based on the difference in SLP anomalies at Tahiti and Darwin. The WPO identified through the application of a principal-component analysis involves changes to the wind pattern in the northwestern and north-central North Pacific.

In NIMR (2011), forecast models were estimated using 30-year data (1981-2010), which were well-fit to the data before 2001, but did not show a good fit after 2001; this was because the patterns of KDD series changed after the middle of 1990s as shown in Fig. 1. To improve predictability, we decided to use 20-year data (1991-2010) for forecast model and 30-year data for classification rules.

As potential predictors, we use 12-month data (January-December in the previous year) of each factor because the KMA usually announces seasonal forecast for KDD in the springtime in February. The set of potential predictors contained as so many as 252 variables (3 regions×12 months×4 factors + 12 months×9 indices). To simplify the process, at first, we selected the secondary potential predictors (Table 6) among 252 variables, which were significantly correlated with KDD. Predictors were named by three components (region, factor, and month). For example, AT7 means the monthly mean of ground temperature over Region A in July in the previous year and AO3 means the Arctic oscillation index in March in the previous year.

The ternary forecast was assessed based on a 3×3 contingency table (see Table 7). We considered the thresholds that had good skill scores for both binary forecast and ternary forecast, since the above normal case was more serious than the normal or the below normal case. In this study, binary forecast was generated from ternary forecast in order to satisfy the concordance between two forecasts with merging the normal and the below normal cases into the normal cases.

In the binary forecast, we considered three skill scores based on a 2×2 contingency table. They were the hit rate (HR2), the probability of detection (POD2) and the false alarm rate (FAR2), and were computed by the following equations using the notation in Table 8.

HR2= (n₀₀ + n₀₁ + n₁₀ + n₁₁ + n₂₂) / N,

POD2= n ₂₂ / (n₂₀ + n₂₁ + n₂₂),

FAR2= n₀₂ /(n₀₂ + n₁₂+ n₂₂)

where N is the number of total cases, HR2 is the exact forecast rate, POD2 is the conditional probability that above normal is forecasted given that above normal is observed, and FAR2 is the conditional probability that above normal is not observed given that above normal is forecasted. It is natural that the higher the values for HR2 and POD2 and the lower value for FAR2, the better the outcomes.

For the assessment of ternary forecast, we compute skill scores based on the 3×3 contingency table. (Burrows, 1991) assigned weights to the multicategorical forecast cases, so that a forecast that is different from the predictand in one category would be a better forecast than a forecast which is different in two categories, and so on. Following Burrows' procedure, we also used weights in Table 7, and computed the hit rate (HR3) and the weighted hit rate (wHR3) defined by the following equations using the notation in Table 7.

HR3= (n₀₀ + n₁₁ + n₂₂)/N,

wHR3= (n₀₀ + 0.5 n₀₁ + 0.5 n1₀ + n₁₁ + 0.5 n₁₂ + 0.5 n₂₂ + n₂₂) / N

It is natural that the higher the values for HR3 and wHR3, the better the outcomes. Note that HR3 is the hit rate in a 3×3 table whereas HR2 is the hit rate in a 2×2 table.

To check the generalization ability of forecast guidance, we performed a leave-one-out cross-validation. The results of model validation were examined by the above skill scores on the basis of the 3×3 contingency table.

Ternary forecast was generated by a two-step process as follows. At first, an ordinal logistic regression model was utilized to generate trinomial distribution (p₀, p₁, p₂) using the following equations.

e₀= exp (b₀ + X); e₁= exp(b₁ + X);

p₀= e₀/(1+e₀); p₁ = e₁/(1+e₁) - p₀;

p₂= 1 - p₀ - p₁

where b₀ and b₁ are constants, and X is a linear combination of final predictors.

Second, the generated trinomial distribution was converted into a ternary forecast. To generate a ternary forecast, we usually selected the category with the highest probability in the generated trinomial distribution. However, the threshold method might be better than the maximum probability method. For example, the skill scores of the maximum probability method were HR2 = 0.85, POD2 = 0.71, FAR2 = 0.17, HR3 = 0.75, wHR3 = 0.87, whereas those of the threshold method are HR2 = 0.90, POD2 = 1.00, FAR2 = 0.22, HR3 = 0.85 and wHR3 = 0.93 (Table 9). In this study, the threshold method was preferred to the maximum probability method. Two thresholds, T1 and T2, were used. A conversion of trinomial distribution (p₀,p₁,p₂) to a ternary forecast (0 = below normal, 1 = normal, 2 = above normal) were made using the followings:

IF p₀ >T1 THEN

forecast = 0;

ELSE IF p₁/(p₀+p₂) > T2 THEN

forecast = 1;

ELSE forecast = 2

Varying two thresholds from 0 to 1 independently, the above skill scores were computed for each pair of two thresholds and compared so as to determine the optimal thresholds. The proposed forecast guidance consisted of the estimated logistic regression model and the conversion algorithm using two optimal thresholds.

We developed the statistical forecast guidance on ternary KDD in the springtime for each month separately. A three-grade ordinal logistic regression model was applied to generate a trinomial distribution using the final predictors given in Table 6. The optimal thresholds were determined based on skill scores. Some results were obtained in the following way.

First, the variable selection was performed by Fisher's scoring method offered by the statistical package called SAS. The final predictors, selected at the significance level 0.05, were the following: BS12, NOI7, and SOI4 for March; CT7, BR5, and NOI1 for April; and CT4, AR12, and AWS11 for May.

Second, the optimal thresholds were determined based on skill scores, which were 0.3 and 0.65 for March; 0.2 and 0.6 for April; and 05 and 0.5 for May.

Third, forecast guidance generated a binary forecast and a ternary forecast together. The proposed guidance for March was as follows:

x = 0.1223× BS12 - 1.8894× NOI7 - 1.918× SOI4;

e₀ = exp(-8.5420 + x); e₁ = exp(-3.7485 + x);

p₀=e₀ / (1 + e₀);

p₁=e₁ / (1 + e₁)-p₀;

p₂ = 1 - p₀-p₁.

IF p₀ > 0.3 THEN

Ternary_forecast=0;

ELSE IF p₁/(p₁+p₂) > 0.65 THEN

Ternary_forecast =1;

ELSE Ternary_forecast =2;

IF Ternary_forecast=0 OR Ternary_forecast=1

THEN Binary_forecast=1;

ELSE Binary_forecast=2.

The proposed guidance for April was as follows:

x = -3.5088× CT7 - 11.1385× BR5 +2.2744× NOI10;

e₀ = exp(85.8864+ x); e₁ = exp(88.0663+ x);

p₀=e₀/ (1 + e₀);

p₁=e₁ / (1 + e₁)-p₀.

P₂ = 1 - p₀-p₁;

IF p₀ > 0.2 THEN

Ternary_forecast=0;

ELSE IF p₁/(p₁+p₂) > 0.6 THEN

Ternary_forecast =1;

ELSE Ternary_forecast =2;

IF Ternary_forecast=0 OR Ternary_forecast=1 THEN

Binary_forecast=1;

ELSE Binary_forecast=2.

The proposed guidance for May was as follows:

x = 0.9213× CT4+41.0877× AR12-1.9697× AWS11;

e₀ = exp(x); e₁ = exp(4.176 + x);

p₀=e₀/(1 + e₀);

p₁=e₁/(1 + e₁)-p₀;

p₀ = 1 - p₀ -p₁;

IF p₀ > 0.5 THEN

Ternary_forecast=0;

ELSE IF p₁/(p₁+p₂) > 0. 5 THEN

Ternary_forecast =1;

ELSE Ternary_forecast =2;

IF Ternary_forecast=0 OR Ternary_forecast=1 THEN

Binary_forecast=1;

ELSE Binary_forecast=2;

Forth, the results of binary forecast were summarized in 2×2 contingency tables (Table 10). According to contingency tables, the computed skill scores were as follows: HR2 = 0.9, POD2 = 1.0, FAR2 = 0.22 for March; HR2 = 0.9, POD2 = 1.0, FAR2 = 0.2 for April; and HR2 = 1.0, POD2 = 1.0, FAR2 = 0.22 for May. Table 11 shows the results of cross-validation for binary forecast as follows; HR2 = 0.85, POD2 = 0.71, FAR2 = 0.16 for March; HR2 = 0.85, POD2 = 0.75, FAR2 = 0.14 for April; and HR2 = 1.0, POD2 = 1.0, FAR2 = 0.22 for May. These scores indicate that each forecast guidance may be reliable and useful.

Fifth, the results of ternary forecast were summarized in the 3×3 contingency tables (Table 12). Based on cross tables, the computed skill scores were as follows: HR3 = 0.85, wHR3 = 0.93 for March; HR3 = 0.85, wHR3 = 0.93 for April; and HR3 = 0.9, wHR3 = 0.95 for May. Table 13 shows the results of cross-validation for ternary forecast and their skill scores were as follows: HR3 = 0.75, wHR3 = 0.88 for March; HR3 = 0.8, wHR3 = 0.9 for April; and HR3 = 0.9, wHR3 = 0.95 for May. These scores indicate that the proposed forecast guidance may be reliable.

Sixth, we applied the proposed models to forecast KDD in the years 2011 and 2012. The results are summarized in Table 14. In 2011, the actual category and forecasted category of KDD were (normal, above normal), (below normal, below normal) and (above normal, above normal) for March, April, and May, respectively. The rate of correct forecast was 67%. In 2012, they are (below normal, normal), (below normal, below normal) and (below normal, normal) for March, April and May respectively. Unfortunately, the rate of correct forecast was only 33% in this case. In April and May 2012, North Pacific high enlarged unusually to the Korean peninsula, which might have suppressed dust transport to South Korea.

To develop a separate ternary forecast of KDD for each month in the springtime, we considered the following strategy. First, we used 20-year data (1991-2010). Second, various skill scores (HR2, POD2, FAR2, HR3, and wHR3) were considered for the assessment of binary and ternary forecasts. Third, we used two thresholds to generate categorical forecast. Fourth, we propose the forecast guidance, which can be used to generate both the binary and the ternary forecasts together, in order to satisfy the concordance between two modes for each month.

The skill scores of proposed forecast guidance were found to be better than before. HR3 scores were 0.85 for March, 0.85 for April and 0.9 for May. We recommend that binary forecast might be more useful than ternary forecast because above normal cases are more important than normal and below normal cases.

In the future, we will investigate the relationships among large-scale climate indices in order to find how they affect Asian dust events.