Journal of the Meteorological Society of Japan. Ser. II
Online ISSN : 2186-9057
Print ISSN : 0026-1165
ISSN-L : 0026-1165
Article: Special Edition on DYAMOND: The DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains
Prediction Skill of the Boreal Summer Intra-Seasonal Oscillation in Global Non-hydrostatic Atmospheric Model Simulations with Explicit Cloud Microphysics
Ryosuke SHIBUYAMasuo NAKANOChihiro KODAMATomoe NASUNOKazuyoshi KIKUCHIMasaki SATOHHiroaki MIURATomoki MIYAKAWA
Author information
JOURNAL OPEN ACCESS FULL-TEXT HTML
Supplementary material

2021 Volume 99 Issue 4 Pages 973-992

Details
Abstract

In this study, we assess the prediction skill of the Boreal Summer Intra-Seasonal Oscillation (BSISO) mode of one-month simulations using a global Non-hydrostatic Icosahedral Atmospheric Model (NICAM) with explicit cloud microphysics and a grid spacing of 14 km. The simulations were run as a series of hindcast experiments every day of August from 2000 to 2014. A total of 465 simulations were run with a 13950-day integration. Using forecast skill scores for statistical measurements, it was found that the model showed an overall BSISO prediction skill of approximately 24 days. The prediction skill tended to be slightly higher (∼ 2 days) when BSISO events began in the initial phases 7 to 1, which corresponded to the re-initiation phase of the BSISO, during which a major convective center over the Philippine Sea decayed and a new convective envelope began aggregating over the western Indian Ocean. The phase speed and the evolution of the amplitude of the BSISO were well simulated by the model with a clear northwest-southeastward tilted outgoing longwave radiation (OLR) structure over the Maritime Continent and the western Pacific. However, the propagation speed was slower during phases 6 and 7, and the amplitude of the BSISO largely decayed during phases 8–1, which is likely to have been associated with the stagnant behavior of the convective cells over the Philippines sea. Based on a regression coefficient analysis using the moist static energy, the stagnation of the propagation over the Philippines was found to be largely attributable to the small southerly background bias in the model over the Philippines. The bias in the large-scale circulation was likely to have been associated with the bias in the moisture field and the background monsoonal circulation. We concluded that the model physics controlling the background fields are important factors in improving BSISO prediction skill.

1. Introduction

In the Northern Hemisphere, pronounced tropical intra-seasonal oscillation (ISO) has significant impacts on the weather and climate in both the tropics and midlatitudes. During boreal summer, the ISO modulates active/break phases of the Asian summer monsoon (Krishnamurti et al. 1992; Ding and Wang 2009) and occasionally produces anomalous precipitation events in the midlatitudes through teleconnection mechanisms (Ding and Wang 2007; Moon et al. 2011, 2013). The summer ISO is conventionally referred to as the boreal summer ISO (BSISO; recently reviewed by Kikuchi 2021). The BSISO mode is often characterized by a northwest–southeastward tilted cloud convective band that initially spans from the north of the Arabian Sea to the east of the Philippines (Wang and Xie 1997) with both northward/northwestward and eastward propagation of the convective envelope (e.g., Yasunari 1979; Chen and Murakami 1988; Wang and Rui 1990; Hsu et al. 2004). These characteristics are clearly distinct from those of the boreal winter ISO mode, called Madden-Julian Oscillation (MJO), which involves the dominant eastward propagation of a convective envelope along the equator. The complex behavior of the BSISO may be attributable to asymmetric large-scale circulation — with respect to the equator, associated with the Asian summer monsoon in the Northern Hemisphere.

The ISO extensively influences local/remote weather patterns, and therefore, forecasting skills for subseasonal to seasonal predictions that range from two weeks to a few months have been intensively investigated by the research community. One recent attempt is the Sub-seasonal to Seasonal (S2S) Prediction Project (Vitart et al. 2007). Although MJO forecasts have been extensively studied, very limited research has been conducted to examine the forecast skill of the BSISO using operational and/or research models. Only a few studies have examined the forecast skills of the BSISO using operational and/or research models (Waliser et al. 2003; Fu et al. 2013; Lee et al. 2015) or the S2S Prediction Project database (Jie et al. 2017; Wang et al. 2018). Further, the prediction skill of the BSISO tends to be lower than the MJO by approximately 5–10 days (Lee et al. 2015; Wang et al. 2018), although it remains unclear why BSISO prediction is more difficult than that of MJO.

Global non-hydrostatic atmospheric models that improve the realism of convection behavior in the tropics and midlatitudes (Stevens et al. 2019) are expected to exhibit better model performance on sub-seasonal time scales. One such model, the Non-hydrostatic Icosahedral Atmospheric Model (NICAM; Tomita and Satoh 2004; Satoh et al. 2008, 2014) can accurately simulate the fundamental features of the ISO, including both MJO (Miura et al. 2007; Miyakawa et al. 2014) and the BSISO (Oouchi et al. 2009; Taniguchi et al. 2010; Yanase et al. 2010; Satoh et al. 2012; Kinter et al. 2013; Nakano et al. 2015), in short-term simulations. In the first long-term experiment (30 years) by NICAM following the atmospheric model intercomparison project (AMIP) protocol proposed by Kodama et al. (2015), it was shown that the fundamental features of the BSISO, including the northwest-southeastward tilted convective band, are well captured, albeit with a weaker amplitude by a factor of ∼ 2 and a relatively slow propagation speed compared with observations (Kikuchi et al. 2017). On an S2S time scale, Miyakawa et al. (2014, 2017) demonstrated that NICAM exhibited accurate prediction up to 26–28 days using the ensemble hindcasts of boreal winter MJO events. Recently, Nakano and Kikuchi (2019) revealed that NICAM is one of the best models for reproducing the climatological seasonal cycle of the ISOs among general circulation models (GCMs) in CMIP5. However, a statistical analysis of the prediction skill of the BSISO by NICAM is yet to be performed.

The main purpose of this study is to provide a quantitative estimate of the prediction skill of the BSISO using NICAM as a first attempt at statistically assessing the BSISO prediction skill using a global cloud system resolving model. To our knowledge, the statistical prediction skill of the BSISO by a global cloud system resolving model (including NICAM) has not been rigorously examined thus far (i.e., by conducting a series of hindcast experiments). Using new metrics based on a phase space that consists of principal components, we also discuss the processes preventing higher prediction for the BSISO in NICAM and how the predictive skill is related to model biases. This may provide valuable insight into how well such a global high-resolution model can successfully predict the BSISO. This article is organized as follows. Model configuration and methodology are described in Section 2. Model results are presented in Section 3. A discussion is provided in Section 4, and Section 5 summarizes the results and provides concluding remarks.

2. Model and analysis methodology

2.1 Experimental design and data

The 2012 version of NICAM (NICAM.12; Satoh et al. 2014) with a horizontal grid resolution of approximately 14 km was used for this study. This is the same model setting used by Nakano et al. (2015). The number for the vertical layer was 38 with a model top of 36.7 km. Moist processes were calculated by the NICAM Single-Moment Water 6 cloud microphysics scheme (NSW6; Tomita 2008) without convective parameterization. The radiative transfer processes for long and short waves were calculated using the radiation package termed mstrnX (Sekiguchi and Nakajima 2008), and the planetary boundary layer processes were calculated using the Mellor–Yamada–Nakanishi–Niino level 2 scheme (MYNN2; Nakanishi and Niino 2004; Noda et al. 2010). The initial atmospheric conditions were derived from the European Centre for Medium-Range Weather Forecasts Interim reanalysis (ERA-Interim; Dee et al. 2011). The initial sea surface temperature (SST) was also sourced from ERA-Interim data. The temporal evolution of the SST was forecasted using a slab ocean model, which was nudged with the daily climatology of SST with a persistent anomaly at the initial time with an e-folding time of seven days. The simulation began at 0000 UTC each day in August from 2000 to 2014 with an integration period of 30 days, providing 465 sets of 30-day timeintegrated data with different initial conditions. Using the same numerical settings, Nakano et al. (2015) demonstrated that the generation of Typhoon Songda (2004) could be successfully forecasted approximately two weeks prior to the event through better reproduction of the evolution of the BSISO and the extension of the monsoon trough.

ERA-Interim data were also used in the analysis of the observed atmospheric fields and outgoing longwave radiation (OLR) at the top of the atmosphere distributed at a 1° spatial and daily temporal resolution (the “SYN1deg-Day product”, Ed4A) using the Clouds and the Earth's Radiant Energy System (CERES) onboard the Aqua and Terra satellites (Doelling et al. 2013). Although we did not intend to comprehensively compare the performance of NICAM with those of other models, for reference, we also used the operational products of medium-range ensemble prediction systems created by the European Centre for Medium-Range Weather Forecasts (ECMWF; Palmer et al. 1993; Buizza et al. 2007), the UK Met Office (UKMO), and the Japanese Meteorological Agency (JMA). These data were obtained from the THORPEX Interactive Grand Global Ensemble (TIGGE) project (Park et al. 2008). The initial dates of the forecast data used were the same as those used in the NICAM simulations (2007–2012). Overall, 15 days of forecasts were simulated for the ECMWF and UKMO, and 9 days for the JMA. It has previously been reported that the ECMWF tends to exhibit the highest prediction skill for the BSISO (Lee et al. 2015; Jie et al. 2017; Wang et al. 2019).

2.2 Calculation of the BSISO index

This section introduces a bimodal ISO index (Kikuchi et al. 2012; Kikuchi 2020) and describes the spatial and temporal variability of the ISO modes. The BSISO mode inherently involves the complex propagation modes directed eastward, northward and northeastward (Lawrence and Webster 2002), and therefore, Eigen techniques represent a common and convenient method to describe the sequential spatial-temporal evolution of the ISOs. The ISO index presented by Kikuchi et al. (2012) is a bimodal index consisting of the BSISO index (summer) and the MJO index (winter). Both the MJO and BSISO modes were constructed by the first and second modes in an extended empirical orthogonal function (EEOF) analysis based on the OLR fields in different seasons. The EEOFs of the MJO mode were constructed from the OLR from December to February, whereas those of the BSISO mode were constructed from the OLR during June to August. In this study, we focused only on the BSISO index. An intra-seasonal time filter was defined as a band-pass filter with cutoff periods of 25 days and 90 days. An EEOF analysis with lags of −10, −5, and 0 days was applied to the intra-seasonal time-filtered CERES OLR data from June to August over the period from 1979 to 2009 in the tropics between 30°S and 30°N. Figure 1 shows the spatial and temporal structure of the obtained EEOF1 and EEOF2 modes. The spatial pattern of the BSISO mode can be captured well by the linear combination of the first and second EEOF modes. By projecting the intraseasonal time-filtered OLR fields with lags of −10, −5, and 0 days onto the first two EEOFs, the principal components (PCs) of the BSISO mode were obtained. Notably, a difficulty regarding the method arose when we used a typical band-pass filter. In particular, the OLR data from after the simulation period that was needed to obtain the time-filtered OLR data were unavailable for the simulations, preventing an analysis of the prediction skill. Instead, the OLR anomaly in this study was obtained by several other processes, following Kikuchi et al. (2012). Specifically, we performed a subtraction of the climatological mean and the three harmonics of the climatological annual cycle of the CERES OLR fields for the period of 2000–2017 and of the 40-day mean of the previous 40 days. Then, we applied a 5-day tapered running mean. Hereafter, we refer to this process as the real-time monitoring (RTM) method. In this study, we first calibrated the OLR field by subtracting a time-evolved OLR bias based on each model, which was averaged independent of the initial date, to distinguish the inherent intra-seasonal variability and the time-evolving OLR bias. The common model bias for each model was calculated by the composite of the differences in OLR between the model's OLR field and the CERES OLR field across all the experiments for NICAM, the ECMWF, the UKMO, and the JMA. This bias correction was originally proposed by Ichikawa and Inatsu (2017), who demonstrated that this method is effective when comparing the prediction skill of the ISO in different models. PCs were then obtained by projecting the OLR anomalies extracted by the RTM method with lags of −10, −5, and 0 days onto the EEOFs. When the OLR fields before the initial day of the simulation were required during the calculation of the PCs, we used the CERES OLR field instead.

Fig. 1.

Spatial-temporal pattern of the first (left) and second (right) modes in the extended empirical orthogonal function (EEOF) analysis based on the OLR fields during June to August. OLR values are multiplied by one standard deviation of the corresponding principal components.

As an example, Fig. 2 shows the time evolutions of the BSISO mode obtained by CERES (black) and by a hindcast simulation by NICAM (red) for the period of 1 to 31 August 2008, using the phase space representation composed of the PCs. BSISO phases of 1–8 correspond to the geographic location of the center of the convective envelope. Figure 3 shows the life cycle of the BSISO mode constructed by EEOF1 and EEOF2, which exhibits good agreement with the discussion in Section 1, with a convective anomaly associated with the BSISO mode appearing in the equatorial central Indian Ocean in phase 1. Then, the northward and eastward propagation of the envelope is pronounced in the Indian Ocean (phases 2–4), which creates the northwest–southeastward-tilted OLR band. The convective anomaly is reinforced around the Philippine Sea in phase 5, and its subsequent northward propagation is dominant in phases 6 and 7. In phase 8, a seesaw of the negative OLR appears between the equatorial Indian Ocean and the Philippine Sea. Therefore, the time evolution of the PCs in the phase space corresponds to the propagation of the projected convective envelope into the EEOF modes. In the following sections, we refer to the rotation speed of the PCs on the phase space as the “phase speed” of the BSISO mode. Therefore, the phase speed inherently reflects the propagation speed of the OLR anomalies both in the eastward and northward directions because the phase speed indicates the speed of the transitions between the defined BSISO phases shown in Fig. 3. Wang et al. (2018) have shown that the time evolution of the PCs given by Kikuchi et al. (2012) clearly tracks the propagation of OLR anomalies both in the eastward and northward directions. Further, phase 0 is defined as the case in which the amplitude of the BSISO is smaller than 1 (r < 1).

Fig. 2.

Phase space representation of the state of the BSISO mode (PC1, PC2) for the period of 1–31 August 2009 by CERES (black) and NICAM (red). The state of the BSISO categorized as eight phases following Kikuchi et al. (2012). The approximate positions of the major convective area are denoted. A star mark shows a location of a pair of (PC1, PC2) at a date, and a circle shows a location of PCs at 1, 16, and 31 August 2009. Please see the explanation described in the main text.

Fig. 3.

A life cycle of the BSISO mode based on averaged OLR anomalies for each phase with an assumption that .

2.3 Verification Scores for the BSISO forecast

Three verification scores for the BSISO forecast were calculated using the bimodal index introduced in Section 2.2, including the root mean square error (RMSE), bivariate correlation (COR), and amplitude error (AERR). The definitions of these were derived by Matsueda and Endo (2011) as follows:   

  
  

Here, PCmodel (t, τ) = (PC1model (t, τ), PC2model (t, τ)) is the PC of the BSISO mode in the model that was initialized at time t − τ and is valid at time t, PCobs (t) = (PC1obs (t), PC2obs (t)) is that in the observation, and θ is the angle between the PC components of the observation and the model on the phase space based on PC1 and PC2 (Fig. 2). The verification scores were calculated using the values of PC1 and PC2 following Eqs. (1)(3). In this study, the prediction skill of the BSISO was measured by the maximum lead time at which COR exceeded 0.6, meaning that the θ on the phase space was smaller than approximately 50° (cos 50° ≈ 0.6; e.g., Miyakawa et al. 2014; Wang et al. 2019).

3. Prediction skill of BSISO

We first addressed the overall prediction skill of the BSISO using NICAM. The RMSEs and CORs were calculated using the NICAM hindcast data with the initial day from 1 to 31 August from 2000 to 2014. Figure 4 shows the averaged RMSEs and CORs as a function of the forecast day. The RMSE increased almost linearly from the beginning of the simulations to approximately the 17th forecast day, but slowed thereafter. The CORs of NICAM reached a score of 0.6 on the 24th forecast day, indicating that the lead time of the prediction skill of BSISO is approximately 24 days in NICAM. Supplemental Figure S1 shows the RMSEs and CORs obtained by NICAM, the JMA, the ECMWF, and the UKMO for 2007–2012. It was difficult to estimate the corresponding lead times of the operational centers from the available TIGGE data because of their limited simulation periods. Although the RMSEs and CORs in NICAM exhibited better scores for 2000 to 2014 than for 2007 to 2012, the NICAM prediction skill was between those of UKMO and ECMWF. NICAM demonstrated similar scores to the operational models for BSISO predictions (Fig. S1). These results are significant when considering that we did not use an initialization technique such as a data assimilation system comprising NICAM and the local ensemble transform Kalman filter, as developed by Terasaki et al. (2015).

Fig. 4.

A time series of (a) Root Mean Square Error (RMSE) and (b) bivariate correlation (COR) on the average of the NICAM hindcast data with the initial day of 1–31 August for the period of 2000–2014.

The question remains as to how the prediction skill in NICAM depends on the initial state of the BSISO (amplitude and phase). Because a better performance for the large initial amplitude of the BSISO is common to other GCMs (e.g., Lee et al. 2015; Wang et al. 2019), a simple ensemble mean of the prediction skill based only on the initial phase may be misleading for the observed BSISO cases where the initial amplitude may not necessarily be uniform among all the phases. To quantify the dependence of the prediction skill on the initial amplitude and initial phase, the prediction skill was estimated by the composite at each 0.5 × 0.5 grid in the PC phase space based on the initial amplitude of the phase in the simulation (Fig. 5a). As reported by previous studies, NICAM has a better prediction skill for BSISO events with a large initial amplitude. In the several cases of the BSISO with a large amplitude (r > 2), the prediction skill could be higher than a lead time of 30 days. Figure 5d shows the prediction skill averaged by the initial phases. The prediction skills from the initial phases 7 to 3 were higher than the average, except phase 2, whereas the minimum of the prediction skill was for phase 4, which had a lead time of approximately 22 days. The model appeared to have a lower prediction skill of the BSISO during phases 4–6, during which the convective center was located over the Maritime Continent.

Fig. 5.

Phase space representation of composites of (a) the prediction skill of the BSISO, (b) the sum of the amplitude prediction error during the experiments, and (c) the number of the cases for the initial condition of the BSISO in August 2000–2014 binned for the initial phase and amplitude of the BSISO. The prediction skill is defined by the period during which COR > 0.6. The dotted bin in (a) indicates a bin which includes the number of the cases more than four. The sum of AERR is calculated only for the experiments during which the amplitude of the BSISO is maintained over unity. (d) The prediction skill of the BSISO (day) and the number of the cases averaged by the initial phases of the BSISO. The red dot line denotes the averaged prediction skill in the all cases.

As to the evolution of the simulated BSISO, Fig. 5b summarizes the tendency of its amplitude. This was measured by the composite of the sum of AERR based on the initial amplitude and initial phases across the period from the initial date to the 25th forecast days, which was approximately the prediction skill score of NICAM (Fig. 4). The composite was performed only via experiments in which the amplitude of the BSISO was greater than 1 throughout the forecast period from 0 to 25 days. The sum of the AERR was negative in the overall phase space, indicating that NICAM tended to decay the BSISO compared with the observations, regardless of the initial amplitude and phases of the BSISO. Kikuchi et al. (2017) reported that the simulated amplitude of the BSISO in the AMIP-type climatological simulation of NICAM was weaker than observation by a factor of ∼ 2. Therefore, the composite in this study might reflect a transition process from the initial condition to the climatological state of NICAM in terms of the BSISO.

We further examined the extent to which the simulated bias is attributable to biases in the phase speed and the amplitude in the phase space (Fig. 6). To confirm how model biases evolve with the integration period, the phase speed and change in amplitude were estimated for each bin (0.5 × 0.5 in the PC phase space) for the entire simulation period (0–30 days) and the latter half of the simulation period (15–30 days). Figure 7 summarizes the average phase speeds for phases 1–8 in both CERES and NICAM. As introduced in Section 2, the phase speed is defined by the angle of rotation in the counterclockwise direction. In the CERES observations (Figs. 6a, b), the phase speed was relatively smaller in phases 3 and 4, in which the convective envelope was located around the Maritime Continent, while large phase speeds were found in phases 6 and 7, in which the convective envelope was over the Philippine Islands. In the NICAM simulation (Figs. 6c, d), the overall distribution of the phase speed was well simulated compared with those for the observations. This can also be confirmed by the ensemble average of the phase speed in each phase (Fig. 7). In addition, NICAM reproduced bins with large phase speeds near r ∼ 1. However, the large phase speeds in phases 6 and 7 in CERES were not well simulated in NICAM, and this bias became large in the composite during the integration period of 15–30 days regardless of the BSISO amplitude (Figs. 6b, d, 7). The change in the amplitude of the BSISO per day on the phase space in the model was also compared with that in the OLR data observed by the CERES (Figs. 8, 9). Although the average decay rate was relatively well simulated in phases 2–4, the simulated BSISO tended to decay faster according to the observations. In particular, the averaged decaying rate of the amplitude in NICAM was large in phases 7–1, during which the convective envelope propagated to the north of the Philippines and a new convective envelope emerged over the western Indian Ocean (Fig. 9). In phase 3, NICAM especially decreased the BSISO amplitude compared with observations in the case of large BSISO amplitude (r > 2; Figs. 9a, c).

Fig. 6.

Phase space representation of the composites of the phases speed for the phase and amplitude of the BSISO using the integration period (a) during 0–30 days and (b) during 15–30 days in CERES, and (c) during 0–30 days and (d) during 15–30 days in NICAM, respectively. The phase speed is defined by the angle of the rotation over the phase space by one day.

Fig. 7.

The averaged phases speed for the phase of the BSISO by NICAM (black) and CERES (red) using the integration period (a) during 0–30 days and (b) during 15–30 days, respectively. The phase speed is defined by the angle of the rotation over the phase space by one day.

Fig. 8.

Phase space representation of the composites of the change of the amplitude per day for the phase and amplitude of the BSISO using the integration period (a) during 0–30 days and (b) during 15–30 days in CERES, and (c) during 0–30 days and (d) during 15–30 days in NICAM, respectively

Fig. 9.

The averaged change of the amplitude for the phase of the BSISO by NICAM (black) and CERES (red) using the integration period (a) during 0–30 days and (b) during 15–30 days, respectively.

To understand the biases in the EOF phase space (6–8 and 1), we examined the composites of the OLR fields. We focused on the OLR fields in phases 6–8 and 1. Therefore, a lag composite analysis was performed on day 0 of phases 5 and 8. As shown in Section 2.2, the EEOF mode was constructed from the filtered observed OLR fields with lags of −10, −5, and 0 days, and the PCs were obtained by projecting the filtered OLR fields with lags of −10, −5, and 0 days onto the first two EEOF modes. The composites of the OLR fields at lags of −10, −5, 0, 5, and 10 days based on phase 5 in CERES and NICAM are shown in Fig. 10. Notably, a lag of 10 days in phase 5 roughly corresponded to a lag of 0 days in phases 6 or 7, in which the phase speed in NICAM was particularly slow compared with the observations (Fig. 6). The composite was performed when day 0 was included during the integration period from 15 days to 30 days. The temporal evolution of the OLR patterns in phase 5 was well simulated in NICAM even two weeks from the initial date of the simulation. This result is significant because it is widely known that the northwest–southeastward tilted rain band associated with the BSISO is difficult to simulate in most GCMs (Sperber and Annamalai 2008; Sabeerali et al. 2013; Sperber et al. 2013; Neena et al. 2017). However, the propagation of the active convective region over the Philippines and the Bay of Bengal appeared to stagnate in NICAM at a lag of 10 days. The composites of the OLR fields based on phase 8 are also shown in Fig. 11. The temporal evolution of the OLR fields during phase 8 was largely influenced by the propagation properties in phases 6 and 7 based on lags of −10 days and −5 days. The active convective region over the Philippines and the Bay of Bengal remained with a lag of 0 days in phase 8, as well as with a lag of 5 days over statistical significance, leading to the collapse of the BSISO structure in NICAM. In addition, the re-initiation of the new convective envelope with lags of 5–10 days in CERES was not clear in NICAM. The incoherent structure of the convective anomaly with lags of 5–10 days in NICAM led to the fast decay of the amplitude of the BSISO in phases 8 and 1 in NICAM (Fig. 8). This is also consistent with the worse prediction skill scores in NICAM from the initial phases of 4–6 because many of the simulations that began from phases of 4–6 passed through phases 6 and 7 within a simulation period of 30 days. According to previous studies based on a linearized atmospheric model (Annamalai and Sperber 2005; Jiang and Li 2005), re-initiation occurs as a response to an anomalous negative heating source over the central and eastern Indian Ocean, corresponding to the positive OLR anomalies in CERES with a lag of −5 days in phase 8 (Fig. 11). The stagnated negative OLR anomalies remain with a lag of −5 days over the Bay of Bengal in NICAM. Therefore, it is expected that the bias of the anomalous convective heating distribution causes the bias in the re-initiation process with a lag of 5 days. In the next section, the possible environmental effect on the slower northward propagation of the convective center over the Philippines and the Bay of Bengal is discussed.

Fig. 10.

Composite of the spatial and temporal time-filtered OLR patterns (W m−2) at day −10, −5, 0, 5 and day 10 on the phase 5 (left) in CERES and (right) in NICAM, respectively. The dotted area indicates that the anomalies are statistically significant at a 95 % confidence level using the Student's t test.

Fig. 11.

The same as Fig. 10, but for on the phase 8.

4. Discussion: Why does the simulated BSISO exhibit slower northward propagation over the Philippines?

The dynamical mechanism for the northward propagation of the BSISO has previously been discussed separately from that of the eastward propagation (e.g., Zhu and Wang 1993; Wang and Xie 1997; Vecchi and Harrison 2002; Jiang et al. 2004; Drbohlav and Wang 2005; Bellon and Sobel 2008; Kang et al. 2010; DeMott et al. 2013). DeMott et al. (2013) showed that the relative importance of the northward propagation mechanism of the BSISO is highly dependent on the location of the convective envelope. In this study, the northward propagation mechanism of the convective envelope over the Philippines was revisited using the prognostic equation of the column-integrated moist static energy (MSE), with the view that the BSISO could be regarded as a “moisture mode” as the MJO (e.g., Neelin and Yu 1994; Sugiyama 2009; Sobel and Maloney 2012, 2013). The column-integrated MSE budget analysis has previously been studied using simulation outputs (Maloney 2009; Wang et al. 2016) and reanalysis datasets (Kim et al. 2014; Yokoi and Sobel 2015). One of the advantages of the MSE budget analysis is that it includes no large cancellation between the vertical advection and other source terms, as in the moisture budget analysis (Yokoi and Sobel 2015).

Following Yokoi and Sobel (2015), the MSE denoted by h is defined as h = s + Lq, where s = CpT + gz, T and z denote the temperature and geopotential height, respectively, q denotes the specific humidity, g denotes the gravitational acceleration, and L denotes the specific heat of vaporization. The prognostic equation of the column-integrated MSE is expressed as   

where is a time-filtered variable using a band-pass filter with cutoff periods of 25 days and 90 days. In addition, u, v, and w denote the horizontal and vertical wind components, R denotes the radiation (calculated as the sum of the forcing caused by the longwave and shortwave radiation), and E and H denote the sensible and latent heat fluxes, respectively. The angle brackets indicate the mass-weighted vertical integration over the troposphere,   
where ps and pt denote the pressure at the surface and the top boundary with pt = 100 hPa and ps = 1000 hPa, respectively.

To investigate the leading mechanism of the northward propagation of the BSISO near the Philippines and the Bay of Bengal, we calculated the regression coefficients of the column-integrated budget terms using Eq. (4) on the column-integrated MSE at a base point of (120°E, 5°N) for 1–31 August from 2000 to 2014 using the ERA-Interim dataset. Figure 12 shows a horizontal map of the regression coefficients of all the terms in Eq. (4) as well as at the other points between at the base point. The regression coefficient of shows the northwest-southeastward tilted structure, which was consistent with the BSISO structure when the convective envelope was located over the Philippines. The positive regression coefficient of radiation has a similar spatial pattern to that of , indicating that the radiation contributed to the maintenance of the MSE anomaly. Over the Indian subcontinent, the regression coefficient of is positive north in relation to that of , suggesting that the phase progression over the Indian subcontinent was mainly caused by the terms of , which is consistent with the MSE analysis performed by Jiang et al. (2018). For the northward propagation over the Philippines, the largest contribution came from the meridional advection term, , which is composed of the meridional advection of the background MSE caused by wind fluctuations of the BSISO and the meridional advection of the MSE anomaly due to background southerlies. This is consistent with Jiang et al.'s (2004) discussion of the northward shift of the low-level moisture associated with the BSISO. The northward propagation over the Bay of Bengal is also largely attributable to , as confirmed by Fig. S2 based on the results of the MSE budget analysis at a base point over the Bay of Bengal (90°E, 10°N). The vertical advection term also contributes to the northward propagation mechanism over the Philippines north of the base point. With regards to the MJO, as the boreal-winter counterpart of the BSISO, the vertical velocity in the eastern part of the convective area exhibits a bottom-heavy profile with shallow convection (Kemball-Cook and Weare 2001; Kikuchi and Takayabu 2004), making positive. The results shown in Fig. 12 imply that a similar contrast might exist between the northern and southern parts of the convective area of the BSISO.

Fig. 12.

Horizonal maps of regression coefficients of the column-integrated MSE budget terms on the column-integrated MSE at a base point of (120°E, 5°N) on August from 2000 to 2014 in the ERA-Interim reanalysis data (W m−2). The terms are column-integrated MSE zonal, meridional and vertical advection ( and ), radiation , and the sum of the sensible and latent heat flux (E + H). The contour shows the correlation coefficient between the column-integrated MSE at a base point of (120°E, 5°N) and at the other points with the contour interval of 0.3. The green circle denotes the base point of (120°E, 5°N).

The improvement of the biases in the tropical mean state leads to improved accuracy in the propagation of the convective cells (Innes and Slingo 2003; DeMott et al. 2015; Weber and Mass 2017). Hence, the reduction of the bias of the background fields is one of the key elements in the realistic representation of the BSISO structure. In this study, we hypothesized that the propagation error of the simulated BSISO in NICAM could be attributed to biases in the background atmospheric conditions related to the northward propagation mechanism. We examined the background column-integrated specific humidity that was averaged for the simulation period in ERA-Interim and NICAM (Fig. 13). In August, the background column-integrated specific humidity increased from the Indian Ocean towards the Eurasian Continent and the Maritime Continent. In NICAM, the column-integrated specific humidity was larger over the western Indian Ocean than in ERA-Interim, whereas this was smaller over the western Pacific and the Bay of Bengal. This means that the zonal gradient of the background specific humidity over the Indian Ocean was smaller in NICAM. The moisture advection was partly realized by the BSISO wind in the presence of the background zonal MSE gradient, and this bias could lead to slower northward propagation over the Indian Ocean, in turn causing fast decay of the BSISO amplitude during phases 7–1. This type of climatological bias in the background moisture fields in NICAM was discussed by Kodama et al. (2015). They suggested that improvements in the model physics, such as the microphysics scheme, the land and ocean model, and the boundary layer scheme, could be key factors in reducing the model biases of the background fields. In addition, improvement of the surface fluxes is also important to achieve a better representation of the ISO (Matsugishi et al. 2020).

Fig. 13.

Horizonal maps of (upper) the column-integrated humidity (kg m−2), the zonal wind velocity (m s−1) at 850 hPa and at 200 hPa (middle), and the meridional wind velocity (m s−1) at 1000 hPa (bottom) averaged for the simulation period in the ERA-Interim and NICAM, respectively. The right column shows the difference of the column-integrated humidity and the zonal wind at 850 hPa, 200 hPa.

The zonal difference in the specific humidity suppresses the convection over the Maritime Continent and enhances the convection over the west of India, which occurs concurrently with the differences in the averaged background zonal wind velocities shown in Fig. 13. As theoretically discussed by Jiang et al. (2004), the bias in the northward propagation may be attributable to the magnitude of the easterly zonal wind shear. However, no large bias existed for the easterly shear simulated in NICAM over the Philippines and the Bay of Bengal, suggesting that the barotropic vorticity mechanism is not the main factor in the slow propagation in NICAM. The background meridional wind component at 1000 hPa is shown in Fig. 13. The southerlies over the Philippine Sea to the Bay of Bengal (90°–130°E), where the easterly shear is weaker, were underestimated in NICAM, at approximately 50 % of those of the ERA-Interim dataset. The northward gradient of the background moisture field was almost negligible around the Philippines, and this means that the propagation speed of the convective envelope caused by the meridional advection was also 50 % of the observation in NICAM, which was approximately consistent with the bias of the phase speed shown in Fig. 5d. The meridional advection of the moisture anomaly by the background southerlies is likely one of the essential factors in the northward propagation mechanism of the BSISO (Figs. 12, S2). Hence, the bias of the background southerly in NICAM is one of the key components for improving the prediction skill of the BSISO.

5. Summary

The prediction skill of the BSISO of a global cloud system resolving atmospheric model (NICAM) was statistically examined using a forecast skill score based on the bimodal ISO index (Kikuchi et al. 2012; Kikuchi 2020). Using the basic verification scores, NICAM exhibited a prediction skill of approximately 24 days for the BSISO. Although it could be misleading to compare our results with those of Miyakawa et al. (2014) because the simulation settings were not the same, the prediction skill of the BSISO was lower than that of the MJO, which is consistent with the conclusions of other previous studies (Lee et al. 2015; Wang et al. 2019). We found that NICAM had a higher prediction skill for the initial phases of 7 to 1, where the convective envelope began to aggregate over the west of the Indian Ocean. In addition, the phase speed and the change in the amplitude of the BSISO were fairly represented in NICAM compared with the observations, although the phase speed was slow in phases 6 and 7 and the change in the amplitude was underestimated in phases 8–1. The fast decay of the PCs in NICAM in phases 8–1, during which the convective center was located over the north of the Philippines, was caused by the stagnant behavior of the convective cells over the Philippines. The propagation and the phase structures while the convective center was over the Maritime Continent and the western Pacific were simulated well, leading to a clear northwest–southeastward-tilted OLR structure in NICAM. Based on an MSE budget analysis, the slower northward propagation over the Philippines could be attributed to the smaller background southerlies in NICAM over the Philippines. The difference in the large-scale circulation was likely associated with the difference in the moisture field and the associated background monsoonal circulation. Therefore, the model physics controlling the background moisture and the concurrent monsoonal circulation are important factors in improving BSISO prediction skill.

In this study, we conducted a series of NICAM simulations for the month of August with a 30-day integration. In future, it would be useful to conduct simulations for other months because the BSISO tends to prevail all the way from June to October (e.g., Kikuchi et al. 2012). This will enable us to assess the prediction skill of the BSISO more rigorously. Further, a set of hindcast simulations with a longer integration period would enable us to conduct an MSE budget analysis using the model's field directly (Fig. 12), which would be useful for discussing the propagation and re-initiation mechanisms in the model. Furthermore, we could then examine the possibility of improving our prediction skill of the BSISO. Given the BSISO timescale, it is widely believed that accurate prediction of the boundary condition is important as is using precise initial conditions for better BSISO prediction. Recent assessments of the prediction skill of the BSISO have been made based on coupled GCM results (Lee et al. 2015; Jie et al. 2017; Wang et al. 2019). Recently, Miyakawa et al. (2017) conducted MJO simulations using a high-resolution atmosphere-ocean coupled model (NICOCO) that coupled NICAM with a global ocean general circulation model, COCO (Hasumi 2006). In a future study, we will be able to examine the effect of the coupling on the forecast skill of the BSISO. Regarding the initial condition issue, many previous studies have reported that the variability of the BSISO is underestimated in almost all existing reanalysis datasets (e.g., Fu and Wang 2004; Jiang et al. 2019). Fu et al. (2009) and Fu et al. (2011) have shown that the prediction skill of the BSISO could be improved by using rectified initial conditions. Finally, better representations of the stratospheric dynamics may also affect the BSISO propagation and initiation through the stratospheric control mechanism of the ISOs (e.g., Son et al. 2017).

Supplements

Supplemental Figure S1 (Fig. S1) shows a time series of the (a) root mean square error (RMSE) and (b) bivariate correlation (COR) on the average of the NICAM hindcast data (black curves) and the forecast data from the ECMWF (red), UKMO (green) and JMA (blue) with the initial day of 1–31 August for the period of 2007–2012.

Supplemental Figure S2 (Fig. S2) shows the same information as Fig. 12, but at a base point of 90°E, 10°N.

Acknowledgments

We would like to thank Dr. Tatsuya Seiki, Dr. Yohei Yamada, Dr. Tomoki Ohno and Dr. Satoru Yokoi for their many useful comments and discussions. We also greatly appreciate the two anonymous reviewers for their constructive comments. This work is supported by the FLAGSHIP2020 within the priority issue 4 (Advancement of meteorological and global environmental predictions utilizing observational “Big Data”) and Program for Promoting Researches on the Supercomputer Fugaku” (Large Ensemble Atmospheric and Environmental Prediction for Disaster Prevention and Mitigation) by MEXT. KK acknowledges the support of NOAA grant NAl70AR4310250 and the support by JAMSTEC through its sponsorship of research activities at the IPRC (JICS). This study was also supported in part by a Grant-in-Aid for Young Scientists from JSPS (19K14801). All numerical experiments were run on the K computer (proposals hp120313, hp130010, and hp140219), and all figures shown in the paper were created using the Dennou Club Library (DCL).

References
 

© The Author(s) 2021. This is an open access article published by the Meteorological Society of Japan under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.
https://creativecommons.org/licenses/by/4.0/
feedback
Top