Implications from Subseasonal Prediction Skills of the Prolonged Heavy Snow Event over Southern China in Early 2008

During early 2008, southern China experienced an exceptional prolonged heavy snow event (PHSE) accompanied by a severe cold surge and freezing rain, leading to massive economic losses and numerous casualties. This unprecedented extreme event severely damaged public infrastructure and caused an inordinate amount of traffic congestion in many provinces in southern China (Zhao et al., 2008).

Given its severity, huge efforts have been devoted to understanding the physical mechanism of this disaster. For example, it is suggested that during this PHSE, the Siberian blocking high over the mid−high latitudes persisted for three weeks at around 65°E, bringing cold, dry air to southern China and causing the extremely cold event (Bueh et al., 2011). Meanwhile, abundant water vapor was persistently transported from the tropical oceans to the subfreezing areas in southern China, leading to the prolonged heavy snow event(Tao and Wei, 2008; Zhou et al., 2009). The water vapor was advected by anomalous southerly winds, which was induced by the tropical diabatic cooling that occurred when the suppressed tropical convection propagated to the Maritime Continent (Hong and Li, 2009; Ma et al., 2011). In addition, the autumn Arctic sea ice (Wu et al., 2011; Li and Wu, 2012), the westerly jet stream (Wen et al., 2009), and the anomalous Tibetan Plateau warming (Bao et al., 2010) may also have played roles in the formation and maintenance of this extreme event. The PHSE in early 2008 caused severe impacts at an unprecedented scale. However, due to the complicated nature of the influence factors and its long duration (exceeding three weeks), accurate prediction of this event was quite challenge.

Forecasting information of such a persistent extreme event not only requires correct climate mean background condition, but also detailed intraseasonal evolution of the event. For the sake of successful disaster mitigation and prevention, accurate prediction of this type of extreme event 10−30 days in advance is imperative and is the main target of the subseasonal to seasonal (S2S) prediction project (Vitart et al., 2017).

Previous studies have evaluated the skill in predicting Northern winter climates on the subseasonal timescale. For example, Li et al. (2017) indicated that even though the S2S model shows encouraging skill beyond two weeks in predicting East Asia cold surges, the forecast Siberian High and associated anomalous northerly winds are considerably weaker than observed. Zhou et al. (2019) evaluated the subseasonal forecast skill of 2-m air temperature in winter over China based on 11 S2S models and suggested the European Centre for Medium-Range Weather Forecasts (ECMWF) model performed best, with useful forecast skill up to four pentads, while the model from China Meteorological Administration (CMA) exhibited useful prediction skill to about three pentads. Lyu et al. (2019) suggested that the Madden−Julian Oscillation (MJO) and North Atlantic Oscillation (NAO) are important subseasonal predictability sources for the winter cold surges in China, but the ECMWF model can only reproduce the relationship between MJO/NAO and cold surges up to four days. Xiang et al. (2020) reported that the ECMWF model exhibit skillful prediction of cold extremes over a great fraction of the Northern Hemisphere with a lead time of 2−4 weeks. They also found that the winter climate modes over Eurasia are more predictable than those of North America. Li et al. (2020) suggested that the ECMWF and CMA models can skillfully forecast Tibetan Plateau snow cover within a 2-week lead time during winter, but they all overestimated the area of Tibetan Plateau snow cover, leading to a cold bias of the surface air temperature. Cui et al. (2021) suggested the 11 S2S models could skillfully forecast the winter surface temperature intraseasonal mode over mid-to-high latitudes of Eurasia with a 10−20-day lead time, and the maximum attainable predictability of this mode can up to 25 days.

However, it remains to be determined: What is the operational S2S models’ prediction skill for the disaster events such as PHSE over southern China in 2008? Are there any universal predictability sources for such disastrous events? These are elusive issues remaining as unresolved. Building upon the current advances in S2S models, this article examines the prediction skill of the PHSE over southern China in early 2008 using three S2S prediction models (ECMWF and two CMA models) and aims to detect the most useful predictability sources for improving the S2S prediction skill of such events. The remainder of the paper is organized as follows: Section 2 introduces the data and methods; section 3 evaluates the skills of the three models in forecasting the PHSE and examines the possible reasons behind their forecast skills. Summary and discussion are provided in section 4.

Both of the S2S databases from ECMWF and CMA are employed in the present study. The ECMWF model is considered as the best operational model (Zhou et al., 2019), while the CMA model is the independently developed Chinese operational model. The CMA includes two different models, CMA1.0 and CMA2.0. The CMA2.0 is operated in the second phase of the S2S Projection, and includes advancements in parameterization, assimilation and initialization schemes over those in CMA1.0 (Liu et al., 2017, 2021). The ECMWF ensemble model has 11 members, while the CMA models have 4 members.

The CMA1.0 model is initialized every day, while the ECMWF and CMA2.0 model are initialized twice a week (every Monday and Thursday). For a fair comparison, we rearranged the ECMWF and CMA2.0 model data to be an everyday-lead based on the method proposed by Yang et al. (2018a). Firstly, we extracted the forecast data from the lead day N−2 to day N+2 of one target day (date M) of the original data array. Secondly, the extracted data are used as the forecast data of N days lead from day M−2 to day M+2 of the new data array, so that the new data are in daily intervals. If the forecast dates in the new data array overlap, the value of the forecast date is treated as the mean of the two values forecast from different initialization dates. This method can maintain consistency between the original data array and the new data array. Note that the new data array is only available after a lead time of 3 days. Thus, only the data from a 3-day lead in ECMWF and CMA are available for comparison.

The observed daily mean surface air temperature at 2 meters altitude (SAT) and precipitation data derived from gauges over China are gridded to a horizontal resolution of 0.25° × 0.25° as provided by the National Meteorological Information Center (CN05.1; Xu et al., 2009; Wu and Gao, 2013). The atmospheric circulation data are from the ERA-Interim reanalysis dataset (Dee et al., 2011), while the OLR (the net outgoing longwave radiation at the top of the atmosphere, OLR) are from the National Oceanic and Atmospheric Administration (Liebmann and Smith, 1996). Forecast SAT, total precipitation, OLR, geopotential height at 500 hPa, mean sea level pressure, and zonal and meridional wind at 850 hPa are derived from the three models. The database consists of 1998−2010 (2005−19) winters from ECMWF and CMA1.0 (CMA2.0) (Vitart et al., 2017). In the present study, the observed anomalies are calculated by subtracting the climatological annual values, while the anomalies in models are calculated by removing their own climatological annual values from the hindcast, computed as a function of both lead time and initialization date.

The PHSE in early 2008 can be divided into four successive phases based on the varying location and strength of PHSE’s subseasonal evolution. Figure 1 shows a time−latitude cross section of the precipitation and SAT anomalies over southern China. It can be seen that the first phase (P1) of the PHSE occurred during 10−14 January, with the precipitation anomalies mainly to the north of 30°N and rapidly decreasing SAT fell below 0°C from 12 January. The region south of 30°N was relatively dry and warm in P1. During the second phase (P2) from 18 to 23 January, the SAT anomaly was around −2°C to −4°C, and the positive precipitation anomaly propagated southwards and was situated around 30°N. From 25−29 January (P3), two positive precipitation anomaly centers appeared over 20°−25°N and 25°−30°N, respectively. The last phase (from 30 January to 3 February, P4) was the coldest with the SAT anomaly below −6°C, and the largest positive precipitation anomaly covered the entire region south of 30°N. The division of this PHSE is consistent with the previous works of Tao and Wei (2008) and Zhou et al. (2009).

Figure 1.  Latitude−time cross section of the SAT anomaly (contours; °C, with a 2°C interval) and precipitation anomaly (shading; mm d⁻¹) averaged over land from 105°−125°E from 7 January to 6 February 2008. The black contour denotes 0°C. The solid black vertical lines divide the four successive phases of the PHSE.

The normalized root-mean-square error (NRMSE) and the pattern correlation coefficient (PCC) are the metrics to assess the models’ performances in forecasting the PHSE. The NRMSE can reflect the skill in capturing the intensity of the predictand (i.e., SAT and precipitation), while the PCC measures the skill in forecasting the spatial pattern of the predictand.

The NRMSE is defined as

and the equation of the PCC is

where N indicates the number of grids in the study area, and $ {R}_{\mathrm{F}\mathrm{I}} $ ($ {R}_{\mathrm{O}\mathrm{I}} $) denotes the value of R in the model forecasting (observation). ${\bar{R}}_{\mathrm{F}}\!=\!(1/N)\displaystyle\sum\nolimits_{I=1}^{N}\!{R}_{\mathrm{F}\mathrm{I}}$ and ${\bar{R}}_{\mathrm{O}}\!=\!(1/N)\displaystyle\sum\nolimits_{I=1}^{N} \!{R}_{\mathrm{O}\mathrm{I}}$ indicate the area-mean of variables. A small NRMSE suggests that the prediction is in good agreement with the observation in intensity, while a high PCC indicates high skill in predicting the spatial pattern of the predictand.

Figure 2 shows the spatial patterns of SAT and precipitation anomalies over southern China during the four individual phases of the PHSE. In P1, the cold SAT anomalies are observed over the middle reaches of the Yangtze River basin, with its center below −3°C, while the southeast coast shows positive SAT anomalies. In P2, the cold SAT anomalies move further to the south and intensify to below −6°C over the southwest part of southern China, while the positive SAT anomalies are greatly weakened and confined to a small region over Fujian Province. In P3 and P4, the cold SAT anomalies are further intensified and dominate the whole of southern China. The cold anomaly center in Guangdong and Guangxi provinces is below −8°C in P4. Note that the cold SAT center gradually moves southward and southwestward with increasing intensity from P1 to P4 (also shown in Fig. 1).

Figure 2.  Spatial distribution of (a) SAT (shading; °C) and (b) precipitation (shading; mm d⁻¹) anomalies over southern China in each of the four phases P1,…P4. The green line indicates the Yangtze River.

Corresponding to the southward propagation and amplification of the cold SAT anomaly, the spatial distribution of precipitation anomalies during the four phases of the PHSE also shows large disparities. P1 and P2 both exhibit a dipole rainfall pattern. In P1, the wet anomalies appear over the middle and lower reaches of the Yangtze River and dry anomalies are observed to the south. In P2, the dipole rainfall pattern shifts southwards, with its wet pole over the lower reaches of the Yangtze River basin and dry pole over much of the area to the south. In P3, the positive rainfall anomaly extends southward and southwestward, which dominates almost the whole of southern China. In P4, the positive rainfall anomaly over the area to the north of the Yangtze River disappears, and the positive rainfall anomaly over the south further intensifies.

Therefore, each phase of this PHSE shows distinct feature, meaning we should not investigate the PHSE’s prediction skill by assessing its time-mean state. Can the S2S models forecast the PHSE’s intraseasonal evolution? Based on the above knowledge of the 2008 PHSE, in the next subsection we assess the skill of the ECMWF and CMA models in predicting the SAT during each of the four phases of the PHSE.

The overall performance of the S2S models in forecasting the PHSE can be evaluated using the NRMSE and PCC with respect to the two key variables (i.e., SAT and precipitation) of the PHSE. Figure 3 shows the NRMSE skills as a function of lead time for each process. In general, the ensemble mean (solid line) displays better skill due to the reduced systematic errors. Thus, we will mainly discuss the results of ensemble mean hereafter. For all the four phases, the NRMSE generally increases with lead time. The NRMSE of SAT (Fig. 3a) are consistently higher than those of precipitation (Fig. 3b) in both CMA models. In ECMWF, the NRMSE of SAT appears promising, with values around 1.0 in P1−P3, but the NRMSE is quite high in P4—especially after a lead time of 14 days (with an NRMSE larger than 2.0). In CMA 2.0, the NRMSE skills for both SAT and precipitation are better than that of CMA1.0 but still lower than that of ECMWF. In general, the NRMSE for SAT of ECMWF is much lower than that of the two CMA models, while it is comparable with CMA for precipitation in each phase.

Figure 3.  NRMSE of (a) SAT and (b) precipitation averaged over the land from 105°−125°E, 20°−35°N in forecasts of the four phases of the PHSE from the three models. The blue, red, orange line denotes the ECMWF, CMA1.0, CMA2.0 respective ensemble mean, and the corresponding shading indicates the spread of the ensemble members for each model.

Figure 4 shows the PCC skill of the PHSE as a function of lead time for each phase, which denotes whether the model can forecast the actual spatial distribution of the predictand. Note that the spread of PCC skill in precipitation is considerable wider than SAT in both models, indicating larger discrepancy among different ensemble members in simulating precipitation pattern. The PCC skill of both CMA models is relatively lower than that of ECMWF in precipitation at all lead times for P1, P2 and P3, and is comparable with ECMWF in P4. The CMA2.0 has a higher PCC than CMA1.0 for SAT, especially for P4. The high PCC skill (larger than 0.5) from a lead time of 3 to 15 days suggests that the three models have promising skill in forecasting the spatial pattern of the SAT (Fig. 4a) and precipitation (Fig. 4b) anomalies in P1, P2 and P3 within a 15-day lead time, although the PCC skill of precipitation forecast by both CMA drops dramatically beyond the 10−15-day period. Note that in P4, the SAT forecast by the three models has very limited PCC skill. Both CMA and ECMWF show no skill (PCC below 0.5) in forecasting SAT even from a very short lead time (Fig. 4a).

Figure 4.  As in Fig. 3 but for PCC.

To show more detailed information on the model performance in forecasting the PHSE, the SAT and precipitation anomaly patterns over southern China in each phase of the PHSE ensemble mean forecasts at lead times of 5, 10, 15 and 20 days are presented in Figs. 5 and 6, respectively. Figure 5 shows the forecasted SAT by ECMWF and CMA. It can be seen that ECMWF and CMA2.0 generally forecast the cold surge during each individual phase. However, in P1, ECMWF underestimates the cold anomalies over the middle reaches of the Yangtze River basin and forecast a much larger area of the warm anomalies compared with the observed one over southeastern coast, especially beyond a 15-day lead. CMA 2.0 shows no skill after 15-day lead. It can also be seen that CMA1.0 fails to forecast the cold SAT even with a 5-day lead, consistent with the high NRMSE in Fig. 3a. None of the three models can forecast the cooling anomaly center over the middle reaches of the Yangtze River after the 15-day lead, resulting in the higher NRMSE (Fig. 3a).

Figure 5.  Spatial distribution of forecast SAT anomaly (shading; °C) over southern China for the four phases of the PHSE (top to bottom) at lead times from 5−20 days (left to right) by ECMWF (upper panels), CMA1.0 (middle panels) and CMA2.0 (lower panels) ensemble mean. The green line indicates the Yangtze River.

Figure 6.  As in Fig. 5 but for the precipitation anomaly (shading; mm d⁻¹).

Compared to the observations, in P2, P3 and P4, all three models underestimate the pronounced negative SAT anomalies over the areas south of the Yangtze River. This is especially the case in P4, which explains the limited PCC skill in this period (Fig. 4a). CMA and ECMWF possess reasonable skill in forecasting SAT anomalies from a lead time of 3 to 20 days in P2 and P3, but unsatisfactory skills in P1 due to missing the cold anomaly over the middle of Yangtze River basin, and in P4 due to the underestimated cold anomaly over the area south of the Yangtze River. In summary, all the models missed the cold anomalies over the middle of Yangtze River basin in P1 beyond the 10-day lead. In P2 to P4, ECMWF always underestimates the cold anomalies south of 30°N. In both of the two CMA models, the forecast cold anomalies center always shifts northwestward, and a warm anomaly even appears over the southeast coast in P2 at the 5- and 20-day lead and in P4 at the 15- and 20-day lead in CMA1.0, resulting in the low PCC skill as shown in Fig. 4a.

Figure 6 displays the forecast precipitation anomalies in each phase. Consistent with the SAT, the pattern of precipitation anomalies in each individual phase can be generally forecast within a 10-day lead. However, beyond the 10-day lead, the forecast skill for precipitation drops dramatically (as shown in Fig. 4b). Specifically, for P1, P3 and P4, ECMWF can basically forecast the pattern of precipitation anomalies within a 10-day lead time, while the dipole rainfall pattern in P2 fails to have been forecast even at the 10-day or longer lead time (also see Fig. 3b). CMA1.0 shows good skill for precipitation with a 10-day lead in P3 and P4 but no skill even beyond a 5-day lead in P1 and P2. In CMA2.0, the wet anomaly is underestimated over middle-lower Yangzi River in P3−P4 at a 10-day lead. In summary, while both models show no skill with respect to precipitation beyond the 15-day lead time, ECMWF clearly outperforms CMA in all phases within a 15-day lead (also see Fig. 4b).

Questions are raised as to why the three models cannot skillfully forecast the SAT and precipitation beyond the 10-day lead time, as well as why they both show very limited skill for SAT in P4. Evaluation of the models’ performances in forecasting the critical circulation systems associated with the PHSE should provide useful clues for further improvement and development of the S2S models. Therefore, the next subsection investigates the capacity of the three models in forecasting these critical circulation systems during the four phases and attempts to explain their forecast skills.

Here we evaluate the ability of the models to forecast the critical circulation systems on the subseasonal time scale to tentatively explore the reasons behind the model errors in SAT and precipitation during the four phases of the PHSE. Firstly, we need to identify the observed critical circulation anomalies that determine the SAT and precipitation anomalies during the PHSE.

Figure 7a shows the observed 500-hPa geopotential height and sea level pressure anomalies over the Northern Hemisphere during the four phases of the PHSE. The 500-hPa geopotential height is characterized by a Siberian High anomaly that persists from P1 to P4 over the mid−high latitudes of the Eurasian continent. The Siberian High is a manifestation of the outbreak of the cold surge, and a precursor of the cold SAT over East Asia (Yang et al., 2018b). Although the Siberian High is weakened from P1 to P4, it moves gradually southward and eastward. Together with the 500-hPa Siberian blocking high anomaly is a positive sea level pressure anomaly which also propagates southeastward, leading to strengthened anomalous northerlies on the eastern side of the surface high-pressure anomaly, which advects cold air into southern China. Therefore, during the southeastward-propagating movement of the Siberian High and the positive sea level pressure anomalies from P1 to P4, southern China becomes colder and colder, and the anomalous cold center shifts gradually from north of the Yangtze River to the south (Figs. 1 and 2a).

Figure 7.  Spatial distribution of (a) 500-hPa geopotential height anomaly (contours; gpm) and sea level pressure anomaly (shading; hPa) and (b) 850-hPa water vapor transport flux (vectors; kg m⁻¹ s⁻¹) and outgoing longwave radiation anomaly (shading; W m⁻²) observed in the four phases during the PHSE. Only the vectors larger than 5 kg m⁻¹ s⁻¹ are shown. The red line outlines the domain of southern China. The OLR positive (negative) anomaly indicates reduced (increased) upward heat flux.

On the other hand, during the whole period of the PHSE, suppressed convection of a strong MJO propagates eastwards from the Indian Ocean to the western Pacific, providing abundant water vapor by way of inducing persistent anomalous southerly flow (Fig. 7b). In P1 and P2, an anticyclonic anomaly is stimulated over the southern part of the Bay of Bengal as a Matsuno−Gill-type Rossby wave response (Gill, 1980) to the northwest of the MJO suppressed convection (diabatic cooling), whilst a cyclonic wind shear is induced over the South China Sea as a Kelvin wave response to the diabatic cooling. The strong southwesterlies to the northwest of the anticyclonic anomaly merge with the northerly flow over the South China Sea and transports the warm, moist air farther into the cold Yangtze River basin, leading to the snow/frozen rain in the Yangtze River basin and relatively dry conditions over the south part of southern China in P1 and P2 (Fig. 2b). Therefore, due to the MJO suppressed convection over the eastern Indian Ocean, the precipitation anomalies present a dipole pattern over southern China. In P3 and P4, the suppressed convection further propagates eastward to the western Pacific, and the Gill-type Rossby wave response in terms of an anticyclonic anomaly is induced (Fig. 7b) to its northwest, merging into the anomalous western North Pacific subtropical High. Therefore, a large anticyclonic anomaly dominates the entirety of the South China Sea and the western North Pacific, and the southwesterly flow along the northwest flank of the anticyclonic anomaly transports abundant water vapor into the cold region of southern China, leading to the snow/frozen rain over most of southern China during P3 and P4 (Fig. 2b).

In summary, during the 2008 PHSE, the characteristics of the Siberian High anomaly in the mid-to-high latitudes determined the intensity and distribution of the SAT over southern China, whereas the features of the MJO controlled the water vapor transport into southern China. Additionally, these two critical circulation anomalies jointly influenced the anomalous western North Pacific subtropical high, and they ultimately determined the intensity and spatial pattern of SAT and precipitation, and eventually the PHSE. Therefore, these two independent and critical circulation systems from both the tropics and mid-to-high latitudes were persistently at play and contributed to the entire intraseasonal evolution of the PHSE.

To track the dynamic sources related to forecast skills of PHSE, the PCC skill for 500-hPa geopotential height over (20°−90°N, 30°−180°E) and OLR over (30°S−30°N, 30°−180°E) were calculated in each phase (Fig. 8). In general, all models show reasonable PCC skills for 500-hPa geopotential height within a 15-day-lead and for OLR within 10-day-lead, while poor skills are found in P4 for CMA1.0. In addition, large spread of PCC skills can be seen between different ensemble members, and the ensemble mean shows superior PCC skill to that of individual member. The PCC skills in forecasting both 500-hPa geopotential height and OLR from the two CMA models are generally lower than that of ECMWF, other than in P2 for 500-hPa geopotential height for CMA1.0. The CMA2.0 possesses a comparable PCC of 500-hPa geopotential height and OLR with CMA1.0, except for the 500-hPa geopotential height in P3.

Figure 8.  As in Fig. 3 but for (a) 500-hPa geopotential height over (20°−90°N, 30°−180°E) and (b) OLR over (30°S−30°N, 30°−180°E).

Figure 9 plots the 5-, 10-, 15- and 20-day-lead 500-hPa forecasted spatial pattern of the 500-hPa geopotential height and sea level pressure anomalies by these three models. At the 5- to 15-day lead time in P1 and P2, the ECMWF and CMA models generally forecast well the Siberian High anomaly at 500-hPa and the corresponding sea level pressure anomalies, which creates the relatively high forecast skill for SAT within a 15-day lead, particularly for ECMWF (Fig. 4a). In P2, the Siberian High anomaly in CMA2.0 is maintained up to a lead time of 20 days, in accordance with higher PCC skill (Fig. 8a) in contrast with the other models. However, in P3, the Siberian High forecast by both CMA models is very weak and shifted northward beyond the 10-day lead, eventually disappearing by the 20-day lead time (15- and 20-day lead) for CMA1.0 (CMA2.0), resulting in the poor skill with respect to SAT (Fig. 4a). In P4, none of the models could reproduce the Siberian High anomalies beyond 15-day lead time in terms of both the intensity and location. The positive sea level pressure anomalies over inland China are greatly underestimated or even missing after only a 5-day lead time, leading to the failure in forecasting the cold SAT anomalies over the south of the Yangtze River in P3 and P4 (Fig. 5). Moreover, because CMA1.0 always overestimates the subtropical western North Pacific high (Fig. 9) as shown in P1 for all the lead times, and P2 and P4 at the 15-day and 20-day lead, the strength of the cold anomalies over southern China cannot be reproduced well. Instead, a warm SAT anomaly is forecast over southern China (see Fig. 5 for P1, P2 and P4 in CMA1.0), which is induced by the overestimated subtropical western North Pacific high because of greater downward solar radiation.

Figure 9.  Spatial distribution of forecast 500-hPa geopotential height (contours; gpm) anomaly and sea level pressure (shading; hPa) anomaly for the four phases of the PHSE (top to bottom) at lead times from 5−20 days (left to right) by ECMWF (upper panels), CMA1.0 (middle panels) and CMA2.0 (lower panels) ensemble mean. The red line outlines the domain of southern China.

Figure 10 plots the 5-, 10-, 15- and 20-day-lead OLR and water vapor flux anomalies at 850-hPa forecast by the three models. In general, all models forecast the signals of eastward propagation of the MJO from P1 to P4. However, the three models gradually fail to forecast the suppressed MJO convection with the increase of lead time, and the bias is obviously much larger in both CMA models than that in ECMWF, which generates considerably lower PCC skill for OLR in CMA (Fig. 8b). Corresponding to the weakened MJO suppressed convection, the southerly water vapor flux is greatly underestimated at the 10- and 15-day lead in P1, P2 and P3 by all models. Note that CMA1.0 forecasts an overly strong western Pacific subtropical high anomaly (Fig. 9), which prevents the water vapor transport into southern China in P1 and P2 at the 10−15-day lead time and in P3 at the 15−20-day lead time, leading to relatively dry conditions in the region. Also note that because both ECMWF and CMA1.0 forecast persistently strong southwestward water vapor flux in P4, they still show skill in forecasting precipitation at least at a 15-day lead (Fig. 4b), regardless of the failure in forecasting the Siberian High and the corresponding cold anomalies over southern China (Fig. 5). For CMA2.0, although the intensity of suppressed convection is stronger than that of CMA1.0 beyond 15-day lead, the suppressed convection is underestimated in the eastern Indian Ocean in P1−P2 and shifts eastward over the western Pacific Ocean in P3−P4 especially at a lead time of 20 days. Moreover, the western Pacific subtropical high forecasted by CMA2.0 is largely underestimated and shifted eastward (Fig. 9) from 15- to 20-day lead in P3 and P4. These factors cause the bias of water vapor flux anomaly (Fig. 10), leading to the dry anomalies over southern China in CMA2.0 in each phase at 15−20-day lead (Fig. 6).

Figure 10.  Spatial distribution of forecast 850-hPa water vapor transport flux (vectors; kg m⁻¹ s⁻¹) and outgoing longwave radiation anomaly (shading; W m⁻²) for the four phases of the PHSE (top to bottom) at lead times from 5−20 days (left to right) by ECMWF (upper panels), CMA1.0 (middle panels) and CMA2.0 (lower panels) ensemble mean. The red line outlines the domain of southern China. Only the vectors larger than 5 kg m⁻¹ s⁻¹ are shown. The red line outlines the domain of southern China.

In summary, all models possess encouraging skill in forecasting the SAT and precipitation during this PHSE over southern China at a 10-day lead time, except for the SAT in P4, which is forecast with very limited skill even at a 5-day lead time. While all models have the capacity to forecast the Siberian High and the corresponding positive sea level pressure anomalies over inland China within a 15-day lead from P1 to P3, they fail to forecast the positive sea level pressure anomalies over inland China in the last phase of the PHSE, leading to the forecast failure for SAT in P4. It is noteworthy that all models poorly forecast the MJO suppressed convection anomalies with increase in lead time, particularly in the CMA models, resulting in the absence of the positive rainfall anomalies in ECMWF from P2 to P4 at the 15- to 20-day lead, in CMA1.0 from P1 to P2 at the 10- to 15-day lead and from P3 to P4 at the 15- to 20-day lead, and in CMA2.0 beyond 15-day lead in every phase.

Given that the PHSE resulted from a combination of cold SAT and positive precipitation anomalies, a successful forecast of this PHSE would require not only an accurate forecast of precipitation (Zheng et al., 2020) but also a correct forecast of SAT. Therefore, based on our evaluation, the skill of ECMWF and CMA in forecasting this PHSE in 2008 is in general only successful within the 10-day lead time.

From the above analysis, the Siberian High and MJO can be considered as the predictable sources of PHSE in 2008. Is the relationship between Siberian High (MJO) and SAT (precipitation) statistically robust from the historical records? To address this question, Singular Value Decomposition (SVD) analyses were conducted between pentad mean 500-hPa geopotential height and SAT anomalies (Fig. 11), and between pentad mean OLR and precipitation anomalies (Fig. 12) in January from 1999−2010. In the first SVD mode, the 500-hPa geopotential height mode mainly features a north-south dipole pattern with positive anomaly centered north of Lake Baikal and negative anomaly over East Asia between 30°−40°N. The corresponding SAT mode is characterized by uniform negative anomalies. These two fields are highly correlated with a correlation coefficient of 0.71, suggesting a robust relationship between the 500-hPa geopotential height over mid-to-high latitudes of Eurasia and SAT pattern over southern China. It can be seen from that the temporal coefficients of geopotential mode and SAT mode both present a peak in late January of 2008, corresponding to the observed Siberian High and cold southern China in P3−P4 in 2008 (Fig. 2a).

Figure 11.  The first two SVD modes between pentad mean (left) 500-hPa geopotential height and (right) SAT in January from 1999−2010. The blue dash (red solid) line denotes the normalized temporal coefficient of the left (right) SVD mode. The black dash lines indicate two standard deviations. The yellow shading masks the period of January 2008. The correlation between the two temporal coefficients is shown in the upper right of the middle panels. The green line in the right panel indicates the Yangtze River.

Figure 12.  As in Fig. 11 but for the first SVD mode between (left) OLR and (right) precipitation.

The 500-hPa geopotential height in the second SVD mode also shows a north-south dipole pattern but with its location shifted westwards about 20°−30°. A zonal geopotential dipole pattern dominates East Asia, with positive anomalies over the East China Sea and negative anomalies over central Asia. The associated SAT anomalies exhibit a northwest-southeast dipole pattern. These two patterns are highly correlated with a correlation coefficient 0.76. The coupled fields reoccurred in P1−P2 in 2008 in the observations as shown in Fig. 7a and reflected in the high value of their temporal coefficients in early January of 2008 (Fig. 11b). Note that the temporal coefficients exceed two standard deviations in both two SVD modes in January 2008 (Fig. 11), indicating how extreme is the PHSE.

The SVD mode for OLR and precipitation shows the relationship between a tropical zonal dipole convection pattern and enhanced precipitation over all of southern China (Fig. 12). The temporal coefficient of this coupled mode shows a peak in late January of 2008, consistent with the observations in P3−P4 as shown in Figs. 2b and 7b.

The SVD analyses implied that the two critical systems, the persistent Siberian High and the propagated MJO, are the key origins for the intraseasonal variations of SAT and precipitation over southern China in January. The abnormal Siberian High and MJO in early 2008 are the direct reasons for the PHSE. Here, from the SVD analyses, we conclude that these two predictability sources for January climate over southern China are not only suitable for one single case (e.g., the PHSE in 2008), but also shed light on subseasonal prediction of normal condition and other extreme events over southern China.

The present study evaluated the subseasonal prediction skills of the ECMWF and two CMA (CMA1.0 and CMA2.0) operational models from the S2S database for the four successive phases (from P1 to P4) of the PHSE in early 2008 over southern China. In general, the three models possess promising skill in forecasting the spatial pattern of the SAT and precipitation anomalies in P1, P2 and P3 within a 15-day lead, although the PCC skill of the precipitation forecasted by CMA drops dramatically after a 10-day lead. However, all three models failed to capture the SAT in P4. The CMA1.0 (CMA2.0) possessed no PCC skill with respect to SAT after a 9-day lead (12-day lead), while ECMWF showed no PCC skill even from the very beginning.

Previous studies have investigated the underlying mechanisms for the PHSE, and in doing so emphasized factors such as the Eurasian westerly jet stream (Wen et al., 2009), the autumn Arctic sea ice (Wu et al., 2011; Li and Wu, 2012) and diabatic heating of the Tibetan Plateau (Bao et al., 2010) in influencing the mean state of the PHSE. Given the strong subseasonal variation of this PHSE in 2008, here we divided the whole PHSE into four individual phase and identified that the persistent Siberian High over the mid-to-high latitudes Eurasia and the eastward propagated MJO from the tropical Indian Ocean to western Pacific were the two critical circulation systems for the whole PHSE.

The reason behind the forecast skill for this PHSE lies in the performance in forecasting these two critical factors. While all models have the capacity to forecast the Siberian High and the corresponding inland positive sea level pressure within a 15-day lead from P1 to P3, they fail to forecast the positive sea level pressure in that last phase of the PHSE (P4), leading to the low forecast skill for SAT in P4. It is noted that all models underestimated the intensity and misrepresented the location of the suppressed convection anomalies of MJO as the lead time increased, leading to the absence of precipitation from P1 to P3 at the 15- to 20-day lead time. Given that the PHSE resulted from the combination of cold SAT and positive precipitation anomalies, a successful forecast of this PHSE would require an accurate forecast of the evolution of Siberian High and MJO eastward propagation. Therefore, based on our evaluation, the skill of ECMWF and CMA in forecasting this PHSE in 2008 is in general only useful within the 10-day lead time.

Whereas MJO is a planetary scale intraseasonal phenomenon, the forecasted OLR anomalies as shown in Figs. 8b and 10 may not only involve MJO but also other tropical waves (e.g., the Kelvin waves). One may doubt the inaccuracy in evaluating the MJO activity simply based on the OLR anomalies over a rectangle domain in the Indian Ocean and western Pacific (30°S−30°N, 30°−180°E). To address this issue, the OLR-based MJO index (OMI) developed by Kiladis et al. (2014) was employed to represent the MJO activity. This OMI index has better capacity for tracking the seasonality of the propagating features of both the MJO and Boreal Summer Intraseasonal Oscillation (Wang et al., 2019). We calculated the bivariate correlation skills of OMI predicted by three models from Jan. 6 to Feb. 6 in 2008. As shown in Fig. 13, both CMA1.0 and CMA2.0 models show reasonable skills (relatively higher than ECMWF) in predicting OMI within a 13-day lead time, whilst the ECMWF model can capture the OMI up to a 20-day lead. The results are similar to that measured by PCC skills as shown in Fig. 8b. Moreover, based on the perfect prediction assumption in OMI (first defined by Wang et al., 2019), we can determine whether the forecasting errors in forecasting MJO activity is from the bias in forecasting its amplitude (COR_a) or its propagation phase (COR_p). It can be seen that the two CMA models have much better capability in predicting the amplitude than the phase of MJO after a week’s lead. While the three models all show comparably promising skills in forecasting MJO amplitude (higher than 0.8) (Fig. 13b), the ECMWF maintained a good skill in forecasting the MJO propagation after 10-day lead, which is the main reason for the overall outperformance in forecasting MJO activity up to 20 days (Figs. 13a and 8b). Note that we used band-pass filtered OMI index here, which may slightly underestimate the MJO prediction skill compared with real-time OMI index (Kiladis et al., 2014; Wang et al., 2019), but this does not influence the results much.

Figure 13.  Bivariate correlation skills based on (a) OMI (COR), (b) OMI with perfect phase (COR_a) and (c) OMI with perfect amplitude (COR_p) in three models from Jan. 6 to Feb. 6, 2008. The blue, red, orange line denotes the ECMWF, CMA1.0, CMA2.0 ensemble mean, and the corresponding shading indicates the spread of the ensemble members for each model. Assuming the phase (amplitude) is perfectly forecast, the COR_a (COR_p) is determined completely by the amplitude (phase) correlation.

In the present study, we conclude that the forecast failure at longer lead times of the PHSE is mainly attributed to the fact that the models cannot reproduce the persistence of the Siberian High and accurately capture the MJO propagation. However, one may argue that the forecast failure could be model dependent due to different reasons (e.g., initial conditions, data assimilation schemes). To support our claim that the capacity in forecasting this PHSE event mainly relies on the models’ skills in forecasting critical circulation systems (i.e., Siberian High and MJO), we calculated the linear correlation between the prediction skills (PCC) of Siberian High (MJO) and SAT (precipitation) over southern China. As clearly shown in Fig. 14, in EC and CMA 2.0 ensemble members, the PCC skill for 500-hPa geopotential height (representing the Siberian High) is significantly correlated with that of SAT over southern China, whilst the PCC skill for OLR (denoting the MJO) is highly correlated with that of precipitation over southern China, indicating that the models’ capability in predicting Siberian High and MJO is indeed the key for forecasting skill of PHSE. However, it is also noted that the correlation between prediction skills of 500-hPa geopotential height (OLR) and SAT (precipitation) is quite low for CMA1.0, suggesting that the forecast error may also relate to the data assimilation scheme. Therefore, both data assimilation schemes (given that the CMA1.0 and CMA2.0 have totally different data assimilation schemes) and capacity in capturing dynamic sources (two critical circulation systems) are responsible for the forecast error of PHSE. Although both CMA models show lower skills than that of ECWMF in forecasting PHSE, it is encouraging that the CMA2.0 is on the right track considering its ability in capturing the observed relationship between the critical circulation systems and regional meteorological variables (Fig. 14).

Figure 14.  Scatter diagram for PCC skills between (a) HGT (20°−90°N, 30°−180°E) at 500 hPa and T2m (southern China) anomalies, and between (b) OLR (30°S−30°N, 30°−180°E) and precipitation (southern China) anomalies from three models. The black, blue, red, and orange curves are linear fits to the total (1368 dots), ECMWF (792 dots), CMA1.0 (288 dots), CMA2.0 (288 dots). Correlation coefficients between PCC skills are shown in the left bottom of the panel in corresponding colors. The “*” (“**”) indicates the correlation is significant at 90% (95%) confidence level.

It is noteworthy that the snow cover or soil moisture are also potential factors of the subseasonal prediction of the PHSE via influencing the land- atmosphere coupling and model’s initialization. The feedback between snow cover and atmosphere may maintain the extreme event (Tao and Wei, 2008; Zhou et al., 2009; Xiang et al., 2020), and the initial condition of snow cover and the corresponding land processes (Jeong et al., 2012; Orsolini et al., 2013; Li et al., 2019) can also interfere with the subseasonal prediction skills. Thus, the effect of land processes on subseasonal prediction of PHSE or such extreme events deserves further investigation.

The SVD analyses revealed a good relationship between two dynamic sources and PHSE. Therefore, an empirical model using large-scale circulation system (predictors, e.g., Siberian High, MJO) can be established to predict regional climates (predictands, e.g., SAT, precipitation). Given that the large-scale circulation generally presents better prediction skills than regional meteorological variables (e.g., SAT and precipitation) in S2S model (Vitart, 2017; also see Fig. 4 and Fig. 8 in this article), especially beyond 10-day lead, the circulation predictors forecast by S2S models can be introduced to make empirical subseasonal prediction. By this dynamical-statistical combined method, the subseasonal prediction skill for SAT and precipitation over southern China is expected to be improved. Currently, we are developing the dynamical-statistical hybrid models for subseasonal prediction of extreme events in China. The results will soon be reported elsewhere.

Acknowledgements. The authors greatly appreciate the professional and earnest review made by the anonymous reviewers which for sure improved the quality of our manuscript. This work was supported by the National Key R&D Program of China (Grant Nos. 2018YFC1505905 & 2018YFC1505803), the National Natural Science Foundation of China (Grant Nos. 42088101, 41805048 and 41875069). Tim LI was supported by NSF AGS-1643297 and NOAA Grant NA18OAR4310298.