HTML

The daily NCEP/NCAR Reanalysis1 geopotential height and wind field data (Kalnay et al., 1996) and the NOAA outgoing longwave radiation (OLR) data (Liebmann and Smith, 1996) are used for the analysis and model verification with a horizontal resolution of 2.5° × 2.5°. The vertical integrated moisture flux is calculated by using the ECMWF reanalysis (ERA)Interim data (Dee et al., 2011) with a horizontal resolution of 1.5° × 1.5°. A daily highresolution (0.25° × 0.25°) gridded rainfall dataset (CN05.1; Wu and Gao, 2013) is used to avoid the uneven spatial distribution of station data. All the observational data cover the time period of 1981–2018 and the climatology is defined as the 30yr average during 1981–2010.

Two dynamic prediction systems—namely, the BCC Climate Prediction System for S2S version 1 (BCCCPSS2Sv1) based on BCC_CSM1.2 and the Global Seasonal Forecasting System (GloSea5) from the UKMO—are used to analyze the representation and predictability of the EAP teleconnection on the subseasonal timescale.
BCC_CSM1.2 is a fully coupled global climate–carbon model (Wu et al., 2014). It consists of the atmospheric component of the BCC Atmospheric General Circulation Model version 2.3 (BCC_AGCM2.3) with a triangular 106 horizontal resolution and 40 hybrid sigma/pressure vertical levels (T106L40), ocean component of the Modular Ocean Model version 4 (MOM4L40), land component of the BCC Atmosphere and Vegetation Interaction Model version 1.0 (BCC_AVIM1.0), and sea ice component of the sea ice simulator. The initial conditions are generated by applying a nudging technique toward the NCEP1 dataset for the atmospheric analysis and threedimensional variational data assimilation method for the oceanic analysis (Jie et al., 2014; Liu et al., 2017). BCC_CSM1.2 is run every day to make a 60day integration for the hindcast period of 1994–2014. Each forecast consists of 4 laggedaverage forecast ensemble members, which are initialized at 0000 UTC on the first forecast day and at 1800, 1200, and 0600 UTC on the previous day.
The GloSea5 model (MacLachlan et al., 2015) consists of the atmospheric component of the Global Atmosphere 6.0 (GA6.0) at N216 and 85 vertical levels with the top at 85 km, ocean component of the Global Ocean 5.0 (GO5.0; Megann et al., 2014) of the Nucleus for European Modelling of the Ocean (NEMO), sea ice component of the Global Sea Ice 6.0 (GSI6.0; Rae et al., 2015), and land component of the Global Land 6.0 (GL6.0; Walters et al., 2017) of the Joint UK Land Environment Simulator (JULES). The hindcasts are initialized on the 1st, 9th, 17th, and 25th of each month for 1993–2015, with the 60day integration and 7 ensemble members based on the laggedaverage forecast method.
This study mainly focuses on the hindcasts initialized during May–September (MJJAS) for the time periods of 1995–2014 in BCC and 1993–2015 in the UKMO subseasonal prediction models. All the hindcast data are horizontally interpolated onto 2.5° × 2.5° grids.

A nontraditional filtering method suitable for realtime prediction similar to that used by Wheeler and Hendon (2004) and Hsu et al. (2015) is applied to observational and model hindcast data to extract the subseasonal signal of the EAP teleconnection. First, the climatological annual cycle is removed by subtracting the 0–3 harmonics of the climatological daily data; second, the interannual variability is removed by subtracting the previous 120day running mean field; and third, the synopticscale variability is removed by applying a 5day running average. Note that the model’s climatology is calculated for each lead day (LD) and the previous 120day mean as the forecast on Day j consists of the preceding 120 − (j + 1) days of observational data and the forecast data from Day 1 to Day (j − 1).
Similar to MJO, the EAP prediction skill is measured by using the correlation coefficient (COR) and rootmeansquare error (RMSE; Kim et al., 2018), which are defined as below:
$$ {\rm{COR}}(\tau) = \dfrac{{\displaystyle\sum\limits_{t = 1}^N {[{a}(t){b}(t,\tau)]} }}{{\sqrt {\displaystyle\sum\limits_{t = 1}^N {[a^2(t)} ]} \sqrt {\displaystyle\sum\limits_{t = 1}^N {[b^2(t,\tau)} ]} }}, $$ (1) $$ {\rm{RMSE}}(\tau) = \sqrt {\frac{1}{N}{{\sum\limits_{t = 1}^N {[{a}(t)  {b}(t,\tau)]} }^2}}, $$ (2) where
$ {a}(t) $ is the observed EAP index (I_{EAP}) on time t,$ {b}(t,\tau) $ is the corresponding forecast for time t with a lead time of$\tau $ days, and N is the number of forecasts.${\rm{COR}}(\tau) = 0.5$ is usually chosen as a threshold for skillful prediction because the large sets of hindcasts will dramatically decrease the significant correlation coefficient criterion (Rashid et al., 2011). Because the EAP index is normalized, RMSE for a climatological forecast (I_{EAP} = 0) is 1.0, which will be taken as the threshold for skillful prediction of RMSE (Lin et al., 2008).This study uses two different approaches to evaluate the EAP predictability. The first approach is based on anomaly correlation metrics (Kim et al., 2014). The predictability of a single member is calculated by the
${\rm{COR}}(\tau)$ between one ensemble member and each of the other ensemble members, and then averaged over subsamples. The ensemble predictability is measured in the same way, but between one ensemble member and the mean of the remaining ensemble members. In general, the predictability of a single member is lower than that of the ensemble mean.The second predictability evaluation method is defined by the signaltonoise metrics (Neena et al., 2014). The EAP signal (S_{j}) on LDj is calculated as below:
$$\left\langle {S_j^2} \right\rangle = \frac{1}{{N \times nk}}\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^{\rm nk} {[{{({\rm EAP}_{ij}^k)}^2}]} },$$ (3) where N is the total number of forecasts, nk is the ensemble size, and
${\rm EAP}_{ij}^k$ represents the forecast EAP index of a particular case i on LDj by a single ensemble member k. The meansquare error (E_{j}) for each LDj is given by:$$\left\langle {E_j^2} \right\rangle = \frac{1}{{N \times {m_1}}}\sum\limits_{i = 1}^N {\sum\limits_{\rm pair=1}^{m_1} {[{{({\rm EAP}_{ij}^{{{k}_1}}  {\rm EAP}_{ij}^{{k_2}})_{\rm pair}^2}}]} } .$$ (4) Here, k_{1} and k_{2} represent the control and perturbed forecasts and m_{1} represents the number of possible control–perturbed pairs. Although a single ensemble member is taken as the control forecast, the perturbed forecast is different for the single member or the ensemblemean predictability estimate. In the former, every other ensemble member is regarded as the perturbed forecast; in the latter, it takes the average of all the other ensemble members. Therefore the predictability is defined as the LD at which the meansquare error becomes as large as the mean signal.
Note that, in some studies (e.g., Scaife et al., 2014), the actual skill that was evaluated by validating predictions against the corresponding observations is referred to as “predictability” and the potential skill based on treating one of the ensemble members as a proxy for the observations is referred as “the perfect predictability”. These concepts are identical in meaning to the “prediction skill” and “predictability” in this paper.
2.1. Observations
2.2. Model hindcasts
2.3. Methodology

In order to extract the spatial pattern of the EAP teleconnection, an empirical orthogonal function (EOF) analysis was applied to the normalized 5day running mean 500hPa geopotential height (H500) anomalies in summer (June–August; JJA) during the time period of 1981–2018 over the region (10°–75°N, 110°–150°E). Following the method described by Hart and Grumm (2001), the standard deviation used for standardization was calculated via a 21day smooth window (from 1981 to 2010) centered on the day being investigated. We took the seasonal variation of the variability of geopotential height into consideration. Figure 1a shows that the leading EOF mode with the explained variance of 18.6% is characterized by a clearly meridional tripole structure with two positive anomaly centers located in the tropical western Pacific (WP) and northern Sea of Okhotsk (OK), and a negative center located in the midlatitudes of East Asia (EA), which is considered to be the canonical positive phase of the EAP teleconnection (Huang and Li, 1987; Wang et al., 2016). This suggests that the EAP teleconnection is also the dominant teleconnection mode over East Asia on the subseasonal timescale during boreal summer.
Figure 1. (a) The first eigenvector of the normalized 5day running mean 500hPa geopotential height (H500) anomalies in JJA 1981–2018, the value of which is multiplied by the square of eigenvalue to represent correlation coefficients (CORs). The geopotential height anomalies in the boxed areas are chosen to define the EAP index (I_{EAP}). (b) The normalized time series of the first eigenvector (PC1; blue lines) and I_{EAP} (red lines). (c) The mean seasonal cycle of the standard deviation of I_{EAP} obtained via a 21day smooth running window during the time period of 1981–2010. VAR denotes variance in (a) and ACC denotes anomaly correlation coefficient in (b).
A daily EAP index (I_{EAP}) is defined according to these three atmospheric activity centers over WP (12.5°–22.5°N, 115°–135°E), EA (35°–45°N, 120°–140°E), and OK (60°–70°N, 120°–140°E).
$${I_{{\rm{EAP}}}} = ({H_{{\rm{WP}}}}  {H_{{\rm{EA}}}} + {H_{{\rm{OK}}}})/3 ,$$ (5) where H_{WP}, H_{EA}, and H_{OK} represent the geopotential height anomalies averaged over WP, EA, and OK, respectively. Note that H_{WP}, H_{EA}, and H_{OK} were normalized by the standard deviation calculated via a 21day smooth running window. The I_{EAP} can be calculated for all seasons, but should be normalized by the standard deviation of JJA during the time period of 1981–2010. The time series of I_{EAP} is highly consistent with that of the first principal component (PC1) during JJA from 1981 to 2018 (Fig. 1b), with a correlation coefficient of 0.93. It confirms that the EAP index defined here is capable of capturing the variability of leading EOF mode. Figure 1c shows the seasonal cycle of the amplitude of I_{EAP} obtained by calculating the standard deviation of I_{EAP} via a 21day smooth running window. The results show that the amplitude of I_{EAP} is generally small during boreal winter and increases around late April, reaching a maximum from late June to August. This unique seasonal cycle implies that the EAP relationship among WP, EA, and OK characterizes JJA in boreal summer. Although the standard deviation of I_{EAP} is also relatively high in late October, this month is usually considered as the transition season, with a climatology distinct from that of boreal summer. Therefore, consistent with previous studies (such as Wang et al., 2016), boreal summer is defined here as extending from May to September (MJJAS).
The EAP indices from 1981 to 2018 are analyzed yearbyyear by wavelet transform and then averaged for 38 yr (Fig. 2a), and the MJJAS daily average spectrum is shown in Fig. 2b. The shading and red dotted line (values > 0) represent the denoised energy spectrum value at the 95% confidence level, corresponding to the significant oscillation period. The largest wavelet power spectrum is in midAugust and has a period of about 65 days. Two significant bands are seen in the oscillation period of the EAP index with a period of 50–70 days throughout MJJAS and a period of 10–30 days from midMay to July and from midAugust to September, peaking at around 60 and 25 days, respectively.
Figure 2. (a) The averaged wavelet power spectrum (contours) of EAP indices during May–September (MJJAS) of 1981–2018. (b) The averaged wavelet spectrum during MJJAS. The shading in (a) and the red dotted line (values > 0) in (b) represent the 95% confidence level for the wavelet spectrum.
Figure 3 shows the lead–lag relations between the EAP teleconnection and some wellrecognized subseasonal dominant modes (such as MJO and BSISO). There is only a weak relationship between the EAP teleconnection and MJO in terms of RMM indices, suggesting a limited influence of MJO on EASM (Lee et al., 2013). By contrast, the EAP teleconnection can, to some extent, be reflected by BSISO, with a maximum simultaneous correlation of 0.38 with BSISO11 (PC1 in Lee et al., 2013) and a maximum correlation of −0.34 when the EAP teleconnection leads BSISO21 (PC3) by 1 day. This result is reasonable because BSISO11 and BSISO21 can also capture the cyclone or anticyclone over the western Pacific. However, compared with the BSISO indices defined by Lee et al (2013), which mainly cover the areas of Indian summer monsoon and western north Pacific summer monsoon, the EAP index focuses on East Asia and reflects the anomalies over middle and high latitudes. In addition, the relatively higher correlations were obtained among the EAP index, East Asia–western North Pacific ISO index (Lin, 2013), and the western Pacific ISO index (Lee and Wang, 2016; Yang et al., 2020), showing that these indices also can reflect the anomalies over mid and highlatitude East Asia (figures omitted).
Figure 3. Lead–lag correlation coefficients between (a) the EAP index with RMM1, RMM2, and itself; (b) the EAP index with BSISO11, BSISO12, BSISO21, and BSISO22; and (c) the EAP index and its time tendency, during MJJAS of 1981–2018. The dashed line is the 95% confidence level for the correlations in (a) and (b).
As suggested by Plaut and Vautard (1994), the time tendency of I_{EAP} (denote as I_{EAP}′) has a high lead–lag correlation with I_{EAP}, which can be calculated by the central finite difference method:
$${I_{{\rm{EAP}}}}^\prime (t) = [{I_{{\rm{EAP}}}}(t + \Delta t)  {I_{{\rm{EAP}}}}(t  \Delta t)]/2\Delta t,$$ (6) where
$\Delta t$ represents the length of time and is equal to one day in this study. The analysis of lead–lag correlations between the EAP index and its tendency (Fig. 3c) shows that I_{EAP}′ tends to lead I_{EAP} by about 4 days for variability in the range of 10–30 days with a correlation of 0.89 and at a lag of 10 days for variability in the 30–90day range with a correlation of 0.91.Given the strong lead–lag behavior of I_{EAP} and I_{EAP}′ (Figs. 4a, b), it is reasonable to establish a twodimensional phase–space diagram of [−I_{EAP}′, −I_{EAP}] to monitor the evolution of EAP events visually, with the state of EAP teleconnection as a point in the diagram. Similar to Wheeler and Hendon (2004), we divide the twodimensional phase–space diagram into eight phases (Fig. 4c) and the implication of each phase can be inferred by following the definition and combined position of I_{EAP} and I_{EAP}′. To be specific, Phases 2–3 (6–7), in which amplitudes of the normalized EAP index are larger than their corresponding normalized time tendency (I_{EAP} > I_{EAP}′), indicate the positive (negative) peak phases of EAP events with a “+ − +” (“− + −”) meridional wave train with the anomalous center of precipitation located over EA (WP and OK). The difference between Phases 2 and 3 (or Phases 6 and 7) lies in the sign of I_{EAP}′, which represents the EAP events in the developing (Phases 2 and 6) or decaying (Phases 3 and 7) state. By contrast, amplitudes of the normalized time tendency are larger than their corresponding normalized EAP index (I_{EAP} < I_{EAP}′) in Phases 8–1 and 4–5, which can be viewed as transition phases from positive to negative events (Phases 4–5) or negative to positive events (Phases 8–1). Figure 4c shows the evolution of EAP activities during MJJAS 2016 in terms of a phase–space diagram. Many of the sequential days trace anticlockwise circles around the origin as a result of the lead–lag relationship between I_{EAP} and I_{EAP}′, with two strong positive and three strong negative EAP events signified by large amplitudes.
Figure 4. Time series of the normalized 5day running mean EAP index and its time tendency in MJJAS of (a) 2016 and (b) 2018. (c) The phase–space diagram of [−I_{EAP}′, −I_{EAP}] in MJJAS 2016, with different months in different colors.
Considering the dual significant spectra of the EAP index (Fig. 2) rather than the lead–lag composite, phase composites for anomalies in H500 and OLR are shown in Fig. 5 to depict characteristics of the evolution of EAP events. The composites are performed only on those cases where the amplitude of [−I_{EAP}′, −I_{EAP}] is > 1.0—that is, (I_{EAP}^{2} + I_{EAP}′^{2})^{1/2} > 1.0. The wave activity flux (WAF) defined by Takaya and Nakamura (2001) is used to diagnose the energy dispersion of Rossby wavelike perturbations, which is independent of the wave phase and parallel to the local group velocity. In the transition phase to positive EAP events (Phase 1), a significant negative anomaly center exists over the highlatitude area of Europe centered near Novaya Zemlya, from which strong Rossby wave activity propagates east and southeast to the OK and EA. The poleward energy dispersion from the subtropical WP is discernible, corresponding to the suppression of convective activity to the east of the Philippines. During the peak phases (Phases 2–3), anomalies over the OK are enhanced by obtaining the upstream energy from Novaya Zemlya, which is weakened itself.
Figure 5. The composite of anomalous OLR (color shading; W m^{−2}), H500 (contour; gpm), and WAF (vector; m^{2} s^{−2}) in the phase–space diagram of [−I_{EAP}′, −I_{EAP}]. The composited OLR anomalies significant at the 95% confidence level and geopotential height anomalies > 5 gpm are shown. The number of days falling within each phase category is given on the top right of each panel.
The depressed convection in WP propagates westward, and the convergence between the poleward wave fluxes from WP and southward wave fluxes from the OK markedly strengthen the negative center over EA. Therefore, the canonical EAPtype meridional tripole formed with a zonally extended structure. In Phase 4, the subtropical depressed convection continues to move westward, and the Rossby energy disperses east and southeast from the OK and EA center, enhancing the anomalies over North Pacific and the Aleutian Islands. New negative anomalies and active convection emerge over Novaya Zemlya and WP, which could be the precursor of the next negative EAP event. The spatial evolution of negative EAP events during Phases 5–8 is nearly the same as the opposite signs of positive events, except for the stronger amplitude of subtropical convection anomalies as well as the slightly weaker OK anomaly center, indicating the greater contribution of subtropical activities to the formation of negative EAP events. In addition, according to the numbers of days falling within each phase shown in Fig. 5, the observed EAP states are basically even distributed, with the maximum discrepancy of 25 days between Phases 1 and 7.
Figure 6 shows the phase composites of the vertically integrated moisture flux and precipitation anomalies in China. In general, the spatial pattern of the moisture flux is strongly modulated in EA, as are the rainfall anomalies over North China and south of the YRV basin during different stages of the life cycle of the EAP teleconnection. In Phase 1 (the transition phase of positive EAP events), an anomalous moisture flux cyclone lies in the Sea of Japan, steering moisture to Northeast China on its northern flank and enhancing local precipitation. In Phases 2–3 (the peak phase of positive EAP events), the anomalous anticyclone over the subtropical western Pacific migrates westward, forming an anticyclone–cyclone pair at low to midlatitudes of EA. The confluence between the northerly winds driven by the anomalous cyclone as well as the southwesterly winds to the northern flank of the anticyclone significantly strengthen precipitation over the YRV basin, consistent with previous studies on synoptic scales (Chen and Zhai, 2015; Wang et al., 2018). The anomalous northerlies over North China weaken the transport of moisture and local precipitation. The lowlatitude cyclone moves westward and strengthens the southwesterly winds over the south of YRV in Phase 4. At the same time, the transport of moisture originating from the Bay of Bengal is reinforced, leading to anomalous rainfall over the eastern Tibetan Plateau and southern YRV. The configurations of negative EAP events during Phases 5–8 also nearly mirror the positive events, except for the stronger amplitude of the anomalous lowlatitude cyclone, which is consistent with the results shown in Fig. 5.
Figure 6. As in Fig. 5, but for the anomalous precipitation rate in China (color shading; mm day^{−1}) and vertically integrated water vapor flux (vector; 10^{4} g m^{−1} s^{−1}). Stipples for precipitation and thick black vectors for water vapor flux mark the results that are significant at the 95% confidence level.