Predictability and Prediction Skill of the Boreal Summer Intra-Seasonal Oscillation in BCC_CSM Model

The boreal summer intra-seasonal oscillation (BSISO) is the predominant sub-seasonal variability over the East Asia (EA) and western North Pacific (WNP) region and critical for seasonal forecast of the EA summer mon -soon. This study examines the theoretically estimated predictability and practical prediction skill of the EAWNP BSISO in the Beijing Climate Center Climate System Model version 2 (BCC_CSM2.0), which is one of the par-ticipants in the Sub-seasonal to Seasonal Prediction Project. The results from the uninitialized free run of BCC_ CSM2.0 show that the model reasonably simulates EAWNP BSISO in terms of its variance, propagation, and structure. Measured by the bivariate correlation (> 0.5) and root mean square error (< 2 ) between the predicted and observed real-time BSISO index, the prediction skill and predictability of EAWNP BSISO are about 14 and 24 – 28 days, respectively. The initial/target strong BSISO cases have a relatively higher prediction skill than the initial/target weak BSISO cases. For the theoretically estimated BSISO predictability, a similar dependence on target amplitude occurs in the model, while no significant dependency on initial amplitude is found. Moreover, diagnosis of the phase dependence reveals that BSISO is less skillful for the prediction starting from active or active-to-break transition phases of WNP rainfall, whereas it is more predictable when prediction is targeting extreme dry/wet phases of WNP rainfall. Finally, systematic errors are found in BCC_CSM2.0 such as the under estimation of BSISO amplitude and the faster phase speed.


Introduction
The boreal summer intra-seasonal oscillation (BSISO) is a leading mode of the East Asia and western North Pacific (EAWNP) summer monsoon system. Unlike its boreal wintertime counterpart, the Madden-Julian Oscillation (MJO; Madden and Julian 1972) that is featured by equatorial eastward propagation, BSISO over EAWNP exhibits significant northward/ northeastward propagation with a periodicity of 10 -90 days (Lau et al. 1988;Hsu and Weng 2001). It significantly modulates the wet and dry spells of the EAWNP summer monsoon and has a great impact on the evolution of Meiyu-baiu frontal systems, the advance/retreat of the WNP subtropical high, and the genesis/intensification of tropical cyclones (Chen et al. 2009;Mao et al. 2010;Hsu et al. 2016). Hence, the skillful prediction of BSISO can fill the gap between the synoptic and seasonal prediction and gain an insight into the sources and errors of extended-range forecast over the East Asian region. Motivated by this understanding, BSISO simulation and forecast using dynamical models has received increasing attention, with the aim of reducing social and economic losses (Brunet et al. 2010).
However, the current modeling and forecasting capability for the space-time characteristics of BSISO remain limited (Lin J.-L. et al. 2008;Sperber and Anna malai 2008;Joseph et al. 2010;Sabeerali et al. 2013;Fang et al. 2017a, b). Past multi-model assessments of the model performance of BSISO indicated that general circulation models have difficulty in capturing the salient features of BSISO by typically exhibiting weak northward propagation or a standing oscillation. The BSISO forecast is even more challenging because prediction on this time scale shows sensitivity to both the atmospheric/ocean initial conditions and the lower boundary conditions. In recent years, sub-seasonal forecast using dynamical forecast models has become an intensive research topic and has been documented in several studies (Waliser et al. 2003;Liess et al. 2005;Vitart and Molteni 2009;Fu et al. 2013;Kang et al. 2014;Lee et al. 2015). Most of these studies found that useful skills of BSISO up to 10 -20 days of forecast lead time have been achieved for most models. However, the forecast skill of BSISO is still far lower than its theoretically esti-mated predictability, which is 25 -40 days or beyond (Ding et al. 2011). Therefore, to further improve the BSISO prediction skills, it is crucial to systematically evaluate present-day dynamical BSISO prediction capabilities and to gain an insight into the weaknesses of dynamical prediction systems.
Similar to the winter MJO (Madden and Julian 1972), the prediction skill of BSISO depends on its initial amplitude, with higher skill occurring for forecasts initialized from larger amplitude (Seo et al. 2005;Fu et al. 2009;Abhilash et al. 2014). The predictability and prediction skill of BSISO were also found to show a strong dependence on its phase. Goswami and Xavier (2003) examined the predictability of rainfall associated with BSISO over the Indian region based on observational data. They found that the skills for the active and break monsoon phases are about 10 and 20 days, respectively, which suggests that the breakto-active transition phase is less predictable than the other phases, known as the "monsoon prediction barrier". Recently, Lee et al. (2015) assessed the potential predictability and prediction skill of BSISO over the Indian monsoon region using six coupled models in the Intra-seasonal Variability Hindcast Experiment (ISVHE) project. They found that while the predictability and prediction skill of BSISO are dependent on the initial phase and month, the degree of dependence varies by model. Some studies found that the multi-model ensemble approach and inclusion of air-sea coupling can extend the BSISO predictability by 2 -7 days (Fu et al. 2007(Fu et al. , 2013Lee and Wang 2016).
While most of the previous work mentioned above focuses on BSISO over the Indian region, BSISO predictability and prediction skill over the EAWNP region are less well known. Recently, Lin (2013) and Lee and Wang (2016) developed a BSISO index based on the first two combined empirical orthogonal function (CEOF) modes of outgoing longwave radiation (OLR) and 850-hPa zonal wind (U850) over the EAWNP region. Similar to the Wheeler-Hendon MJO index (Wheeler and Hendon 2004), by using this index, we can isolate the BSISO signal from the forecast and observational data without applying a temporal filter. This method, therefore, is viable for real-time EAWNP BSISO monitoring and prediction. Lee and Wang (2016) also applied this approach to hindcast datasets from the ISVHE project to evaluate the forecast skill of EAWNP BSISO. The results indicate that EAWNP BSISO is less predictable than BSISO over the Indian or the entire Asian monsoon region. The multi-model mean of the BSISO predictability and prediction skill over EAWNP is about 33 -37 and 8 -10 days, whereas those over the Asian or the Indian monsoon region can reach up to 40 -45 and 9 -17 days, respectively.
In this study, this approach will be applied to a set of hindcast experiments conducted by the Beijing Climate Center Climate System Model version 2.0 (BCC_CSM2.0). BCC_CSM2.0 is one of the participants in phase 5 of the Coupled Model Intercomparison Project (CMIP5) and the Sub-seasonal to Seasonal (S2S) Prediction Project. Some studies have reported the basic performance in MJO simulation and prediction for the BCC_CSM2.0 model (Zhao et al. 2015;Liu et al. 2017). They found that the models can capture salient characteristics of MJO in terms of its intensity, periodicity, propagation, and structure, although the MJO intensity is underestimated and the power of westward propagation is overestimated. Meanwhile, analysis based on hindcast data from S2S indicated that BCC_CSM2.0 shows useful MJO prediction skill up to 21 -22 days, which is comparable to the performance of several operational S2S models. However, up to now, its BSISO predictability and prediction skill over the EAWNP region have not been examined, which served as the primary motivation for the present study.
The rest of this article is organized as follows. Section 2 introduces the model, reanalysis data, and the BSISO characteristics in free simulations. Section 3 shows details of the hindcast and verification method. The BSISO predictability and forecast skill over the EAWNP region are discussed in Section 4, and the results are summarized in Section 5.

Model
The model used in the present study is the BCC_ CSM2.0 one developed at the Beijing Climate Center (BCC), China Meteorological Administration, which consists of fully coupled components of the atmosphere, ocean, ice, and land. The atmospheric component in BCC_CSM2.0 is the BCC Atmospheric General Model version 2.1 (BCC_AGCM2.1) with a T106 horizontal resolution and 26 hybrid sigma/pressure layers in the vertical direction . Based on the Community Atmosphere Model version 3 at the National Center for Atmospheric Research, BCC_AGCM2.1 is developed by incorporating many modifications to the dynamical framework and model physics parameterizations such as deep cumulus convection, adiabatic adjustment, snow cover, and sea surface heat flux (Wu et al. 2008;Wu 2012). The oceanic component is the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model version 4 (MOM4_L40) with 40 levels in the vertical direction and a horizontal resolution of 1° × 1° poleward of 30°S -30°N gradually descending to 0.33° between 30°N and 30°S. The sea ice component is the GFDL Sea Ice Simulator (SIS), and the land model is the BCC Atmosphere and Vegetation Interaction Model version 1.0 (BCC_AVIM1.0) (Ji 1995). The model components are coupled without any flux correction. BCC_CSM2.0 has been used for short-term climate prediction in BCC, and the performance is reliable (Liu et al. 2014(Liu et al. , 2015.

Dataset
In this study, the zonal and meridional winds at 850 hPa from the NCEP-NCAR reanalysis and the OLR from the National Oceanic and Atmospheric Administration (NOAA) (Liebmann and Smith 1996) for the period of 2000 -2013 are utilized to calculate the aspects of BSISO and climatological features over the EAWNP region for the model validation. The observational precipitation and sea surface temperature (SST) data (2000 -2013) used in this study are the Tropical Rainfall Measuring Mission (TRMM; Huffman et al. 2007) 3B42 rainfall dataset version 7 with 0.25° × 0.25° spatial resolutions and daily SST from the NOAA optimum interpolation SST (OISST) version 2 (Reynolds et al. 2007), respectively.

Free simulation of EAWNP BSISO
Realistic BSISO simulation provides a solid foundation for effective BSISO prediction. Therefore, to reveal the fidelity of BCC_CSM2.0 in simulating EAWNP BSISO, a free integration was performed over a 20-year period with the same greenhouse-gas external forcing as that in the CMIP5 historical simulation before conducting the hindcast experiment. The data from the last 16 years of the free run described are used to represent the climatological and intra-seasonal variability behavior in the model. To extract the BSISO component, a 20 -70-day band-pass filter is applied to both the model simulation and observation data.
A realistic background state in the model is of crucial importance to the model's representation of BSISO features. Figure 1, therefore, firstly presents the simulated and observed mean OLR and 850-hPa wind field during May-September over the EAWNP region. BCC_CSM2.0 can reasonably simulate the cyclonic low-level circulation, which is related to the South China Sea monsoon trough. The expansive anti-cyclonic circulation over the WNP region associated with the WNP subtropical high is also well reproduced (Fig. 1b). The model captures the band of low OLR extending from WNP westward across the South China Sea and the Indo-China Peninsula, which corresponds to the WNP summer monsoon, but OLR is slightly overestimated over the eastern part of the Bay of Bengal (Fig. 1b). Figure 2 shows the simulated and observed May-September mean precipitation and SST over the EAWNP region. The SST climatology in the model is comparable to the observations, especially over the warm pool region south of 20°N, where both the observed and modeled SST exceed 28°C (Fig. 2a). The spatial pattern correlation coefficient between the simulation and the observation is 0.92. However, the model tends to underestimate SST over the WNP region with a bias as large as approximately 1°C (Fig.  2b). This cold bias is common in CMIP5 models (Song and Zhou 2014). Moreover, the model captures the gross pattern of mean precipitation, but dry biases are found over the East Asian region, and wet biases appear over most of the WNP region.
The geographic distribution of simulated and observed BSISO intensity over the EAWNP region is shown in Fig. 3. The BSISO intensity is represented by the variance of the 20 -70-day filtered U850 anomalies during May-September. In observations, the BSISO intensity exhibits maxima located over the Philippine Sea, the South China Sea, and the eastern Bay of Bengal. Another secondary maximum is also evident over the sea south of Japan (Fig. 3a). It is found that the model patterns of the variance look largely similar to observations despite some deficiencies in the magnitudes. The model overestimates the BSISO intensity over the Bay of Bengal and the South China Sea, while it underestimates it over the east Philippine Sea and the East Asian region (Fig. 3b). These biases are similar to those in BCC_AGCM2.1 (Fang et al. 2017b), which is the atmospheric component of BCC_CSM2.0.
The spatial structure of BSISO over EAWNP represented by the first two CEOF modes of the daily anomalies of 20 -70-day filtered OLR and U850 averaged between 90°E and 150°E is displayed in Fig. 4. For the observations, the first CEOF mode features enhanced convection centered near 15°N. To the north (south) of this convection center, easterly (westerly) anomalies prevail. As for the second CEOF mode, both convection and wind patterns are in close quadrature with CEOF1 (Figs. 4a,b). This leading pair of CEOFs jointly represents 60.5 % of the total variance, similar to the analysis of Lin (2013). The simulated first two leading modes explaining 42.2 % of the filtered variance show convection and low-level circulation patterns matching well with the observational counterpart, although the magnitude of both the OLR and U850 anomalies over the East Asian region is overestimated (Figs. 4c,d). The corresponding prin-cipal components (PCs) are also in quadrature with PC1 leading PC2 by 9 and 11 days for the simulation and observation, respectively (figure not shown).
Given that northward propagation is a fundamental characteristic of EAWNP BSISO, Figure 5 depicts the lag correlation of 20 -70-day filtered U850 and precipitation along 122 -135°E with time series of the 20 -70-day filtered precipitation anomalies averaged over the Philippine Sea and the South China Sea (10 -20°N, 110 -130°E). The model well captures the observed northward propagation of BSISO-related U850 and precipitation anomalies from the equatorial western Pacific to East Asia. The phase relationship is also reasonably reproduced, with easterly (westerly) anomalies leading (lagging) the BSISO convection center by 3 -5 days (Fig. 5). However, compared with the observations, the BSISO signals appear to be weaker over the equatorial western Pacific and stronger over East Asia. It should also be noted that the phase speeds of BSISO in the model are slightly faster than the observations. Figure 6 further presents the composite BSISO life cycle using precipitation and 850-hPa wind anomalies. Following Lin (2013) and Lee and Wang (2016), a BSISO event can be divided into eight phases according to the phase map of PC1 and PC2 associated with the CEOF1 and CEOF2 modes of EAWNP BSISO. For each phase, only the days satisfying the following two criteria are selected for the composite: (a) the BSISO phase angles are within the phase; (b) the BSISO amplitudes are greater than 1.5 (64 % of all summer days for the observation versus 60.7 % for the model). The observed BSISO cycle is featured by a northward (or northwestward) propagation of a southwest-northeast tilted convection band from the equatorial Pacific to East Asia. During the northward propagation process, enhanced (suppressed) convection takes place in the southern (northern) portions of zonally elongated anomalous cyclones (Fig. 6a). As mentioned in previous research (e.g., Hsu and Weng 2001;Mao et al. 2010), such a circulation-convection pattern is consistent with the theoretically analytical solution for an off-equatorial heating source (Gill 1980). BCC_CSM2.0 shows obvious northward propagation of the intra-seasonal rain band in spite of regional differences in magnitude. The alternation of low-level cyclonic and anti-cyclonic anomalies occur-   ring over the WNP region is particularly well captured in the model. However, the rain band of the model is zonally oriented rather than tilted, and the BSISO signal over the Bay of Bengal is not properly simulated as seen in phases 1, 7, and 8 (Fig. 6b). This may be associated with a weaker eastward-propagating MJO during the boreal summer, as revealed from the examination of MJO in four versions of BCC_CSM by Zhao et al. (2015). They also found that the eastward-propagating MJO signal shows a faster phase speed and is especially weaker when MJO propagates into the Maritime Continent and the East Indian Ocean.

Hindcast experiments
The realistic BSISO simulation in BCC_CSM2.0 seen above motivates us to further conduct prediction experiments. In the prediction system, the atmospheric component is initialized by adopting a fast nudging technique that forces the model solution toward the analysis. The atmospheric fields, including winds, temperature, geopotential height, and surface pressure, obtained from the National Center for Environmental Prediction's Final Operational Global Analysis data (NCEP FNL) are nudged for all vertical levels with a 450 s relaxation time scale. To reduce the coupling shock at the beginning of the forecast experiment, an initialization was also introduced for the underlying surface atmospheric state through a similar nudging scheme. The initialized variables include the sea level pressure, surface air temperature and wind, downward surface longwave and shortwave radiation from the NCEP FNL analysis, and precipitation rate from the BCC merged precipitation dataset with a 0.5° × 0.5° horizontal resolution at a 3-h time interval . The ocean initialization uses data assimilation products made available by BCC Global Ocean Data Assimilation System version 2 (BCC_GODAS2.0; Zhou et al. 2016), which is based on MOM4_L40 and a three-dimensional variational data assimilation system. Multiple sources of observational sea level anomalies and SST are real-time assimilated in BCC_ GODAS2.0. Details about the assimilation system and the associated data descriptions please refer Zhou et al. (2016).
Since the purpose of this study is to explore the BSISO prediction skill over the EAWNP region, hindcast runs were carried out on the 1st, 6th, 11th, 16th, 21st, and 26th of each month during 2000 -2013 from May to September. To reduce the uncertainties when initializing the prediction system, for each initial date, a daily four-member time-lagged ensemble was adopted with atmosphere initial conditions at 0000, 0600, 1200, and 1800 UTC. Therefore, in total, we made 420 (14 years × 5 months × 6 cases) forecast cases during the 14 years; each date has four members, so that we totally have 1680 sets of forecast integrations. For each run, the model was integrated for 60 days. More details of the hindcast experiment can be found in Liu et al. (2017).

Hindcast methodology
Following Lin (2013) and Lee and Wang (2016), the prediction skill of BSISO over the EAWNP region was examined by using the real-time BSISO indices. The EAWNP BSISO indices are the first two PC time series (PC1 and PC2) of the CEOF of the observational daily OLR and U850 anomalies averaged over 90 -150°E during May-September from 2000 to 2013. Before the CEOF calculation, the slow annual cycle (the annual mean and the first three harmonics of climatological annual variation) and the interannual variability (time mean of the 120 days immediately preceding each day) are first removed. Following Neena et al. (2014) and Liu et al. (2017), the 14-year climatology fields of forecasted OLR and U850 as a function of both the starting date of hindcast and the lead day were first removed from the total fields to obtain the hindcast anomaly fields. Then, similar to the procedure for observations, the running mean of the preceding 120 days (the corresponding observational fields are appended before the forecast starting date) for each day was further removed to avoid the influence of interannual variability.
The real-time multivariate BSISO indices, which are the PC indices, are extracted by projecting the above observational and predicted field anomalies onto the two observational CEOF modes mentioned above. Then, both the hindcast and observed PCs are normalized by the standard deviation of the observed PCs. The amplitude is defined as (PC1 2 + PC2 2 ) 1/2 , while the phase angle is denoted as tan −1 (PC2/PC1). The initially strong and weak cases are classified based on the observed BSISO amplitude. A strong BSISO case is identified when its amplitude is larger than 1.5 (37 % of total observation cases) and a weak case for those equal to or less than 0.8 (41 % of total observation cases) during the boreal summer from 2000 to 2013.
The prediction skill is assessed by calculating the bivariate anomaly correlation coefficient (ACC), the bivariate root mean square error (RMSE), and the phase angle error (PAE) (Lin H. et al. 2008; Lee et al. (1) Here, t, τ, and N are the initial condition time, forecast lead/lag time in days, and the total number of forecasts, respectively; a 1 ( t ) and a 2 ( t ) are the observed PC1 and PC2 at time t, and b 1 (t, τ) and b 2 (t, τ) are the corresponding forecasts at time t for lead/lag time τ. A common benchmark for useful prediction skill is the lead time when the bivariate ACC is larger than 0.5 and the bivariate RMSE is less than 2 for a climatological forecast (Rashid et al. 2011). PAE is calculated to examine the propagation speed error in predictions; a negative PAE indicates that the predicted propagation is slower than the observation. Details about the assessment method are provided by Xiang et al. (2015).

Predictability and prediction skill of BSISO over EAWNP
In this section, we lay emphasis on the overall pre-dictability and prediction skill as well as their dependence on the BSISO amplitude and phase. Similar to prediction skill, the predictability is also measured by calculating the bivariate ACC between each ensemble member (regarded as "truth" based on the assumption that the model is perfect) and the ensemble mean of the other three members (regarded as "perturbations") and then averaging the ACC for all the ensemble subsamples. This indicates that the predictability is the upper limit of the prediction skill theoretically. Thus, the difference between the predictability and the prediction skill will offer room for further improvement of the BSISO prediction by reducing the model error and improving initial conditions in the current prediction system.

Overall prediction skill
In this section, the overall prediction skill is first investigated. Figure 7 gives the bivariate ACC and RMSE of the BSISO prediction over the EAWNP region for the ensemble mean and four individual ensemble members. The bivariate ACCs for both the single-member and four-member ensemble mean drop rapidly during the first 10 days, and then drop gradually on the following 10 -40 days (Fig. 7a). Taking the ACC in excess of 0.5 as a criterion, the single-member prediction skill (dashed lines) is about 10 -13 days. As expected, the ensemble mean ACC with the prediction skill of 14 days is superior to that from individual members. The skill is similar to the skill of BSISO over the Indian monsoon region in coupled models that participated in the ISVHE project (Lee et al. 2015;Lee and Wang 2016). Similarly, the bivariate RMSE displays rapid growth during the first 10 days followed by a period of slower error growth (Fig. 7b). The superiority of the ensemble mean is also prominent, particularly for lead times beyond 10 days. Hereafter, in this study we mainly focus on analysis of the four-member ensemble mean results. The ensemble spread defined by the standard deviation of ensemble member hindcasts of PC1 and PC2 relative to their corresponding ensemble mean is also examined (Fig. 7b). It is well known that the ensemble spread tends to equal the error of the ensemble mean in a perfect model forecast system (Kim et al. 2014). However, the ensemble spread shown in Fig. 7b is much smaller than the ensemble mean RMSE for all lead times. This indicates that the current prediction system is underdispersive, which needs to be overcome to further improve the BSISO prediction skill.

Dependency of predictability and prediction skill
on BSISO amplitude Previous studies have indicated that the prediction skill and predictability depend strongly on the initial or target BSISO/MJO amplitude (Fu et al. 2013;Wang et al. 2014). Figure 8a shows the prediction skill measured by bivariate ACC for forecasts initialized with weak and strong BSISO cases. The prediction skills for the initially strong cases are systematically higher than the initially weak cases from the first day of forecast, with a prediction skill up to 16 days for the strong cases versus 13 days for the weak cases. This result demonstrates the dependence of prediction skill on the initial amplitude, probably owing to the organized (disorganized) convection in the initial condition in the strong (weak) BSISO. As for the predictability, the potential skill is effective for lead times of about 24 -28 days for all cases, initially strong cases, and weak BSISO cases, and the differences among them are small (Fig. 8b). While this implies that the predictability of BSISO does not depend on the initial amplitude, it indicates a skill gap of more than 10 days that could be overcome by improving the current prediction system. The above results indicate that the model exhibits high prediction skill starting from an existent BSISO. To examine the predictability and prediction skill of BSISOs prior to their occurrence, Figs. 8c and 8d further show the bivariate ACC for forecasts targeting strong and weak BSISOs relative to forecast lag days. Here, day 0 denotes the occurring time for the target BSISO events, whereas a negative number represents the lag days before the BSISO occurrence. It is clear that both the predictability and prediction skill show stronger dependence on target BSISO amplitude, with much higher predictability and prediction skill for target strong cases than weak cases. The prediction skill is about 18 days for target strong cases and only about 6 days for target weak cases (Fig. 8c), while the predictability is about 30 -33 and 14 -18 days for target strong and weak cases, respectively (Fig. 8d). This suggests that a larger gap between the practical prediction skill and the potential predictability can be filled for the target strong cases.

Dependency of predictability and prediction skill on BSISO phase
The phase dependence is another key feature of BSISO predictability and prediction skill. To analyze such dependence, Fig. 9 shows the prediction skill of BSISO as a function of lead/lag time and initial/target phase. For each phase, only cases with the initial/target phase angle within a given phase are used. Here, each (strong) phase has about 50 (20) cases. It is shown that the differences in bivariate ACC among various initial phases become distinct beyond the lead time of about 10 days, characterized by high skill in phases 4 -7 with a peak for initial phase 5 and relatively low skill during initial phases 1 -3 (Fig. 9a). The contrast is more appreciable for initially strong cases, with skill up to 25 days for hindcasts initialized at phase 5, but only about 9 days for phase 3 (Fig. 9b). This implies that BSISO prediction is more skillful for initial convection in East Asia, the equatorial Pacific, and the Bay of Bengal, but skill is lower for prediction starting from active or active-to-break transition phases of WNP rainfall (see Fig. 6).
With regard to the variation of prediction skill for the target phase (Figs. 9c, d), relatively higher skill occurs for the forecasts targeting BSISO in phases 1 -2 (about 17 days) and 5 -6 (about 15 days), while a sharp decrease appears in phases 4 and 8 (about 11 days) for the lag time after 10 days (Fig. 9c). Similar results are found for target strong cases with prediction skill enhanced for most of the phases during the lag time of 10 -20 days except for target phase 4. Note that the skill decrease in target phase 4 corresponds to the low skill for forecasts from initial phases 1 -2, as shown in Fig. 9b. Phases 1 -2 and 5 -6 correspond to the extreme phases of WNP summer rainfall, and high skill in the forecasts targeting BSISOs in these phases is probably associated with the stronger observed BSISO amplitude over this region, as seen in Figs. 3 and 6. Figure 10 shows the predictability and differences between predictability and prediction skill in bivariate ACC. The predictability of BSISO is significantly lower in initial phases 2 -3 than in other phases beyond the lead time of about 20 days (Fig. 10a). The results are more distinct for initial strong cases (Fig. 10b). This result is consistent with the phase dependence of prediction skill shown in Figs. 9a and 9b, indicating that the model struggles to correctly forecast propagation of the rain band from WNP to the East Asian region. The differences between potential predictability and real prediction skill are apparent for all the initial phases, and show the largest values at initial phases 2 and 6 at lead times of 5 -25 days and at phase 4 as the lead time increases beyond 25 days for the initially strong cases. As for the target phase, relatively low predictability is found for target phase 3, similar to that for the initial phase (Fig. 10c). For the target strong cases (Fig. 10d), the predictability is significantly increased for all phases but with a higher predictability exceeding 35 days at target phase 3, which is opposite that for all cases as shown in Fig.  10c. This indicates that the predictability of EAWNP BSISO is more dependent on its amplitude than its target BSISO phase. The difference between predictability and prediction skill is notable for all targeting BSISO phases and is especially large at longer lag times at phases 3 -4 and 6 -7 for the target strong cases (Figs. 10c, d). This implies that the overall BSISO prediction skill can be enhanced in specific BSISO phases by decreasing the forecast error in the current prediction model.

Evaluation of BSISO amplitude and propagation
prediction To examine the amplitude and propagation features, the change of the predicted BSISO amplitude and PAE as a function of forecast lead day and initial phase for the initially strong cases are given in Fig.  11. It is found that the BSISO amplitude is always underestimated, and the mean amplitude for the first 20 days is about 22 % weaker than the observation (Fig.  11a). This is a common drawback for current operational model systems shown in many previous studies (Rashid et al. 2011;Wang et al. 2014;Neena et al. 2014;Xiang et al. 2015). For individual initial phases, the underestimated amplitude during the first 20 days is evident for all phases, especially for phases 1 and 8 (Fig. 11c), when the observed amplitude is relatively high and the BSISO is well organized. As for the PAE, positive values are found for most leading days with a mean value of 2.6° averaged over the first 20 days (Fig.  11b). This indicates that the predicted propagation speed of BSISO is faster-than-observed, which is in agreement with the model bias of propagation speed in BSISO simulation (see Fig. 5) and MJO prediction (Liu et al. 2017). Similarly, PAE varies considerably with initial phases with positive error occurring for most phases, especially for phase 6 (Fig. 11d).
To further assess the BSISO amplitude and propagation, Figure 12 shows composite anomalies of OLR and U850 averaged along 125 -135°E for both observation and prediction for initially strong BSISOs. Here, we compare the results for initial phases 1 and 5, since phase 1 is the phase where hindcasts exhibit low prediction skill and a relatively large gap between the predictability and prediction skill, while phase 5 shows the opposite (Figs. 9, 10). For phase 1, the hindcasts have some capability to predict the northward propagation signals, and the amplitude of U850 is also realistically reproduced (Figs. 12a, b). However, the positive OLR anomaly is not captured over the equatorial Pacific at lead times of 1 -10 days and almost disappears over the WNP region beyond the 10-day lead time (Fig. 12b). This discrepancy might contribute to the low prediction skill of BSISO in this phase, as shown in Fig. 9a. For phase 5, although the amplitudes of both U850 and OLR signals are still underestimated, the northward propagation of the negative OLR signal from the equatorial Pacific to WNP as well as its lead-lag relationship with U850 are better predicted (Fig. 12d). Note that a faster propagation of the U850 anomaly is predicted over the WNP region in both phases, especially at lead times of 1 -20 days. This is possibly the reason for the fast phase speed of the predicted BSISO, as shown in Figs. 11b and 11d.

Summary and discussion
BSISO, with a period of 20 -70 days, exerts a crucial role in weather and climate variability over the EAWNP region and can offer a key signal for the extended-range forecast of the East Asian summer monsoon. In this study, we analyze the predictability and prediction skill of the EAWNP BSISO in BCC_ CSM2.0, which participates in the S2S Prediction Project based on a set of hindcasts for the period 2000 -2013. The hindcast experiments were initialized with a 5-day interval, and the four-member ensemble forecast runs during the hindcast period. BSISO is diagnosed using the index taken as the PCs of the two Fig. 10. Predictability (bivariate ACC, contours) and differences between predictability and prediction skill (bivariate ACC, shadings) as a function of different initial/target phases (x axis) and forecast lead/lag days (y axis) for (a, c) all cases and (b, d) initial/target strong cases. The contour interval is 0.1.
leading CEOFs of OLR and U850 averaged along 90 -155°E. The bivariate ACC and RMSE are used to measure the prediction skill and predictability. Although BCC_CSM2.0 underestimates both the mean and intra-seasonal variability of the East Asian summer monsoon rainfall and overestimates them over the WNP region, it simulates the dominant modes and the associated northward propagation of EAWNP BSISO reasonably well. The useful BSISO prediction skill is about 14 days, while the predictability is around 24 -28 days in BCC_CSM2.0 as estimated by the bivariate ACC dropping to 0.5 and the RMSE increasing to 2 . This two-week skill is comparable to the performance of dynamical models that participated in the ISVHE project (Lee and Wang 2016). A comparison between the ensemble spread and error in hindcasts indicates that the ensemble prediction system is underdispersive (Fig. 7), implying that optimization of the ensemble prediction strategy is needed for the current model prediction.
The BSISO prediction skill and predictability obviously depend on its amplitude. Relatively higher prediction skill is found for the initial/target strong BSISO cases than for the initial/target weak BSISO cases. For the potential predictability, while a similar dependence is found for the target BSISO amplitude, it does not show any sensitivity to the initial BSISO amplitude. We also examine the prediction skill and predictability of BSISO with respect to its initial/target phase. For the initial phase, the model tends to show relatively higher skill for predictions from phases 4 -7 when the convection is initially located over East Asia, the equatorial Pacific, and the Bay of Bengal, and lower skill for phases 1 -2 when the convection develops in the South China Sea and the WNP region. Similarly, analysis on the dependence of predictability indicates that BSISO is less predictable when the BSISO-related rain band propagates northward from WNP to East Asia. As for the target phase, while no significant dependence of BSISO predictability is found, prediction is skillful when the target phase is the dry (phases 1 and 2) or wet (phases 5 and 6) phase of WNP rainfall, probably due to the larger BSISO amplitude in these phases. The difference between the predictability and prediction skill of EAWNP BSISO is distinctive for all the initial/target phases, indicating that there is still much room for further improvement in model physical processes and initialization schemes.
Some common model problems in BSISO/MJO prediction also occur in our study. For example, the predicted BSISO amplitude is always underestimated during the entire forecast period with a dramatic drop at the beginning of the prediction, which is a common deficiency in MJO/BSISO prediction (Seo et al. 2005;Rashid et al. 2011;Xiang et al. 2015). It is probable that the underestimated BSISO amplitude is partially ascribable to the model's weak BSISO variability in the free run, especially over the East Asian region as shown in Fig. 3. The sharp drop in the amplitude at the beginning of the prediction is probably associated with the initialization shocks, which result from the imbalance of the coupled system due to the utilization of separate ocean and atmosphere analyses in the initialization conditions. This needs to be more carefully examined in future work.
Moreover, many previous studies have found that the active-to-break monsoon transition associated with BSISO is more predictable than the break-toactive transition (Goswami and Xavier 2003;Fu et al. 2013). However, the results in our study, in contrast, demonstrate that the active-to-break transition of the WNP summer monsoon is less predictable. The fasterthan-observed BSISO propagation in our study is also contrary to many previous studies concerning the MJO prediction (Kim et al. 2014). This indicates that these features vary from model to model, which may be mainly associated with the insufficiency in model physics. Finally, previous studies have demonstrated that the multi-model ensemble approach produces better prediction quality as compared to any single model (Neena et al. 2014;Lee et al. 2015). Thus, to yield more insight into the model deficiency and development of the ensemble generation strategy, future work will try to assess the forecast skill of intraseasonal variability in a multi-model framework from the S2S Prediction Project.