HTML
--> --> -->Using near-equatorial station data, Madden and Julian (1971, 1972) discovered that the ISO is a low-frequency signal of anomalous convection and circulation, which propagates eastward along the equator. Follow-up studies (Madden and Julian, 1994; Zhang, 2005) reported that the ISO has strong seasonality: though the eastward-propagating ISO exists all year round, it is the strongest during the boreal winter. In contrast, the eastward-propagating ISO at the equator is much weaker during the boreal summer. Instead, a northward-propagating ISO appears in the Asian summer monsoon region, which influences the onset and active-break periods of the Indian monsoon (Yasunari, 1980; Webster et al., 1998; Annamalai and Slingo, 2001), and causes severe weather events (such as heavy precipitation, typhoons and heat waves) in East Asia and northwestern Pacific monsoon regions (Maloney and Dickinson, 2003; Mao and Wu, 2006; Hsu et al., 2011, Hsu et al., 2016, Hsu et al., 2017; Li et al., 2015). To distinguish this northward-propagating ISO from the equatorial eastward-propagating ISO, which is often called the Madden-Julian Oscillation (MJO), the ISO in the Asian summer monsoon region is referred to as the boreal summer intraseasonal oscillation (BSISO; Lee et al., 2013).
The main predictability of extended-range or subseasonal forecasts is via the ISO (Waliser, 2006). Nowadays, both statistical models and dynamic models can predict the MJO and BSISO to some extent. It was found that some statistical models can successfully forecast the MJO as far as 20-30 days in advance, when using the real-time multivariate MJO index defined by (Wheeler and Hendon, 2004) (Seo et al., 2009; Kang and Kim, 2010; Cavanaugh et al., 2015; Zhu et al., 2015). With the improvement of model physics in dynamic models and the development of data assimilation techniques and computing resources, atmospheric and coupled ocean-atmosphere general circulation models can forecast the MJO two to three weeks in advance (Lin et al., 2008a; Rashid et al., 2011; Wang et al., 2014). In recent years, the capability of predicting the MJO by using the prediction model of the European Centre for Medium-Range Weather Forecasts (ECMWF) has been improved year after year; it can now effectively predict the MJO 27 days in advance. (Xiang et al., 2015) used the coupled atmosphere and ocean model of the Geophysical Fluid Dynamics Laboratory (GFDL) to study its MJO forecast capability. The authors indicated that the GFDL coupled model achieved similar success of predicting the MJO 27 days in advance as the ECMWF model. In comparison, model forecast skills for the BSISO are relatively weak. (Fu et al., 2013) evaluated the skills of predicting BSISO precipitation in the summer monsoon regions by using the models of the ECMWF and the National Centers for Environmental Prediction (NCEP), as well as by using the coupled model of the University of Hawaii. They found the pattern correlation coefficient (PCC) of model results and observation dropped significantly to below 0.5 after one to two weeks. (Lee et al., 2015) analyzed six forecast models that participated in the Intraseasonal Variability Hindcast Experiment. They reported that the forecast skill of the real-time BSISO index (Lee et al., 2013) by a single model was only 5-15 days, and that by an ensemble could be as long as 10-30 days.
To improve the forecast skill of the ISO and to fill the gap between medium-range and seasonal forecasts (Brunet et al., 2010), the World Weather Research Programme and the World Climate Research Programme, under the World Meteorological Organization, jointly organized a project for subseasonal-to-seasonal (S2S) forecasting (Vitart et al., 2017). One of the tasks of the S2S project is to evaluate and diagnose the sources of errors in ISO simulations by surveying the ISO forecast capability of operational models. The goal is to improve S2S forecast skills by numerical models. At present, 11 institutions provide their forecast models' outputs, including the China Meteorological Administration (CMA), the ECMWF, the Japan Meteorological Agency (JMA), the NCEP of the U.S., the Australian Bureau of Meteorology (BoM), the UK Met Office (UKMO), the Centre National de Recherche Meteorologiques (CNRM) of France, the Institute of Atmospheric Sciences and Climate of the National Research Council of Italy, the Hydrometeorological Centre of Russia, the Korea Meteorological Administration, and Environment and Climate Change Canada (ECCC). These S2S model outputs are available for downloading from the website (http://s2s.cma.cn/dataset). A recent paper by (Jie et al., 2017) assessed the forecast skill of BSISO using an S2S database that contains the hindcast outputs from 10 operational models. Their results showed that the skill of a single model forecast is around 6-17.5 days. Some models (such as ECMWF, UKMO, BoM, CNRM, and ECCC) show a remarkable increase in their forecast skills (up to 10-24.5 days) when using the ensemble mean. However, only a slight improvement was found for the ensemble mean of the NCEP, CMA and JMA models (Jie et al., 2017).
The S2S prediction model of the CMA is based on version 1.1 of the Climate System Model of the Beijing Climate Center (BCC_CMS1.1), which participates in phase 5 of the Coupled Model Intercomparison Project. Following the hindcast experiment design of the S2S project, BCC_CSM1.1 was used to conduct daily forecasts for the future 60 days during 1994-2014. In this study, we refer to this hindcast dataset as the prediction results of the BCC S2S model. (Liu et al., 2017) investigated the forecast skills of the MJO during winter by the BCC S2S model. They noted that the model could provide skillful forecasts for up to 16 days in advance; in other words, the model could predict the winter MJO 16 days in advance with accuracy in terms of spatial pattern, temporal evolution and propagation. The prediction skill of monsoon ISO by the BCC S2S model can be found in (Jie et al., 2017), who evaluated all 10 operational models in the S2S project. As it will also be shown later in this study, the BCC S2S model can only predict the 30-90-day northward propagating BSISO about 10 days in advance. This skill is relatively low among the operational models. To improve its forecast capability, it is important to identify the processes that inhibit the northward propagation in the BCC S2S model at longer forecast leads. This will help us understand the key factors governing the BSISO propagation in some models.
The paper is organized as follows. In section 2, we introduce the BCC S2S model, its outputs, and the methods used to evaluate the model capability in BSISO prediction and to diagnose the model errors in capturing the northward-propagating BSISO. In section 3, the forecast skills of BSISO strength, spatial modes and propagation are analyzed. In section 4, the main mechanisms that promote the northward propagation of BSISO are examined at different forecast lead times to understand what caused the degraded BSISO forecast skill at longer lead times. A summary and discussion are given in section 5.
2.1. BCC S2S model
The BCC S2S model is based on BCC_CSM1.1, with the inclusion of atmosphere-ocean-ice-land coupling (Wu et al., 2014). The atmospheric component of BCC_CSM1.1 is the BCC Atmospheric General Model, version 2.1 (Wu et al., 2010), with a horizontal resolution of T106 and 40 layers in the vertical direction. The land component is the BCC Atmosphere and Vegetable Interaction Model, version 1 (Ji, 1995). The oceanic component is the GFDL Modular Ocean Model, version 4 (Griffies et al., 2005), with gradually increased resolutions (1° to 1/3°) from 30°S or 30°N to the equator and 40 levels in the vertical direction. The sea-ice component is the GFDL Sea Ice Simulator (Winton, 2000), which has the same resolution as that in the oceanic component.The S2S hindcast experiment designed by the S2S project is to carry out the prediction by initializing on each day from 1 January 1994 to 31 December 2013 and integrating the model for 60 days. To reduce the uncertainty of initial conditions, for each initial date the strategy of four members using initial conditions at 0000, 0600, 1200 and 1800 UTC is adopted. Since the ensemble mean generally has higher skill, we focus on the outputs of the ensemble mean from these four members. Daily averaged variables are used in this study.
Because the hindcast experiment of S2S prediction was conducted each day for the future 60 days during the period of 1994-2013, the model outputs contain a continuous distribution of hindcasts on all dates from 1994 to 2013, at all lead times from 1 to 60 days. Thus, a time series of model data at any fixed lead time can be formed to compare against the observation for understanding the model forecast skill at that lead time.
2
2.2. Observational data and reanalysis products
The observational and reanalysis datasets used for model evaluation include: (2) daily outgoing longwave radiation (OLR) data (Liebmann and Smith, 1996), which are used as a proxy for convection, provided by the National Oceanic and Atmospheric Administration (NOAA); (3) the ERA-Interim dataset (Dee et al., 2011) of the ECMWF; (4) the NCEP-NCAR reanalysis (Kalnay et al., 1996); and (4) the NOAA High-resolution Blended Analysis of Daily Sea Surface Temperature (SST; Reynolds et al., 2007). To reduce the uncertainty of reanalysis products, we averaged the ERA-Interim and NCEP-NCAR fields for horizontal wind components, 3D specific humidity and surface sensible heat flux (SHF). The resolution in the reanalysis datasets and model outputs are unified on a 2.5°× 2.5° grid over eight layers (1000, 925, 850, 700, 500, 300, 200, and 100 hPa). The summer season is defined as May to October.2
2.3. Evaluation metrics
To better evaluate the real-time forecast capability of the BSISO, we adopt the multi-variate (MV) EOF analysis proposed by (Lee et al., 2013) to derive the real-time BSISO index from different lead forecasts. Our evaluation focuses on the spatial pattern and temporal evolution of the real-time BSISO index in the BCC S2S model output. As in (Lee et al., 2013), we first remove the climatology (the first three harmonics) and the interannual variability (using the running mean of the last 120 days) from daily OLR and 850-hPa zonal wind (U850) fields in the Asian monsoon region (10°S-40°N, 40°-160°E). We then normalize and decompose the two anomaly fields, using the MV-EOF analysis. The first two modes (EOF1 and EOF2) can effectively represent the north-northeastward propagation of the 30-90-day BSISO. Using the time series (PC1 and PC2) of these two modes, we can then divide the BSISO lifecycle into eight phases (Wheeler and Hendon, 2004; Lee et al., 2013), as in Table 1. When (PC12 + PC22)1/2>1.0, it is considered as an active BSISO event. The third and fourth EOF modes represent the northwestward-propagating signal with a shorter period of 10-30 days (Lee et al., 2013). In our model assessments, the filtering and EOF analysis were carried out separately for each lead time. (Jie et al., 2017) noted that the operational models generally show lower skills in predicting this 10-30-day mode. Our analysis also finds that the skill of the 10-30-day BSISO prediction of BCC S2S model is about nine days (not shown), similar to the range of weather forecasts. Therefore, only the 30-90-day BSISO mode is investigated in this study.To quantitatively evaluate the simulation and forecast skills of the 30-90 day BSISO, we calculate the PCC, root-mean-square error (RMSE) and bivariate anomaly correlation coefficient (ACC) of PC1 and PC2. These metrics have been used widely in evaluating MJO/ISO prediction (Lin et al., 2008a; Xiang et al., 2015; Liu et al., 2017). The definitions of these metrics are given below: \begin{eqnarray*} {\rm PCC}&=&\frac{\sum_{j=1}^M\sum_{i=1}^N[(a_{i,j}-\bar{a})(b_{i,j}-\bar{b})]} {\sqrt{\sum_{j=1}^M\sum_{i=1}^N(a_{i,j}-\bar{a})^2\sum_{j=1}^M\sum_{i=1}^N(b_{i,j}-\bar{b})^2}} ;\\ {\rm RMSE}&=&\sqrt{\frac{1}{MN}\sum_{j=1}^M\sum_{i=1}^N(a_{i,j}-b_{i,j})^2} ;\\ {\rm ACC}(\tau)&=&\frac{\sum_{t=1}^T(a_{1t}b_{1t}+a_{2t}b_{2t})} {\sqrt{\sum_{t=1}^T(a_{1t}^2+a_{2t}^2)}\sqrt{\sum_{t=1}^T(b_{1t}^2+b_{2t}^2)}} . \end{eqnarray*} Here, a is observation and b is model output; i and j refer to longitude and latitude points, respectively; M and N indicate the numbers of longitude and latitude points, respectively. An overbar in the definition of PCC represents an area average. The subscript t is time, and τ is the forecast lead time. All the time points in May-October during the period of 1994-2013 were included for the ACC calculation. The subscripts 1 and 2 are for the first and second variables (such as PC1 and PC2), respectively.
2
2.4. Diagnosis of BSISO propagation
Several mechanisms associated with atmospheric internal dynamics and air-sea (land) coupling processes have been proposed to explain the northward propagation of the BSISO [see reviews by (DeMott et al., 2012) and (Wang, 2012)], which are summarized in Fig. 1. In observations, the positive vorticity and moistening in the lower troposphere appear to the north of BSISO convection (Hsu and Weng, 2001; Jiang et al., 2004; Tsou et al., 2005; Bellon and Sobel, 2008), favoring northward development of the convection. Two mechanisms may account for the generation of positive vorticity to the northern flank of the convection over the Indian Ocean (Jiang et al., 2004; Bellon and Sobel, 2008). According to (Jiang et al., 2004), the interaction between background vertical wind shear and the meridional gradient of anomalous baroclinic divergence can induce a positive anomaly of barotropic vorticity (ζ'+) to the north of the convection. This process can be written as \begin{equation} \frac{\partial{\zeta}'_+}{\partial t}\propto\bar{u}_T\frac{\partial{D}'_-}{\partial y} , \ \ (1)\end{equation} where an overbar represents the seasonal-mean component and a prime indicates the anomalous component related to the BSISO, which is extracted from the 30-90-day Lanczos filtering (Duchon, 1979). In Eq. (2), ζ'+ indicates the barotropic mode of the vorticity anomaly; $\bar{u}_T$ is the seasonal-mean vertical wind shear, defined by the difference in zonal wind between 200 and 850 hPa $(\bar{u}_200-\bar{u}_850)$; and D'- represents the baroclinic mode of divergence (anomalous divergence between 200 and 850 hPa). Accompanied by the enhanced convective heating anomaly, anomalous ascending motion appears in the mid troposphere and a divergence (convergence) anomaly can be observed in the higher (lower) level (D'->0). The maximum D'- is located at the convective center; thus, the meridional gradient of baroclinic divergence is negative (positive) to the north (south) of the convection. Under the easterly shear over the Asian summer monsoon region, a positive (negative) barotropic vorticity may be induced in the northern (southern) part of the convective center.
(Bellon and Sobel, 2008) emphasized the role of vorticity advection by the background meridional flows in the generation of a positive vorticity anomaly to the north of the convection. The process can be written as \begin{equation} \frac{\partial{\zeta}'_+}{\partial t}\propto-\bar{v}_T \frac{\partial{\zeta}'_-}{\partial y} , \ \ (2)\end{equation} where $\bar v_T$ is the seasonal-mean vertical wind shear of meridional wind between 200 and 850 hPa $(\bar v_200-\bar v_850)$ and ζ'- represents the baroclinic mode of vorticity anomaly (difference in the vorticity anomaly between 200 and 850 hPa). Over the BSISO convective area, the positive vorticity anomaly at the low level is coupled with the negative vorticity anomaly at the high level, forming a baroclinic vorticity couplet. As the background monsoonal northerly (southerly) in the higher (lower) troposphere works on the baroclinic couplet of the vorticity anomaly, a positive (negative) vorticity anomaly is induced over the northern (southern) region of the convection.
The leading phase of planetary boundary layer (PBL) moisture is also favorable for the northward propagation of the BSISO, as it destabilizes the atmospheric column (Hsu and Weng, 2001; Kemball-Cook and Wang, 2001; Jiang et al., 2004; Hsu and Li, 2012). The PBL moistening can result from the PBL moisture convergence related to the low-level vorticity anomaly and the meridional advection of moisture. Both the seasonal-mean flows and perturbation wind and moisture fields can affect moisture advection. To diagnose the major processes, we decompose the specific humidity and meridional wind fields into the summer-mean component (with an overbar) and the BSISO-related perturbation part (with a prime). Thus, three processes associated with mean flow and perturbation interactions can cause meridional advection of moisture, as expressed by Eq. (4) below: \begin{equation} -\left(v\frac{\partial q}{\partial y}\right)'=-\left(\bar{v}\frac{\partial{q}'}{\partial y}\right)' -\left({v}'\frac{\partial\bar{q}}{\partial y}\right)'-\left({v}'\frac{\partial{q}'}{\partial y}\right)' , \ \ (3)\end{equation} where q is the specific humidity. In addition to the internal atmospheric dynamics, the surface warming induced by anomalous surface heat fluxes to the north of the convection also plays a role in the northward propagation via destabilization of the lower atmosphere (Webster, 1983; Kemball-Cook and Wang, 2001; Fu et al., 2003). (Webster, 1983) pointed out that increased SHF over the South Asian peninsula could reduce the stability of the low-level atmosphere and favor the northward propagation of convection near the equator. (Kemball-Cook and Wang, 2001) and (Fu et al., 2003) found that the leading phase of the SST anomaly related to anomalous evaporation or latent heat flux (LHF) also favors induction of BSISO northward propagation.
In section 4, all these processes are diagnosed and compared using the observation and model outputs at different lead times. This will help us identify the key factors governing the northward propagation of the BSISO in the BCC S2S model.

Figure 3 compares the BSISO modes predicted by the model against the observation using the MV-EOF (Lee et al., 2013). For the observations, the EOFs and PCs are derived from a filtered time series that represents the evolution of the "real world." For the data at a lead time of n days, the time series represents the state after the model has been integrated for n days from an observed initial state, and is thus influenced by the varying initial states and the model evolution from those initial states. In Fig. 3a, the observed first mode (EOF1) shows the convective center (dark blue areas) located over the equatorial Indian Ocean, while there is anomalous descending motion (dark red areas) over the SCS and WNP. A cyclonic anomaly and an anticyclonic anomaly appear over the northern flanks of the convective center and the descending motion, respectively (Fig. 3a), which favor the northward propagation of the BSISO convection (Jiang et al., 2004; Tsou et al., 2005). The observed second mode (EOF2) of the BSISO represents the northeastward propagation of the convection and circulation shown in EOF1. The convection over the equatorial Indian Ocean has moved to the Indian subcontinent and the Bay of Bengal, while the descending motion over the SCS and WNP moves to the east of Taiwan Island and its strength is weakened. With the new anomalous descending branch appearing over the equatorial Indian Ocean, a clear wave train oriented southwest-northeast is formed (Fig. 3h). In the model, the Lead 1d signal of the BSISO propagates toward the northeast (Figs. 3b and i); however, as the lead time increases, the wave train of BSISO convection is oriented along the east-west direction. Both predicted EOF1 and EOF2 modes move eastward, without the northward development seen in the observation (Figs. 3c-g and 3j-n). Note that the patterns of EOF modes do not change much after a longer integration (such as Lead 11d), which resemble the inherent BSISO modes in the climate simulations of BCC_CSM1.1 (not shown). The results suggest that the loss of forecast skill with longer lead times is likely due to the model adjusting toward its intrinsic modes. It is also found that the prediction skills of BSISO-related convection in terms of PCC and normalized RMSE are lower than those of circulation (U850), suggesting the BCC S2S model has more difficulty in capturing convection anomalies than circulation anomalies over the Asian monsoon region.

The ACC of the forecast PC1 and PC2 changes with time (Fig. 4a). Using ACC = 0.5 as the threshold——namely, the model is considered skillful when the ACC is larger than 0.5——the results show that the forecast skill of the model drops as lead time increases. The results show that the BCC S2S model remains skillful when the lead time is about 10 days; the ACC skill is decreased to below 0.5 when the lead time is greater than 11 days. Figure 4a compares the ACC of PCs derived from each lead-time simulation against the observed PCs. Some studies obtained model PCs by projecting simulated anomalous fields onto the observed EOF modes (Fu et al., 2013; Kim et al., 2014). To confirm the model skill, we also use the projection method to obtain the simulated PCs and then computed their ACC. It appears that the ACC also drops below 0.5 after 11 days (Fig. 4b), consistent with the results in Fig. 4a. We compare the individual ACCs for PC1 and PC2, respectively, and find that the ACC of PC1 (12 days) is slightly higher than that of PC2 (9 days) (not shown).

The lead-lag correlation coefficient of PC1 and PC2 (black curve in Fig. 5) shows that PC1 leads PC2 by about 12 days in the observation, while their simultaneous correlation coefficient is zero. This indicates the two modes with a 90° phase difference represent a propagation mode (Lee et al., 2013). The model results show that PC1 leads PC2 by 12-13 days. At longer than Lead 11d (red dashed curve in Fig. 5), the correlation is smaller than 0.2. In addition, at Lag -5d, the correlation between the two EOF modes changes from zero to negative, showing that the model fails to simulate the propagation observed. As the lead time increases, the lead-lag correlation coefficient decreases, and the model fails to forecast the 90° phase difference observed (dashed curve in Fig. 5). This indicates the model forecasts fail to maintain the observed PC1 and PC2 evolutions from the initial state.

To confirm that the BCC S2S model can only predict the BSISO within 10 days in advance, we carry out spectral analysis of OLR at Lead 11d and compare it with the observation and the results of Lead 4d, which show skillfully predicted BSISO. In Fig. 6, the observed power spectra of PC1 and PC2 show clear 30-90-day signals of the BSISO (Figs. 6a and d). The 30-90-day signals are significant in the outputs of Lead 4d, although their variances are weaker than the observed (Figs. 6b and e). At Lead 11d, the model is unable to predict the convection amplitude of the BSISO (Figs. 6c and f).

As for the propagation features, we compare the life cycle of BSISO convection for observation, Lead 4d and Lead 11d, based on the phase composites of active BSISO days [( PC12+ PC22)1/2>1]. Note that the phase composites for model outputs were constructed using the PC1 and PC2 predicted by the model at different forecast leads. Unsurprisingly, the model does a better job for a short lead time (such as Lead 4d) with a distinct northward propagation over the Asian monsoon region (middle rows in Fig. 7), similar to the observation (left-hand rows in Fig. 7). As seen in the observation, the convection center of the BSISO is located over the equatorial Indian Ocean during Phase 1. The anomalous descending branch is over the Bay of Bengal, the SCS and the WNP, while there is a relatively weak convection area south of Japan, showing a wave train along the northeast-southwest direction. This wave train moves northward and northeastward over time. The anomalous convection originally over the equatorial Indian Ocean moves to the Indian subcontinent during Phases 2 and 3, reaches the SCS during Phases 4 and 5, and gradually weakens thereafter. An anomalous descending branch can be seen forming over the equatorial Indian Ocean during Phases 4-6 and then moving north/northeastward. As the forecast lead time increases, the BSISO structures become weakened and less organized. At Lead 11d (right-hand rows in Fig. 7), the BSISO signal shows a weak strength compared to the observed and moves mainly eastward, suggesting that the northward propagation component is rather weak. Again, these features resemble the intrinsic model BSISO. The low ACC skill, which represents the capability in describing the first two PCs of MV-EOF analysis and propagating mode, may be attributable to the model failure of forecasting the northward propagating component at Lead 11d and longer lead times.


According to (Jiang et al., 2004) and (Bellon and Sobel, 2008), the positive vorticity is generated by the meridional distribution of baroclinic divergence $(\partial D'\_/\partial y)$ interacting with the background easterly shear $(\bar{u}_T)$, as well as the summer-mean northerly shear $(\bar{v}_T)$ working on the meridional gradient of the baroclinic vorticity anomaly $(\partial \zeta'\_/\partial y)$. Figure 9 shows the diagnostic results of these two processes in observation and model predictions. The model is skillful at predicting the background easterly and northerly shears over the Asian monsoon region even at a longer lead (Lead 11d), while it tends to underestimate the meridional gradients of baroclinic divergence and vorticity modes associated with the BSISO (Fig. 9). The small values of $\partial D'\_/\partial y$ and $\partial \zeta'\_/\partial y$ may result from weak BSISO variability in the BCC S2S model (Figs. 1, 6 and 7).

Through reducing atmospheric stability, the leading phase of PBL moistening is another factor promoting the northward propagation (Hsu and Weng, 2001; Jiang et al., 2004; DeMott et al., 2012). The increased moisture in the lower troposphere is observed to the north of the convection as it propagates northward in the observation (Fig. 8d), and in the results of Lead 4d (Fig. 8e). In the Lead 11d prediction, the moisture anomaly is relatively weak (Fig. 8f). This weakened moistening is partly related to weak moisture convergence (Fig. 10c) induced by the underestimated low-level vorticity anomaly (Fig. 8c). In addition to moisture convergence (Figs. 10a-c), the increased moisture comes from meridional moisture advection (Jiang et al., 2004), although this effect (Figs. 10d-f) is smaller than moisture convergence (Figs. 10a-c). In the Lead 11d hindcast, there is negative moisture advection to the north of the convection (Fig. 10f), hindering the northward development of the convection.

We further analyze the causes of biases of moisture advection using the terms on the right-hand side of Eq. (4). The results show that the negative meridional moisture advection results mainly from the advection of anomalous moisture by the mean flow [-($\bar v(\partial q'/\partial y)$)'], which leads to a large negative moisture anomaly to the north of the convection at Lead 11d (Fig. 11). These terms at Lead 4d show a small negative anomaly, which are offset by a larger positive anomaly induced by -[v'($\partial\bar q/\partial y$)]'. Based on the geographical distributions of summer-mean flow and the 30-90-day moisture anomaly, we find that the model reasonably predicts the monsoonal southerly in the lower troposphere but fails to forecast the meridional gradient of the BSISO-related moisture anomaly (not shown).

The warm SST anomaly (induced by reduced LHF) and increased SHF are conducive to convection development through destabilizing the lower troposphere (Webster, 1983; Kemball-Cook and Wang, 2001; Fu et al., 2003). The BCC S2S model is able to produce the positive SST and SHF anomalies over the northern region of the BSISO convection. The amplitudes of SST and SHF anomalies, however, decrease as lead time increases (Fig. 12). Therefore, the effect of surface heat fluxes is too weak to drive the northward propagation at the lead time of 11 days.

The model underestimates the variability of the BSISO signal, even at the lead time of 1 day (Fig. 2), which deteriorates quickly as the lead time increases. The evaluation of the first two MV-EOF modes of the northeastward-propagating 30-90-day BSISO (Lee et al., 2013) shows that the northward propagation component weakens as the lead time increases, resulting in biasedly strong eastward propagation (Fig. 3). The patterns of the first two EOF modes at longer lead times resemble those of the BSISO modes in the long-term (20 years) simulation by the model, suggesting that the loss of prediction skill is related to the BCC model adjusting toward its intrinsic BSISO mode. Analysis of the time series (PC1 and PC2) of the first two modes reveals that the model can forecast the lead-lag correlation correctly ~10 days in advance (Figs. 4 and 5). After Lead 11d, the 30-90-day signal becomes insignificant (Fig. 6) and the northward propagation is no longer obvious (showing an eastward propagation feature instead) (Fig. 7). (Jie et al., 2017) evaluated the forecast skill of the BSISO in 10 operational models that participated in the S2S project and found that these models can predict BSISO events up to 6-24.5 days in advance. Compared to other operational models, the forecast skill of the BCC S2S model is relatively weak.
To understand the factors that inhibit the BSISO prediction skill, the main mechanisms responsible for the northward propagation proposed by some observational studies were diagnosed using the outputs of the BCC S2S model. The insignificant northward propagation at Lead 11d arises from the weak anomalies of cyclonic vorticity, PBL moistening and surface heat fluxes to the north of the convection (Figs. 8 and 12), which reduce the northward propagation of the BSISO over the tropical Indian Ocean. The model can predict the summer-mean easterly and northerly shears associated with the generation of the leading phase of barotropic vorticity (Jiang et al., 2004; Bellon and Sobel, 2008), while it fails to forecast the BSISO-related anomalies of divergence and vorticity disturbances at a longer lead time (Fig. 9). Although the monsoon circulation is well predicted, the model's capability in reproducing observed meridional gradient of moisture anomaly is also limited. The biased distribution of the moisture anomaly causes negative meridional moisture to the north of the convection at Lead 11d (Figs. 10 and 11). Moreover, the amplitudes of warm SST and SHF decrease as the forecast lead time increases (Fig. 12). All these weak anomalies of moistening and warming in the lower troposphere seem to work against the generation of enhanced instability in the atmospheric column over the northern flank of the BSISO convection, leading to the lack of northward propagation at Lead 11d.
In summary, the limited skill of BSISO prediction by the BCC S2S model may be attributable to the defects in forecasting the amplitude and propagation features of BSISO signals. Compared to the BSISO circulation, the model shows lower skill in predicting the convection associated with the BSISO. Studies have shown that adjustment and improvement of parameterization schemes for convection in the model are the keys for improving the simulation and prediction of BSISO moisture and convection (Lin et al., 2008b; Kim et al., 2014; Jiang et al., 2015). Based on the diagnostic methods used in this study, we plan to test model sensitivity to parameterization schemes and to the key parameters in these schemes in terms of BSISO simulation and forecast skills using the BCC S2S model. Applying the same diagnostic methods to other S2S models (Vitart et al., 2017) with various forecast skills (Jie et al., 2017) can also advance our understanding of the factors controlling the BSISO prediction in state-of-the-art operational models. (Liu et al., 2017) found that the MJO prediction skill of the BCC S2S model increased from 16 days to 22 days when improved initial states were used. It is also part of our ongoing work to examine the effects of initial conditions on the BSISO prediction skill.