Impacts of atmospheric and oceanic initial conditions on boreal summer intraseasonal oscillation forecast in the BCC model

In this paper, we evaluate the capability of the Beijing Climate Center Climate System Model (BCC-CSM) in simulating and forecasting the boreal summer intraseasonal oscillation (BSISO), using its simulation and sub-seasonal to seasonal (S2S) hindcast results. Results show that the model can generally simulate the spatial structure of the BSISO, but give relatively weaker strength, shorter period, and faster transition of BSISO phases when compared with the observations. This partially limits the model’s capability in forecasting the BSISO, with a useful skill of only 9 days. Two sets of hindcast experiments with improved atmospheric and atmosphere/ocean initial conditions (referred to as EXP1 and EXP2, respectively) are conducted to improve the BSISO forecast. The BSISO forecast skill is increased by 2 days with the optimization of atmospheric initial conditions only (EXP1), and is further increased by 1 day with the optimization of both atmospheric and oceanic initial conditions (EXP2). These changes lead to a final skill of 12 days, which is comparable to the skills of most models participated in the S2S Prediction Project. In EXP1 and EXP2, the BSISO forecast skills are improved for most initial phases, especially phases 1 and 2, denoting a better description for BSISO propagation from the tropical Indian Ocean to the western North Pacific. However, the skill is considerably low and insensitive to initial conditions for initial phase 6 and target phase 3, corresponding to the BSISO convection’s active-to-break transition over the western North Pacific and BSISO convection’s break-to-active transition over the tropical Indian Ocean and Maritime Continent. This prediction barrier also exists in many forecast models of the S2S Prediction Project. Our hindcast experiments with different initial conditions indicate that the remarkable model errors over the Maritime Continent and subtropical western North Pacific may largely account for the prediction barrier.


Introduction
The boreal summer intraseasonal oscillation (BSISO) and Madden-Julian oscillation (MJO) are two major modes of atmospheric intraseasonal variability over the tropics. Different from the MJO that is robust in the boreal winter with eastward propagation along the equator, the BSISO prevails in the boreal summer and often propagates northward over the tropical Indian Ocean to the western Pacific region. The BSISO is closely related to the weather and climate events in Asia or globally: the generation and reinforcement of tropical cyclones over the western Pacific (Wang and Zhou 2008;Kikuchi and Wang 2010), the outbreak and retreat of Asian summer monsoon (Hoyos and Webster 2007;Rajeevan et al. 2010), and the variability of extra-tropical circulation and precipitation (Ding and Wang 2007;Moon et al. 2013;Hsu et al. 2016). Therefore, the BSISO is an important target in climate modeling and forecasting.
The simulation of the BSISO, however, bears various problems, e.g., the unrealistic BSISO structure, propagation, and intensity due to the misrepresented essential model physics and air-sea coupling (Waliser et al. 2003;Fu and Wang 2004;Kim et al. 2008;Lin et al. 2008bLin et al. , 2011. The state-ofthe-art climate models are slowly and steadily improved in capturing the basic characteristics of the BSISO (Sabeerali et al. 2013;Fang et al. 2017), which leads to steady progresses in dynamic forecasts of the BSISO to some extent. It is found that the intraseasonal oscillation over the western Pacific monsoon region and Indian monsoon region can be predicted with a lead time of 1-2 weeks in several climate models (Lin 2013;Suhas et al. 2013). Particularly, the European Centre for Medium-Range Weather Forecast (ECMWF) model can provide a useful skill of nearly 20 days for BSISO prediction (Lee et al. 2017). Moreover, the forecast skill can be further enhanced by multi-model ensemble forecasts Neena et al. 2014;Lee et al. 2015). Although tangible progresses have been made, the intraseasonal forecasting by dynamic models is still challenging due to its remarkable sensitivity to various model settings, such as model physics, initial conditions, and air-sea coupling (Seo et al. 2005;Fu et al. 2008Fu et al. , 2011Ling et al. 2014;Jie et al. 2017).
The prediction skill and predictability of the BSISO depend on the initial state of various climate components and their persistence. As a phenomenon with the time scale between synoptic weather and seasonal climate, the BSISO is remarkably affected by both the atmospheric initial conditions and underlying surface forcing. By assessing the sensitivity of BSISO forecast to atmospheric initial conditions, it is found that the forecast skills of the BSISO vary when initialized with different atmospheric reanalysis data (Fu et al. 2009. Previous study revealed that the tropical intraseasonal oscillation is sensitive to different sea surface temperature (SST) forcing, especially to that varying on the intraseasonal time scale (Pegion and Kirtman 2008). Moreover, lower boundary SST state is important for the structure and propagation of the BSISO (Liu et al. 2015). Accurate SST is thought to be essential for maintaining a higher predictability of the BSISO, and skills of the ocean-atmosphere coupled forecasts and uncoupled forecasts may be comparable under the same SST conditions (Fu et al. 2008). It is also indicated that the SST initial conditions can exert more influence on the predictability of tropical intraseasonal oscillation than the atmospheric initial conditions .
In recent years, the Sub-seasonal to Seasonal (S2S) Prediction Project, which provides S2S hindcasts (up to 60 days) and real-time forecast products from 11 operational and research centers, has been jointly established by the World Weather Research Program (WWRP) and the World Climate Research Program (WCRP) . One of the key foci of the S2S Prediction Project is the predictability of intraseasonal variability. A series of studies have examined the performances of multiple S2S forecast models in predicting MJO, BSISO, and corresponding teleconnected influences, demonstrating remarkable differences among various models Vitart 2017;Zhou et al. 2018). Results have partially shown that the Beijing Climate Center (BCC) S2S forecast model, which participates in the S2S Prediction Project, exhibits relatively low skill in MJO and BSISO forecasts. Furthermore, the BCC S2S model's limited forecast skill of the BSISO is thought to be associated with its inability to describe the physics for northward propagation of BSISO at relatively longer lead time (He et al. 2019). To reduce the deficiency in S2S forecasts, Liu et al. (2017) improved model initial conditions and enhanced MJO forecast skill by using the BCC S2S model. Therefore, it is worthwhile to investigate whether or not and to what extent the BSISO forecast can be improved by the optimization of initial conditions, and to further explore the relationship between BSISO forecasts and the reliability of model's simulation and initialization. It is also necessary to verify the validity of the optimization of initial conditions in improving S2S forecasts, and to understand different impacts of model initial conditions and internal physics on BSISO forecasts.
The rest of this paper is organized as follows. Details of the model, experiments, validation data, and methods are given in Section 2. In Section 3, we present the BSISO simulation by the BCC S2S model, and in Section 4, we show the BSISO prediction skill of the model and its sensitivity to initial conditions. Conclusions are provided in Section 5.

Model and data
The model used in this study is the BCC Climate System Model (BCC-CSM) version 1.2, which participates in the S2S Prediction Project ). In the fully coupled BCC-CSM, the atmospheric component is the BCC Atmospheric General Circulation Model (AGCM) version 2, with a T106 triangular truncation in the horizontal direction and 40 layers in the vertical direction. The land component is the BCC Atmosphere and Vegetation Interaction Model version 1.0. The ocean and sea ice components are the National Oceanic and Atmospheric Administration (NOAA) Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model version 4 (MOM4) and Sea Ice Simulator, respectively. The BCC model has been used in short-term climate operational prediction with reasonable skills in forecasting major climate variability at seasonal and sub-seasonal time scales (Liu et al. 2015a(Liu et al. , 2018. The hindcast dataset of the BCC model for the S2S Prediction Project is used in this study, and the hindcasts started on 1st, 6th, 11th, 16th, 21st, and 26th of each month during 2000-2013 are selected. Each hindcast makes a 60-day forecast integration and includes four ensemble members, with a 6-h time lag of atmospheric initial conditions for two neighboring members. For the S2S hindcasts, the atmospheric initial conditions (i.e., surface pressure, air temperature, zonal and meridional winds) are from the NCEP Reanalysis 1 (NCEP-R1; Kalnay et al. 1996), and the oceanic initial conditions (i.e., temperature and salinity) are from the multi-level temperature and salinity fields of the BCC Global Ocean Data Assimilation System (BCC-GODAS). A long-term coupled initialization run is conducted to assimilate the above observations by using nudging strategy, which produces the model initial conditions for all hindcasts. More details about the BCC S2S hindcast experiments can be found in Liu et al. (2017).
Two additional sets of hindcast experiments are conducted to improve model's sub-seasonal forecast skill. Experiment I (EXP1) uses the similar schemes as the S2S hindcast experiments, except that the atmospheric initial conditions are improved by introducing more reliable atmospheric reanalysis data, i.e., the National Centers for Environmental Prediction (NCEP) FNL (final) Operational Model Global Tropospheric Analyses to replace the NCEP-R1. On the basis of EXP1, Experiment II (EXP2) introduces the daily Optimum Interpolation Sea Surface Temperature version 2 (OISST; Reynolds et al. 2007) to upgrade the SST initial conditions; in other words, EXP2 uses improved atmospheric and SST initial conditions. From sea surface to the depth of 30 m, a linear transition of initial fields from OISST to BCC-GODAS is adopted to keep the vertical continuity of the ocean state, and to amplify the downward impact of reliable SST initial conditions. With improved initial conditions, both EXP1 and EXP2 use the same hindcast dates as the S2S hindcast experiments. In addition, a 15-year free run without the assimilation of observations is conducted to evaluate the performance of the model itself in simulating the BSISO.
For model evaluation, we use daily outgoing longwave radiation (OLR) and wind fields from the ERA-Interim reanalysis (Dee et al. 2011), and daily SST field from the NOAA OISST (Reynolds et al. 2007). The observed daily anomalies of OLR and 850-hPa zonal wind (U850) are obtained by removing the slow annual cycle (annual mean and the first three annual harmonics of the 15-year average) and interannual variability (average of the previous 120 days). This processing method is mostly used for the verification of real-time intraseasonal forecast, although it may be limited in filtering out the interannual and long-term variability (Lyu et al. 2019a). For the results from the S2S hindcast, EXP1, and EXP2, the daily anomalies are defined by removing the lead-dependent model climatology and the mean of anomalies over the previous 120 days. To better describe the BSISO features as in Lin (2013), multivariate Empirical Orthogonal Functions (EOFs) are calculated for the observed anomalies of OLR and U850 averaged over 90°-150°E from May 1 to September 30, and the BSISO index is then defined by the two leading PCs (PC1 and PC2). The amplitude and phase angle of the BSISO are expressed as (PC1 2 + PC2 2 ) 1/2 and tan −1 (PC2/PC1), respectively. The forecasted PCs are further computed by projecting the forecasted intraseasonal anomalies onto the observed EOFs. Following the definition in Lin et al. (2008a), the bivariate anomaly correlation (BAC) and root mean square error (RMSE) between the forecasted and observed PCs are computed to measure the BSISO forecast skill.

BSISO characteristics in the model simulation
Given that the capability of BCC-CSM in simulating the key characteristics of BSISO is highly influential to its performance in predicting the BSISO, an uninitialized 15-year free run is conducted to show the representation of the BSISO in the model control simulation. Figure 1 presents the two leading EOF modes of daily OLR and U850 anomalies averaged over 90°-150°E. In the observation, the first EOF mode shows a strong positive center of convection anomaly near 15°N, and the second mode exhibits suppression and enhancement of convection near 10°N and 25°N, respectively. In these two modes, the maximum westerly (easterly) wind anomalies are about 5°to the south of the maximum (minimum) convection anomalies. The simulated EOF modes agree well with the observation, exhibiting a pattern correlation coefficient of 0.95 (0.91) between the observed and simulated EOF1 (EOF2) mode. The first two leading modes explain about 47% and 38% of the total variance in the observation and simulation, respectively. The results suggest that Power spectral features of the first two PCs in the observation and simulation are shown in Fig. 2. The observed PC1 and PC2 are characterized by a typical intraseasonal period, showing a dominant period of 30-70 days. In the simulation, the power spectrum of PC1 has no significant peak at the intraseasonal time scale, while that of PC2 shows two bumps near the 40-day and 20-day time scales, respectively. Compared with the observation, the spectral variances are clearly underestimated at 40-100-day time scale, but are somewhat overestimated at 20-30-day time scale in the simulation, indicating that the northward-propagating BSISO is featured by a weaker strength and faster speed in the model. Figure 3 displays the lag correlation between PC1 and PC2. In the observation, maximum (minimum) correlation between corresponding PCs is found when PC1 leads (lags) PC2 by 10 days, which are similar to the results of Lin (2013). In the model, when PC1 leads PC2, the positive correlation is basically consistent with that in the observation. However, when PC1 lags PC2, the negative correlation is obviously different from the observation, showing a lower correlation and a quicker transition with lag time. It suggests that the simulated BSISO propagates faster and shows an unreasonable phase evolution when transitioning from EOF2 to EOF1, which partially corresponds to the overestimated variance at high frequency (Fig. 2).
Composites of the OLR and 850-hPa wind anomalies for different BSISO phases are shown in Fig. 4, only BSISO events with larger-than-one amplitude being selected. In the observation, the convection initiates from the tropical Indian Ocean and Maritime Continent (phase 1), and propagates northeastward to the South China Sea as well as the tropical western Pacific (phases 2-4). With further northward propagation, the convection signal finally reaches southeastern China and the subtropical western North Pacific (phases 7-8). The model can basically reproduce the northwardpropagating convection from the Indian Ocean to southeastern China, but with some apparent deficiencies. Particularly, the simulated wet anomalies at phases 1-5 and dry anomalies at phases 6-8 are considerably weaker than the observation. The underestimate is especially remarkable over the tropical Indian Ocean and Maritime Continent. In addition, from phase 6 to phase 8, a faster-than-observed propagation of convection anomaly over southeastern China and the western North Pacific and a more northward extension of convection anomaly are found, which generally agree with the results shown in Fig. 3. The weaker-than-observed BSISO variability over the tropics, especially over the equatorial Indian Ocean, is also a common deficiency for state-of-the-art climate models, including most models participated in the Coupled Model Intercomparison Project Phase 5 (Sabeerali et al. 2013). Fig. 2 Power spectra of the first two PCs of the combined EOF analysis for observations (solid red) and control simulation (solid blue). a PC1; b PC2. The rednoise spectrum at 95% confidence level is shown by dashed line

BSISO forecasts and impacts of improved initial conditions
Several sets of hindcast experiments are conducted using the BCC S2S forecast model. The hindcasts for the S2S Prediction Project suffer severe limitations due to the problems of imperfect model and unreliable initial conditions ). Thus, EXP1 and EXP2 are carried out to improve the model forecast skill by upgrading the atmospheric and both atmospheric and SST initial conditions, respectively. Figure 5 shows the forecast skills of SST and U850 in the first 2 weeks of forecasts, and the differences of skill among different hindcasts. In the first week, with optimized atmospheric and ocean initial conditions in the model, remarkable forecast skills of SST and U850 are found over most areas in EXP2 (Fig. 5a1, b1). In the second week, the forecast skill of SST remains high over the extra-tropics but relatively lower over the tropical Indian Ocean and western Pacific (Fig. 5a5), while that of U850 is considerably low over most areas (Fig. 5b5), denoting an obviously higher predictability of SST than that of atmospheric circulations. During the first 2 weeks of forecasts, SST forecast skill is little increased over the tropics from the S2S hindcasts to EXP1 (Fig. 5a2, a6), while U850 forecast skill is enhanced over many regions, especially over the southeastern Indian Ocean and southwestern Pacific (Fig. 5b2, b6). Moreover, during the first 2 weeks, the SST forecast skill is remarkably improved over most areas in EXP2 over EXP1 (Fig. 5a4, a8), while the U850 forecast Fig. 5 Anomaly correlations between observations and forecasts for a1, a5 SST and b1, b5 850-hPa zonal wind averaged during the first and second weeks of forecast. Also shown are the differences of correlation skills a2, a6, b2, b6 between EXP1 and S2S hindcasts; a3, a7, b3, b7 between EXP2 and S2S hindcasts; and a4, a8, b4, b8 between EXP2 and EXP1. Stippling represents the statistical significance of correlation above the 99% confidence level skill basically remains unchanged in the first week, but is slightly improved over the tropical western Pacific in the second week (Fig. 5b4, b8). It indicates that for climate variability forecasts over the western Pacific, the improvement of atmospheric initial conditions exerts more impacts on atmospheric circulation forecasts rather than on SST forecasts. The optimization of SST initial conditions can impact SST forecasts, and meanwhile affect the atmospheric circulation forecasts beyond the lead time of 1 week. Figure 6 shows the overall forecast skills of the BSISO in various hindcast experiments. The BAC and RMSE between forecasted and observed PCs are computed to measure the BSISO forecast skill. Taken BAC = 0.5 as the threshold of useful skill, the forecast skill of about 9 days is found in the S2S hindcasts, and the skill increases to 11 and 12 days in EXP1 and EXP2, respectively (Fig. 6a). Also, taken the RMSE = 1.414 as the upper limit of useful skill, similar improvements are seen as those of BAC (Fig. 6b). The forecast skills in the experiments with improved atmospheric and oceanic initial conditions are comparable to those in several state-of-the-art models (Lee and Wang 2016;Jie et al. 2017). These results reveal a 3-day increase of forecast skill from the S2S hindcasts to EXP2, demonstrating the impacts of optimized atmospheric and oceanic initial conditions. It confirms the importance of accurate initial conditions in BSISO forecasts revealed by previous studies (Fu et al. 2008(Fu et al. , 2009Lee et al. 2015). In addition, from EXP1 to EXP2, the BSISO forecast skill remains almost unchanged in the first week but is improved in the second week. This corresponds to the skill difference of U850 over the tropical western Pacific between the two sets of hindcasts as shown in Fig. 5, indicating the gradually emerging impacts of ocean initial conditions beyond the lead time of 1 week.
The sensitivity of forecast skill to initial conditions may vary with different definitions of the BSISO. The BSISO index defined by Lin (2013) is used in this study, namely, the zonal average of OLR and U850 over 90°-150°E is examined. By using another BSISO index ), the  Results show that the forecast skill is about 9 days in the S2S hindcasts, which increases to 10 and 12 days in EXP1 and EXP2, respectively (figures omitted). Although the former index, which shows a smoother spatial feature, is more easily to be depicted by the model than the latter index, comparable skills are found for both indexes. Comparison with the results shown in Jie et al. (2017) indicates that the forecast skill of EXP2 is comparable to those of the models in the S2S Prediction Project. Figure 7 shows the interannual variation of BSISO forecast skill in the first 2 weeks of forecasts during 2000-2013. From the S2S hindcasts to EXP1 and EXP2, the overall forecast skill of the BSISO shows a slight improvement in both the first and second weeks, indicating the positive effects of the optimization of initial conditions. Meanwhile, the forecast skill exhibits significant interannual variation. For example, the forecast skill was extremely low in 2010, while it was relatively higher in 2003 and 2004 than that in the other years. From the S2S hindcasts to EXP1 and EXP2, the features of monthly evolution of forecast skill (Fig. 8) reveal that the useful skill increases to 2-3 weeks during May-August 2003 several times. However, relatively low skills remain unchanged in most summers, and especially low values of less than 5 days are found during May-July 2010. The failure of initialization optimization in 2010 may be caused by the obviously low predictability of the BSISO itself, in addition to the uncertainty in the interaction between model errors and initial errors. The above skill differences might be closely related to the interannual variation of BSISO amplitude (Liu et al. 2015b), given that the forecasts by various models are always more skillful for BSISO events with stronger amplitude than those with weaker amplitude (Lee et al. 2015;Jie et al. 2017). In this study, we find that clear interannual differences are also featured by different skill degrees for the three hindcast sets, and different amplitude of skill declines from the first week to the second week in those years, altogether revealing complicated impacts of atmospheric and oceanic initializations. Figure 9 illustrates the dependence of BSISO forecast skill on initial/target phases and lead times in various hindcasts. We again take BAC = 0.5 as the threshold of useful skill. In the S2S hindcasts, the forecast skills range from 8 to 12 days for different initial phases, with the minimum for initial phase 6 and relatively high skills near initial phases 8 and 1 (Fig. 9a). It indicates that the BSISO is easy to predict when it initiates from the equatorial Indian Ocean, but difficult to forecast when it starts from the South China Sea and western North Pacific. From the S2S hindcasts to EXP1, the forecast skills are obviously improved for many initial phases except phases 4-6 (Fig. 9b). Moreover, the forecast skills from initial phases 1 and 2 are further improved in EXP2, with the maximum skill near the 18-day lead time. However, useful skills at the lead time less than 10 days for forecasts near initial phase 6 are hardly changed (Fig. 9c). Besides, the forecast skill varies with target phase, with relatively higher values near target phases 5-8 and lower values near target phases 3-4, corresponding to skill peaks near initial phases 8, 1, and 2, and to skill valley near initial phase 6 ( Fig. 9e, f). The results indicate that the BSISO is difficult to forecast when the convection signal dissipates over the western North Pacific and reinitiates from the Indian Ocean and Maritime Continent, but can be skillfully predicted when the convection signal propagates from the Indian Ocean into the western North Pacific in both EXP1 and EXP2. Previous study has shown that the most significant forecast skills of the BSISO near initial phases 1 and 8 are found in two thirds of the models in the Intraseasonal Variability Hindcast Experiment (ISVHE) Project (Lee et al. 2015), whereas the lowest forecast skills of the BSISO mostly appear around target phases 3-5 for 10 models in the S2S Project . Although the forecast skill and its sensitivity to initial and target phases are always model-dependent, these results to some extent suggest a common prediction barrier of the BSISO as depicted in many of the stateof-the-art climate forecast models.
Given the most significant skill enhancement for initial phase 2 and the smallest improvement for initial phase 6, Fig. 10 shows the time-latitude composites of forecasted OLR and U850 anomalies for these two initial phases. For initial phase 2, the enhanced convection anomaly in the observation propagates northward from the equator to the South China Sea, and the suppressed convection anomaly propagates from the South China Sea to the Northwest Pacific (Fig. 10a). In the S2S hindcasts, the intensity of convection anomaly is apparently underestimated, and the northward propagation of convection shows a faster-than-observed decay with lead time (Fig. 10c). The intensity and propagation of convection anomalies in EXP1 and EXP2 are obviously improved due to the improvements of initial conditions (Fig. 10e, g). For initial phase 6, the negative convection anomaly in the observation propagates from the equator to the South China Sea, while the positive convection anomaly propagates over the South China Sea (Fig. 10b). In the S2S Fig. 9 BSISO forecast skill (bivariate anomaly correlation) as a function of lead time and BSISO phase for strong initial or target cases (amplitude > 1.0) in a, d S2S hindcasts; b, e EXP1; and c, f EXP2. The forecast skill is stratified by initial phase (top panels) and target phase (bottom panels) Fig. 10 Time-latitude composites of outgoing longwave radiation (shadings; units: W/m 2 ) and 850-hPa zonal wind (contours; units: m/s) anomalies for the BSISO initial phase 2 (left column) and phase 6 (right column) in a, b observations; c, d S2S hindcasts; e, f EXP1; and g, h EXP2 Fig. 11 Time-latitude diagram of predictability error of outgoing longwave radiation (contours; units: W/m 2 ) and 850-hPa zonal wind (shadings; units: m/s) anomalies for the BSISO initial phase 2 (left column) and phase 6 (right column) in a, b S2S hindcasts; c, d EXP1; and e, f EXP2 hindcasts, the intensity of convection anomaly is apparently weaker than the observation, and the convection enhancement shows a longer duration over the South China Sea with a more northward propagation over the subtropical areas (Fig. 10d). In EXP1 and EXP2 (Fig. 10f, h), although the magnitude of convection enhancement is improved at short lead time to certain extent, the abovementioned deficiencies generally maintain, indicating that improved initial conditions have no obviously positive effect on the forecast of BSISO propagation from initial phase 6. It is suggested that the intrinsic model problem to sustain dry phases near the equator, as shown in Fig. 4, might have contributed to the unskillful forecasts of the BSISO from phase 6 significantly.
For the two initial phases mentioned above, the variations of predictability errors and forecast errors of U850 and OLR with lead time are shown in Figs. 11 and 12. The forecast error is defined as the RMSE between forecasts and observations, and the predictability error is estimated by computing the RMSE between each perturbed member and the control member and then averaging the results of the three ensemble subsamples. These definitions are similar to the forecast skill and predictability measured by BAC in Kim et al. (2014). Based on a perfect-model assumption, the predictability error is used to measure the error features caused by the uncertainty of initial conditions, and its difference with the forecast error can partially represent the impacts of model errors. For predictability errors of the forecasts (Fig. 11), the differences among ensemble members are always small at short lead time, and gradually increase with lead time, thus showing large spreads at longer lead time. The increase of predictability errors due to the uncertainty of initial conditions is obviously faster in the S2S hindcasts than those in EXP1 and EXP2 within 10-day lead time, denoting more reliable initial conditions adopted by the latter two. The predictability errors become remarkable near 15°N and 40°N beyond the lead time of 10-15 days, whereas the error magnitude is generally comparable among various forecasts, suggesting that an upper limit of potential predictability may be reached. The variations of forecast errors show different features (Fig. 12). In EXP2, the forecasts from initial phases 2 and 6 are both featured by relatively remarkable forecast errors over the tropical and subtropical western Pacific beyond 1-week lead time, with a maximum center near 10-20°N at 10-20-day lead time (Fig. 12a, b). Especially, the error magnitude for forecasts from initial phase 2 is obviously smaller than that from initial phase 6. The differences between forecast errors and predictability errors (Fig. 12c, d) further indicate that associated with the northward propagation of the BSISO, the model error evolves with lead time and latitude, presenting a much larger amplitude and longer duration for forecasts from initial phase 6 than from initial phase 2. Particularly, two maximum error centers are located over the Maritime Continent and the subtropical western North Pacific, respectively. These results agree with those shown in Figs. 1, 2, 3, and 4, in which a significantly underestimated convection near the Maritime Continent and a clearly faster-thanobserved propagation of the BSISO with a gradual subsiding convection anomaly over the western North Pacific are simulated. All these factors lead to a rapid decline of forecast skill for the propagation of the BSISO beyond the 10-day lead time from initial phase 6. It demonstrates that the model deficiency other than unreliable initial conditions should largely account for the abovementioned BSISO prediction barrier.

Conclusions
In this study, we evaluated the capability of the BCC model in simulating and forecasting the BSISO. We also explored the impacts of improved initial conditions on BSISO prediction skill by the model. Our main findings are summarized here.
The model simulates a faster-than-observed propagation of the BSISO from the equatorial Indian Ocean to southeastern China and the subtropical western North Pacific, with an apparently underestimate in the intensity of convection over the Indian Ocean and western equatorial Pacific. The useful forecast skill of the BSISO is about 9 days in the S2S hindcasts, which is increased to 12 days in EXP2 due to the improved atmospheric and oceanic initial conditions. The improved Fig. 12 Time-latitude diagram of a, b forecast error and c, d differences between forecast error and predictability error of outgoing longwave radiation anomalies (contours; units: W/m 2 ) and 850-hPa zonal wind anomalies (shadings; units: m/s) for the BSISO initial phase 2 (left column) and phase 6 (right column) in EXP2 forecast skill in EXP2 is comparable to the skills in several state-of-the-art dynamic forecast models (Lee and Wang 2016;Jie et al. 2017). Note that from EXP1 to EXP2, the BSISO forecast skill remains almost unchanged in the first week of forecast, but is slightly improved in the second week, indicating the gradually emerging influence of the ocean initial conditions beyond the lead time of 1 week.
The BSISO forecast skill significantly depends on initial and target phases. In the S2S hindcasts, the model shows higher skills in initial phases 8 and 1 when the convection anomaly develops over the tropical Indian Ocean, and shows lower skill in initial phase 6 when the convection anomaly is initially located over the South China Sea. From the S2S hindcasts to EXP1 and EXP2, the overall forecast skill is enhanced for most initial phases except phase 6 due to the improved initial conditions. In EXP1 and EXP2, the propagation from the tropical Indian Ocean and Maritime Continent to the South China Sea is obviously improved, but the propagation from southern China and the subtropical western North Pacific is still difficult to predict. Further exploration reveals that the rapid increase of model errors over the Maritime Continent and the subtropical western North Pacific near initial phase 6 is responsible for the low forecast skill of the BSISO and its insensitivity to the change of initial conditions.
The above results reveal that the BSISO forecast skill can be improved with the introduction of more accurate atmospheric and oceanic initial conditions. However, there are still some obvious deficiencies in the BSISO prediction due to the imperfection of the model itself. On the one hand, the fasterthan-observed propagation is predicted in various hindcasts by the model, which is attributed to the short periodicity of the BSISO in the model as shown in Fig. 2. Given that the slowerthan-observed northward propagation is also found in other models, it is suggested that this feature is highly modeldependent Fu et al. 2013). On the other hand, the model shows apparent difficulties in predicting the BSISO propagation from southern China and the subtropical western North Pacific, which also occurs in most models in the S2S Project . In this case, the convection over the western North Pacific is in an active-to-break transition, whereas the convection from the tropical Indian Ocean and the South Asian monsoon region to the Maritime Continent is in a break-to-active transition. Several studies have shown that the break-to-active transition is less predictable than the active-to-break transition of monsoon over the South Asia region (Goswami and Xavier 2003;Fu et al. 2013), which is consistent with the finding of this study. Therefore, for the prediction barrier of the BSISO during its northward propagation, which is a common deficiency in many climate models, more efforts should be made in the future for further improving the forecast skill of the BSISO.
In addition to the BSISO, the MJO is also one of the targets for improving the S2S forecast in this study, given that the MJO is a dominant mode of intraseasonal atmospheric variability and exerts remarkable modulations on the tropical and extra-tropical weather and climate events (e.g., Kikuchi and Wang 2010;He et al. 2011;Moon et al. 2013;Hsu et al. 2016;Lyu et al. 2019b). Liu et al. (2017) showed that the useful skill of MJO forecast by the BCC model was about 16 days in the S2S hindcasts, and increased to 18 days by optimizing atmospheric initial conditions and 22 days by optimizing both atmospheric and oceanic initial conditions. This indicates that the BSISO and MJO can both be well predicted by the model with updated initialization schemes, although these two climate phenomena differ in active period and propagation feature. At the second phase of the S2S Project, a BCC model version with enhanced resolution and upgraded initialization schemes will be used to further improve sub-seasonal forecast. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.