Improved intraseasonal variability in the initialization of SINTEX‐F2 using a spectral cumulus parameterization

A newly developed spectral cumulus parameterization (spectral scheme) was implemented in the Scale Interaction Experiment‐Frontier version 2 (SINTEX‐F2) seasonal prediction system to improve intraseasonal variability in the system initialization. A simple sea surface temperature (SST) nudging scheme using different SST data and restoring times was used to initialize the system, and the initialized atmosphere obtained from both the original convection scheme (Tiedtke scheme) and the new spectral scheme was evaluated against observational data. It was found that that climatology and variability simulated by the spectral scheme were comparable to those simulated by the original scheme. In addition, the intraseasonal variability represented by the Madden–Julian oscillation (MJO) was better simulated by the spectral scheme than the original scheme. An analysis of the structure of the organized convection revealed the successful simulation of low‐level shallow convection before the peak of the organized convection by the spectral scheme when compared with the observation, a result lacking in the original scheme simulation. In addition to the positive qualitative results, a statistical and quantitative analysis showed that the spectral scheme captured the MJO‐related variability better than the original scheme. In conclusion, the prediction system using the spectral scheme is expected to improve seasonal predictions for seasonal variability whose evolution is affected by intraseasonal variations.


| INTRODUCTION
Although the ability of recent seasonal prediction systems to capture climate variability is improving (e.g., Weller and Cai, 2013;Bellenger et al., 2014), greater prediction skill is still required even for the most recent seasonal prediction systems. The reason for this is that predictions using numerical models remain imperfect, and more accurate prediction is needed for stakeholders to mitigate the risks of climate variability, which influences natural disasters, various resources and human health (Sahu et al., 2012;Iizumi et al., 2013;Takaya et al., 2014). Therefore, it is necessary to understand the limits of the predictability of the system and to improve the prediction skill by updating numerical schemes employed in the system. To improve seasonal predictions, past studies have developed improved general circulation models (GCMs) of the system (including updates of physical parameterizations) and initialization techniques (e.g., Hudson et al., 2011;Merryfield et al., 2013;Saha et al., 2014;Johnson et al., 2019).
The Scale Interaction Experiment-Frontier (SINTEX-F) prediction system was collaboratively developed by the European Union and Japan (Luo et al., 2005a(Luo et al., , 2005b for seasonal predictions based on an earlier climate model (SINTEX, Gualdi et al., 2003;Guilyardi et al., 2003). This prediction system was found to provide good predictions of interannual variability originating from the coupledmode phenomena in the tropics, such as the El Niño Southern Oscillation (ENSO, Luo et al., 2008b), and the Indian Ocean Dipole (Saji et al., 1999;Luo et al., 2008a). However, there remains uncertainty in the system, such as poor predictions of extratropical variability, and furthermore, the model sometimes failed to predict distinct events in the tropics (Doi et al., 2016). SINTEX-F2 was developed to further improve the prediction skill of SINTEX-F, employing higher resolution and improved physical processes for both the atmospheric and oceanic components of the coupled GCM (CGCM; Masson et al., 2012). While the new system improved seasonal predictions compared with the preceding system, especially for mid-latitudes, the overall prediction in the tropics did not significantly improve (Doi et al., 2016) without an additional assimilation technique for the ocean component (Doi et al., 2017).
Improved representation of atmospheric, oceanic, and coupled-mode variability in the simulation is considered to be critical for enhanced prediction skill. Uncertainty in the atmospheric model is generally more significant than the ocean model, due to the large uncertainty existing in the parameterization of clouds (Baba, 2019). One of the significant cloud-related tropical variability is the Madden-Julian oscillation (MJO, Julian, 1971, 1972) which has an intraseasonal timescale. Several past studies suggested that high-fidelity in simulating the MJO is important for predicting seasonal variability, that is, the MJO is influential in the initiation, development, and termination of the interannual variability (e.g., Takayabu et al., 1999;Rao and Yamagata, 2004;Hendon et al., 2007;Rao et al., 2009;Zhang, 2013;Baba, 2020b). It is therefore considered that a more accurate simulation of MJO features may contribute to improved seasonal predictions. However, simulating reasonable features of the MJO with realistic climatology remains a challenging task even for the current climate models (Hung et al., 2013;Baba, 2019;Baba and Giorgetta, 2020). The reason for this is that the MJO involves various cloud activities consisting of different cloud types during the development that are challenging to resolve in climate models (e.g., multi-cloud structure named as "self-similarity," Kiladis et al., 2009;Riley and Mapes, 2011).
To simulate cloud activities, a well-tuned or improved convection scheme is required. In general, improvements in the simulation of the MJO using only tuning of the parameters are difficult to achieve (Hung et al., 2013). Moreover, the tuning methods and the tuned parameters are generally considered as model dependent; for example, an identical convection scheme indicates different performances for MJO in different GCMs (Crueger et al., 2013(Crueger et al., , 2018. Past studies analysed roles of each physical process or developed new parameterizations to realize better MJO simulation (e.g., Fu and Wang, 2009;Zhang and Song, 2009;Crueger et al., 2013;DeMott et al., 2015;Liu et al., 2018;Yang and Wang, 2019). Improving representation of shallow convection was found to contribute to better MJO simulation (Zhang and Song, 2009;Liu et al., 2018;Yang and Wang, 2019), but general validity of each developed parameterization has not been verified in other GCMs. Although atmosphereocean coupling was also found to improve MJO simulation, the improvement depends on parameterizations of atmospheric model (Crueger et al., 2013;DeMott et al., 2015).
For an improved representation of convective clouds regardless of the GCMs, Baba (2019) recently developed a new spectral cumulus parameterization that can model unresolved clouds better than other some existing schemes. Previous studies have revealed that this convection scheme can realize similar improvement in different GCMs in terms of climatology and climate variability, especially for the MJO (Baba, 2019(Baba, , 2020cBaba and Giorgetta, 2020). SINTEX-F2 requires initialization before the prediction, so improving the fidelity of the intraseasonal variability by using this scheme in the initialization (and across the prediction) may further improve the prediction skill of SINTEX-F2, particularly in the tropics. These improvements involve adaptations to the convection scheme in the initialization, so it is necessary to examine the validity of the replaced convection scheme before conducting prediction, because a change in the convection scheme has a significant influence on the climatology of the model (Kim et al., 2011;Hung et al., 2013). The purpose of the present study is to evaluate the simulated intraseasonal variability in the initialization of SINTEX-F2 before the prediction when using the newly developed convection scheme of Baba (2019). The new convection scheme is implemented in SINTEX-F2, and the effects of climate variability during the initialization are evaluated for the purpose of ensuring that the basic performance of the model is retained. The statistical features of the MJO-related variability in the initialization which can influence the evolution of the interannual variability (Baba, 2020b) are examined in detail.
This article is structured as follows. The details of the model and experimental setup, including the details of SINTEX-F2 and the convection schemes used in this study, are described in Section 2. The results and discussion are presented in Section 3, which first evaluates the climatology and tropical variability simulated by the new scheme and then analyses the intraseasonal variability. For the intraseasonal variability, a qualitative analysis using MJO diagnostics, composite, and statistical analyses are presented. Section 4 provides a summary and conclusions.
2 | MODEL AND EXPERIMENTAL SETUP

| SINTEX-F2
SINTEX-F2 is a seasonal prediction system consisting of a CGCM. The atmospheric GCM (AGCM) coupled in SINTEX-F2 is ECHAM5 (Roeckner et al., 2006) which has a horizontal resolution of 1.125 (T106) and 31 vertical layers. Compared to SINTEX-F, the vertical resolution of the atmospheric component is increased from 19 to 31 levels, and the physical parameterizations used in the model are updated from ECHAM4 (Roeckner et al., 1996). The ocean GCM (OGCM) coupled with the AGCM is the Nucleus for European Modelling of the Ocean (NEMO, Madec, 2008), which employs a higher horizontal resolution than the previous version of the OGCM in SINTEX-F; the present OGCM employs the ORCA05 grid (0.5 resolution) rather than the ORCA2 grid (2 resolution). Furthermore, the updated model includes a sea-ice model. The AGCM and OGCM are coupled by the Ocean Atmosphere Sea Ice Soil version 3 coupler (Valcke, 2013) in SINTEX-F2. The atmospheric and oceanic fluxes are exchanged by the coupler every 2 hr and no flux correction is applied.
For the initialization, SINTEX-F2 employs a simple sea surface temperature (SST) nudging scheme that uses different SST data and restoring times. The two types of optimum interpolation SST (OISST) data were used for the SST nudging, that is, daily (Reynolds et al., 2007) and weekly (Reynolds et al., 2002) datasets. The restoring time is defined by negative feedback values to the surface fluxes (−2,400; −1,200; and −800 W m −2 ), and each value corresponds to 1, 2, and 3-day restoring times, respectively. Using the two types of SST data and the three feedback values, six ensemble members were generated. The application of different vertical mixing coefficients for the ocean component, which was previously considered to generate ensemble spread was ignored here, because the change in the coefficient was found to generate an insignificant ensemble spread (Doi et al., 2016(Doi et al., , 2017. Using two different convection schemes, the six ensemble members were run for 15 years starting on January 1, 1982, and the model results for the last 14 years were analysed. The prediction run can be conducted by starting a nudging-free run from the first day of each month (e.g., run from January 1 to June 30 for a 5-month prediction) using the initialized condition; however, the prediction run was not performed here, since this study only focuses on the initialization. The convection schemes employed here are the Tiedtke scheme (Tiedtke, 1989) and the spectral cumulus parameterization of Baba (2019) (denoted as spectral scheme hereafter). Details of these schemes are described in the following Section 2.2. Table 1 summarizes the experimental cases, which are roughly classified by the different convection Schemes (TK or SP). For example, within each case, the ensemble members are referred to as TK1dd (TK scheme, 1-day nudging using daily SST data) or TK3dw (TK scheme, 3-day nudging using weekly SST data). Except for the difference in the convection scheme, all other parameterizations and their parameters were identical across the cases.

| Tiedtke scheme
The original convection scheme for ECHAM5 is the Tiedtke scheme (Tiedtke, 1989). The Nordeng's modification which employs convective available potential energy (CAPE) closure for deep convection (Nordeng, 1994) is adopted for the original scheme. For simplicity, the modified scheme is referred to here as Tiedtke scheme. This scheme is a cloud mass flux-based convection scheme which is characteristic in terms that it assumes predefined bulk cloud types (deep, shallow, and midlevel) and switches parameterizations in each grid column depending on the identification criteria (Möbis and Stevens, 2012). Therefore, while the scheme can consider several cloud types, the basic structure is dependent on the presumed cloud types identified by the employed criteria, and the coexistence of different cloud types in a grid column cannot be expressed. Past study have shown that the scheme implemented in ECHAM can simulate the MJO to some extent (Crueger et al., 2013); however, the performance is known to be model dependent (Crueger et al., 2018). The author's previous study revealed that because the parameterizations of each cloud type in this scheme vary greatly due to the tuning of the entrainment, continuous and seamless cloud developments were poorly simulated, leading to a failure in simulating the organized convection of the MJO (Baba and Giorgetta, 2020). Tiedtke scheme also has been used with different formulations in Integrated Forecast System (IFS) of European Centre for Medium-range Weather Forecasts (Bechtold et al., 2008(Bechtold et al., , 2014. The modified convection scheme employed modified entrainment model, criteria for cloud types, and convective closure, and Hiron et al. (2013) reported that relative humidity dependent formulation could improve MJO simulation. Although IFS-type Tiedtke scheme is newer one, the present purpose is to evaluate the spectral scheme comparing with the original scheme of SINTEX-F2, so the newer scheme and related modifications are not considered in this study.

| Spectral scheme
To further improve the convective cloud properties simulated by the model, the spectral scheme of Baba (2019) was implemented in SINTEX-F2. This spectral scheme was developed to simulate convective clouds behaviour that is not well resolved in GCMs by the existing schemes. The main novel parameterization of this scheme is the entrainment parameterization which can be applied to various cloud types and which is given by and 1 2 where ε u is the updraft entrainment rate, B u is the incloud buoyancy, δ u is the updraft detrainment rate, w u is the cloud mean updraft velocity, and z is the vertical coordinate. C 1 = C 2 ≈ 0:2 , a=C 1 + 0:5 , and b =C 2 +0:75 are model parameters determined by a statistical analysis on the convective cloud properties from a cloud-resolving model simulation. In this convection scheme, the updraft convection is expressed using the spectral approach (Arakawa and Schubert, 1974;Chikira and Sugiyama, 2010). Thus, 14 different cloud types are defined in each grid column, characterized by the different cloud base updraft velocities. Equations (1) and (2) are solved with the updraft budget equations for each cloud type in which the vertical profiles of the cloud mass fluxes for each cloud type are computed. The bulk updraft properties are given by a summation of the properties of all cloud types. In contrast to the updraft, the spectral approach is not employed for the downdraft calculation, and instead, the downdraft is calculated by a bulk cloud mass flux. Other details such as the budget equations, cloud microphysics, and the numerical procedures for the convection scheme are described in the appendix of Baba (2019). In the present implementation, a shallow convective closure, which was found to be valid for improving the scheme's capability in simulating intraseasonal variability, was adopted as well as a water budget correction (Baba, 2020c). The author's previous studies found that this convection scheme could improve the tropical variability, such as atmospheric response to ENSO-and MJO-related variability, despite the convection scheme being employed in different GCMs (Baba, 2019(Baba, , 2020cBaba and Giorgetta, 2020). In addition, the scheme has recently been found to be able to reproduce the diurnal cycle of convection over the Maritime Continent better than other existing parameterizations (Baba, 2020a), and to have improved fidelity for global tropical cyclone activity in an AGCM (Baba, 2021).

| RESULTS AND DISCUSSION
3.1 | General performance

| Climatological error
Since the convection scheme is significantly responsible for the biases of GCMs (Hung et al., 2013), the climatological errors simulated by the model using the spectral T A B L E 1 Configurations for the experimental cases. The results from the ensemble mean of each case are denoted as TK and SP, respectively scheme should be evaluated before analysing the performance of the model in simulating the intraseasonal variability. The climatological error was evaluated using standardized error analysis (Reicher and Kim, 2008) as performed in our previous studies (Baba and Giorgetta, 2020;Baba, 2020c). To evaluate the climatology, various observational and reanalysis data were used in the present study, as shown in Table 2. These data, including other observational and reanalysis data, are hereafter referred to as observation for simplicity. Computed standardized errors for the model results are compared in Figure 1. Slightly larger errors are found for SP than for TK for many variables, and the errors in the outgoing longwave radiation (OLR), temperature and zonal wind are larger than for TK ( Figure 1a). However, other variables such as outgoing shortwave radiation (OSR), precipitable water, and meridional wind for SP have smaller errors than those for TK. In particular, the error in the meridional wind for SP is significantly improved over the error for TK. The total relative error of SP compared to TK is 1.17 (i.e., climatological error of SP is 17% larger than that of TK), and thus the errors of both cases can be regarded as comparable. The greater climatological error for SP may be due to the unchanged parameters of physical parameterizations other than convection, which were instead tuned in the original convection scheme. Analysing the errors for different SST nudging data, it was found that the temperature error is significantly dependent on the differences in the SST data, that is, the use of daily SST data could reduce the temperature error in SP ( Figure 1b). This suggests that the climatological temperature error for SP is more sensitive to the SST data than for TK. It should be noted that the differing restoring time did not cause a clear difference in the climatological errors (not shown here).

| Mean state
The simulated annual mean precipitation of the two cases and the observation are compared in Figure 2. TK overestimated precipitation in most of the Pacific and Indian tropical convergence zones. The simulated precipitation is much improved for SP, with a similar amount of precipitation found as in the observation. TK also simulated an unrealistic double intertropical convergence zone over the Indian ocean; this precipitation distribution is more realistic for SP. However, SP underestimated the precipitation in the convergence zone of the Indian ocean as compared with the observation. The inferior feature of precipitation in this area may be caused by the overestimated precipitation in the western Pacific which is commonly seen in both TK and SP. While SP reduced the overestimated precipitation in the tropics, the overestimation in the western Pacific could not be improved, so that the precipitable water over the Indian ocean decreases. Further improvements are necessary for SP to better simulate precipitation distribution over the Indian ocean, and this may be realized by tuning parameters other than the convection scheme.
The energy budget from the radiative fluxes at the top-of-atmosphere (TOA) must be evaluated when the convection scheme is replaced. This is because the TOA energy budget depends on convection and it determines the total energy budget for GCMs, leading to overall warming or cooling trends for the globally averaged temperature. Annual mean distributions of OLR and net shortwave radiation (SW) are compared in Figure 3. Both cases are found to simulate qualitatively similar annual mean OLR and net SW compared with the observation (not shown here), but their biases are different. TK simulated smaller OLR in the subtropics, while SP simulated larger OLR in the subtropics and tropics, than the observation. TK simulated smaller OSR over land and larger OSR over the Southern ocean, while SP simulated larger OSR in the tropics. This also suggests that SP (TK) simulated lower (higher) cloud cover than the observation in tropical ocean, and thus further tuning of the cloud cover may be necessary especially for SP to reduce the radiative bias. The resulting annual and global mean TOA energy budgets are compared in Figure 4. Both cases showed sufficiently small TOA energy budgets  compared with the observation (Loeb et al., 2009). Originally, TK has a positive imbalance in the TOA budget, and this bias does not change significantly between the ensemble members. On the other hand, SP showed a negative imbalance and smaller magnitudes than TK for all ensemble members. Based on these results, although SP simulated a larger radiative bias, it can be concluded that SP simulated comparable climatological radiative properties at the TOA.

| Interannual variability
The ENSO has an interannual timescale and has the most significant impact on seasonal variability, so it is necessary to validate the model's capability to simulate atmospheric response to the ENSO before discussing the variability of shorter timescales. To statistically evaluate the atmospheric response, the regressed monthly precipitation anomaly onto   because this case originally overestimated tropical precipitation as shown in Figure 2. Wind stress acting onto the ocean surface is a major driver of sea surface current in tropics, and its response to the ENSO must be correctly simulated even if the convection scheme is replaced. The response of the simulated wind stress to the ENSO is compared in Figure 6, with regression adopted as for the precipitation. TK overestimated the westerly wind stress in the western north Pacific, while SP underestimated the westerly wind stress in the central equatorial Pacific. However, the results indicate that both cases could successfully simulate qualitative global wind stress response to ENSO.

| MJO diagnostics
To qualitatively evaluate the fidelity of the simulated MJO, MJO diagnostics (Waliser et al., 2009) were applied to the last 14 year of the model results and the observational data. Following the methods of our past studies (Baba, 2019(Baba, , 2020c, for the MJO diagnostics, GPCP (Huffman et al., 2001), Advanced Very High Resolution Radiometer (AVHRR, Liebmann and Smith, 1996), and National Center for Environmental Prediction/National Center for Atmospheric Research reanalysis (Kalnay et al., 1996) were used for the daily precipitation, OLR, and wind velocity, respectively. Convectively coupled equatorial waves (CCEWs) have intraseasonal or more shorter timescales, and they are known to be the dominant tropical atmospheric variability at these timescales. The features of the simulated CCEWs were analysed using a wavenumber-frequency diagram (Wheeler and Kiladis, 1999) for daily precipitation. It was found that TK has a weaker MJO signal than the observation and SP, as seen in the weak power of wavenumber 1 to 3 and a frequency slower than 30 days in the symmetric components (Figure 7). The corresponding wavenumber and frequency ranges for the SP indicated a sufficiently strong signal. Moreover, TK failed to simulate the coupled equatorial Kelvin wave frequency seen in the observation and SP, as demonstrated by the peak power being located at slower frequencies. In contrast, TK has improved spectra for the anti-symmetric component as it simulated the features of the mixed-Rossby gravity wave better than SP, as indicated by its peak power being closer to the observation.
The lag-cross correlation of the daily precipitation and zonal wind anomalies are compared in Figure 8. From the observation, an organized convection occurs in the Indian ocean and moves eastward from the Indian to Pacific oceans, as the precipitation and zonal wind anomalies appear. SP simulated these features well, although the eastward propagation is slightly shorter than the observation. While TK simulated the presence of the organized convection in the Indian ocean, which is indicated by the precipitation anomaly, its eastward propagation is much weaker and shorter than the observation and SP. TK also simulated a weaker zonal wind anomaly than the observation and SP.
The principal components (PCs) corresponding to each MJO phase were calculated using a combined empirical orthogonal function (EOF) analysis for the OLR and wind anomalies. Then, using PCs for estimating the MJO index (defined as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PC 2 1 +PC 2 2 q where PC 1 and PC 2 are the first two leading PCs) and phase (Wheeler and Hendon, 2004), the MJO lifecycle composite was computed. The results are compared in Figure 9. The observation shows that a negative OLR anomaly of MJO appears in the Indian ocean with a convergent wind, and that it propagates towards the western Pacific. TK simulated the OLR anomaly; however, its convective strength, represented by the OLR anomaly, is much weaker than in the observation, and the phase does not correctly correspond to the observational phase, as the location of the peak OLR anomaly is displaced from the correct location. SP simulated a strong OLR anomaly, but it is stronger than the observation. Additionally, in contrast to TK, SP correctly simulated the phase of the MJO compared with the observation, as the locations of the peak OLR anomaly at each phase correspond well to those of the observation.

| Composite analysis
Differing features of the MJO were found for TK and SP, with SP qualitatively outperforming TK based on the MJO diagnostics. However, the reasons for this outperformance remain unknown. A composite analysis of the organized convection during MJO development is useful for clarifying the difference in the structure of the MJO resulting from the convection schemes (e.g., Chikira, 2014; Baba and Giorgetta, 2020; Baba, 2020c). To detect MJO peak time and location during the simulation, a MJO detection algorithm based on the minima of the OLR anomaly as used in Chikira and Sugiyama (2013)  variations. The anomaly was then band-pass filtered for zonal wavenumbers 1-5 and Days 20-100, and the average was taken for latitude between 10 S and 10 N. The following conditions were used to select the centre location of the MJO using the processed OLR anomaly. The minimum of the OLR anomaly must be located between 60 E and 150 W, and the minimum must propagate at least 60 . During the propagation, the minimum must be less than −0.7σ where σ is the standard deviation (SD) of the filtered anomaly. For the detection of MJO in the observation, AVHRR was used in the following analyses. In past studies, the behaviour of the MJO has been analysed using moist static energy (MSE), since the behaviour is known to be dominated by the MSE budget (e.g., Maloney, 2009;Kiranmayi and Maloney, 2011;Sobel et al., 2014;Wolding and Maloney, 2015). Sobel et al. (2014) analysed the MSE budget during the development of the MJO as where hi represents vertical integration (between 1,000 and 100 hPa), c p is the specific heat at constant pressure of air, T is the temperature, L v is the latent heat of evaporation, q is the specific humidity, t is the time, v is the horizontal velocity, ω is the vertical velocity, p is the pressure, E is the latent heat flux, H is the sensible heat flux, and R is the radiative flux. Here, the MSE is defined by where g is the gravitational acceleration, z is the height from a reference level, and s is the dry static energy. Equation (3) neglects gz for the time derivative and horizontal advection; however, only a small difference is known to result from this omission (Sobel et al., 2014). The source term of the radiative flux R is calculated from the difference between the net surface and net TOA radiative fluxes. The advection tendencies in Equation (3) for the observation were calculated using daily ERA5 data as done in Sobel et al. (2014). The validity of the use of reanalysis data to compute advection tendencies was confirmed by the preceding study, and the present results are consistent with the results of that study as shown later. Based on the MJO detection algorithm, the composited MSE anomalies were computed for the observation and each case ( Figure 10). Here, the anomaly is defined as the difference of the value from its average during all relative days. The composited MSE anomaly shows that a positive MSE anomaly occurs from Day −20 to Day +10 (Day 0 corresponds to the time when maximum MJO anomaly appears), and that the anomaly profile trend is tilted with altitude. These features are consistent with the results of preceding studies (e.g., Maloney, 2009;Cai et al., 2013). It should also be noted that the MJO is a moisture-mode variability, and thus the MSE anomaly corresponds well to the moisture anomaly during the MJO development (Chikira, 2014;Baba and Giorgetta, 2020). Only the results for TK lack a low-level MSE anomaly increase around Day −20 to Day −10, and shows an upright MSE anomaly profile against the relative days. The tendencies of the column-integrated MSE were composited and their variations during relative days are compared in Figure 11 to clarify the contributions of each term in Equation (3). First, it was found that the MSEs of the observation and SP have quantitatively similar variations, while the MSE for TK is lower value throughout the relative days. The source terms contributing to MJO maintenance are known to be surface heat and radiative fluxes (Sobel et al., 2014;Wolding and Maloney, 2015). Of these fluxes, differences are observed in the surface heat fluxes during the relative days, indicating that TK shows larger surface heat flux than the observation, while that of SP is smaller. The larger surface heat fluxes of TK are assumed relate to larger mean precipitation shown in Figure 2. Although the heat fluxes differ, their anomalies are smaller than the anomalies in radiative flux, which is a main amplifier of the organized convection by longwave radiative heating (del Genio and Chen, 2015;Wolding and Maloney, 2015;Baba and Giorgetta, 2020). The remaining terms are advection tendency terms that contribute to eastward propagation which trends presented here are qualitatively similar to those of preceding study (Sobel et al., 2014). In particular, horizontal rather than vertical advection plays a dominant role in the eastward propagation. It is found that TK and SP resulted in similar horizontal advection trends, so the contribution of the terms are considered equivalent. On the other hand, vertical advection tendencies between TK and SP Relative day F I G U R E 1 1 Comparison of the composited column-integrated moist static energy (MSE) and the tendencies of each term in Equation (3) during the relative days. The horizontal and vertical advection tendencies are column-integrated. ERA5 and AVHRR are used as the reference vary widely. Large anomalies in the advection around Day 0 were found for the observation and SP, while the anomaly for TK was smaller. Focusing on the beginning of the emergence of difference in the vertical advection tendency, it started to appear between Days −20 and −10. Past studies have noted that the vertical MSE advection is associated with vertical moisture advection (Chikira, 2014;Wolding and Maloney, 2015), so it is assumed that vertical transport by shallow convection was lacked in TK between Days −20 and −10, then its accumulated influence caused the large difference in the tendency around Day 0, that is, MJO anomaly.
The transitions of the vertical profiles of the composited advection tendencies are compared in Figure 12. It was found that horizontal (vertical) advection tendencies act to decrease (increase) the MSE at the altitude lower than mid-level. In particular, the vertical advection increases the low-level MSE throughout the relative days. TK simulated a lower MSE increase in the low-level, while SP simulated a larger MSE increase which is similar to the observation. This difference is considered to relate to the MSE anomaly difference in the low-level as shown in Figure 10. Furthermore, TK simulated a weak MSE decrease above the mid-levels, while SP simulated a strong MSE decrease which may contribute to the tilted MSE anomaly trends. The vertical advection tendency in the MSE was further analysed by decomposing the term into terms of dry static energy (s =c p T +gz ) and specific humidity (L v q). It is clearly shown that the MSE decrease is caused by vertical advection of dry static energy, while the MSE increase is caused by that of humidity ( Figure 13). In general, MJO anomaly appears as a large-scale environment change due to the motion of organized convection, and the tendency does not represent amount of condensation, but represents influence of motion of organized convection. The tendencies of dry static energy and specific humidity then represent adiabatic cooling in mid-level and moistening in low-level by upward transport induced by organized convection, respectively (Chikira, 2014;Sobel et al., 2014). Considering these facts, humidity increase is responsible for the low-level MSE increase. According to the preceding studies, the low-level humidity increase occurs in the advance of the MJO convection (Kiranmayi and Maloney, 2011), stemming from the occurrence of shallow convection (Baba and Giorgetta, 2020;Baba, 2020c). Thus, the presence of shallow convection before the organized convection is considered to cause the difference in the MSE trends between TK and SP. -2,000 -3,000 -4,000 -5,000 0 F I G U R E 1 3 As Figure

| MJO score
Accurate prediction of the MJO requires more advanced ensemble and assimilation techniques than the present simple SST nudging scheme (Hudson et al., 2011(Hudson et al., , 2013Vitart, 2014;Kim et al., 2018), and thus accurate prediction of the MJO is challenging for the present model. However, the statistical features of the simulated MJO are important, since it may contribute to the evolution of the interannual variability as described above. To quantitatively validate the statistical fidelity of long-range simulated MJO, the MJO score proposed by Crueger et al. (2013) was calculated as in the previous studies (Baba and Giorgetta, 2020;Baba, 2021). The MJO score MJO sc is given by where R mean is the mean ratio of eastward to westward propagation power R within the MJO wavenumberfrequency ranges (20-100 day period and wavenumbers 1-3), and F OLR is the total explained variances of the first two leading modes of the combined EOF for the OLR. Note that the propagation power R is calculated for the OLR, precipitation, zonal wind at 850 hPa (U850), and zonal wind at 200 hPa (U200). R mean and F OLR represent the eastward propagating feature and strength of the MJO, respectively. F OLR , R mean , and MJO sc of the model results are compared with the observation in Figure 14.
The results show that the strengths of the MJO as represented by F OLR are similar (as is the ensemble spread of F OLR ), while their eastward propagating features differ. TK indicates smaller R than SP, except for R U200 , giving a poorer MJO score for TK; however, the score is better than previous ECHAMs (Crueger et al., 2013). This indicates that the spectral scheme can further increase the fidelity of the simulated MJO over that of Tiedtke scheme from a statistical perspective, especially for the eastward propagation.

| Strength, time, and location
The strength of the MJO alternatively can be measured by the value of MJO index (Wheeler and Hendon, 2004) which is given by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PC 2 1 +PC 2 2 q as described in Section 3.2.1. To compare the simulated MJO strength with the observation, the time variations are plotted in Figure 15. The observed MJO strength shows seasonal peaks in boreal and austral winter, and the peak of boreal winter is stronger than that of austral winter, as shown by the preceding studies (Zhang, 2013). While both schemes partly captured the peaks of the observed MJO strength, the correlations with the observation are not particularly high. For example, the correlations for the ensemble mean and running mean MJO index between the observation, TK, and SP are estimated as 0.30 and 0.18, respectively. This indicates that the corresponding MJO strength as estimated by the MJO index is challenging to model when using only the SST nudging scheme. Originally, the timescale of atmospheric variability is shorter than the oceanic variability, so it is considered that an additional assimilation technique for the atmospheric component is necessary for the present model to obtain better agreement with the observation for the MJO index, which is estimated from the PCs.
Hovmöller plots of MJO centre locations are shown in Figure 16 to compare the simulated MJO occurrence with the observation. Here, the centre locations of the MJO were estimated from the OLR anomalies as described earlier. Since each Hovmöller plot from the simulation is not normalized by the number of ensemble members, the comparison indicates whether the possible MJO-related variability is captured by the ensemble members of each case (the correspondence considering the number of ensemble members will be shown in Figure 17). The comparison indicates that TK sometimes failed to simulate the OLR anomalies of the MJO, that is, MJO-related variability of some events was missed by TK. Although SP sometimes simulated MJOs even when not present in the observation, this scheme successfully captured the majority of MJO events in the observation via its ensemble members. These results imply that the spectral scheme more correctly considered atmospheric variability originating from the MJO in the initialization, even with the present simple SST nudging scheme. Finally, to quantitatively compare the strength and location of the simulated MJO with the observation, density plots for the OLR intensity and longitude are presented in Figure 17. The observation suggests some peaks of the OLR intensity at different longitudes. While TK partly captured some of the peaks at the different longitudes, the density is insufficient, which indicates that a strong MJO signal was not simulated in the correct locations. In contrast, SP simulated most of the strong OLR intensities in the corresponding locations. Moreover, the density of SP is closer to the observation than TK. Therefore, SP can simulate statistically accurate MJO variability in terms of both strength and location, and the performance is better than that of TK.

| SUMMARY AND CONCLUSIONS
A newly developed spectral cumulus parameterization (spectral scheme, Baba, 2019) was implemented in the SINTEX-F2 seasonal prediction system, and the validity for the initialization was examined. Seasonal prediction skill is considered to be influenced by the representation of intraseasonal variability, and thus the fidelity of the MJO, especially its statistical features which may contribute to improved seasonal prediction skill, were the focuses of the study. Experiments were conducted using two different convection schemes, and six ensemble members generated by the SST nudging scheme, with the two cases denoted as TK (original scheme, Tiedtke, 1989) and SP (spectral scheme), respectively.
Simulated climatology was evaluated first since the replacement of the convection scheme is considered to significantly change the model's climatology. From the results of the standardized error analysis, it was found that the relative climatological error of SP increased compared with TK, but the total errors of TK and SP were comparable. A comparison of the annual mean precipitation revealed that SP could reduce the overestimation of precipitation in TK. However, several known biases remained and the tropical convergence zone in the Indian ocean was weakened in SP. The annual mean TOA radiative properties indicated that while TK could simulate distributions better than SP, the results were qualitatively similar, and the total global energy imbalance for the both cases was found to be sufficiently small. F I G U R E 1 7 Comparison of density plots of longitude versus normalized outgoing longwave radiation (OLR) intensity (anomaly) for detected Madden-Julian oscillation (MJO) events. The value of density is given by the count of MJO passing though the sampling area defined by a 2 interval for longitude, and a 1 interval for the normalized OLR anomaly. The OLR anomaly was normalized by its SD σ. The ensemble mean was taken for the model results. The density plots from the model are indicated by shading, and the observation are indicated by solid lines The simulated atmospheric variability was examined for the ENSO which is known to cause the largest interannual variability in the tropics. The regressed monthly precipitation anomaly onto Niño 3.4 revealed that the overestimated convection response of TK was reduced in SP. Even though the convection scheme was replaced, the qualitative response of the wind stress to the ENSO remained unchanged, as found by the regressed wind stress anomaly. Therefore, the performance of the model in terms of interannual variability is expected to be comparable regardless of the convection scheme.
The qualitative features of the intraseasonal variability were evaluated using MJO diagnostics. Wavenumberfrequency spectra of tropical precipitation indicated that TK could not simulate the MJO-related power well or the coupled equatorial Kelvin waves well. SP simulated the powers better than TK, yet a weaker mixed-Rossby gravity wave was observed. Lag-cross correlations for OLR and zonal wind anomalies clearly indicated that the eastward propagating features of the MJO as simulated by TK were insufficient, but were greatly improved by SP. MJO lifecycle composites further suggested that, based on the composite, SP simulated improved MJO features in terms of strength and location of each phase.
The difference in the structure of the organized convection of the MJO was investigated using composite analysis. The composited MSE profile indicated that TK simulated an upright MSE anomaly along with MJO development, which differed from those of the observation and SP. The observed tilted MSE anomaly implied that the MSE anomaly began to increase from the lowlevel and later reached at the upper level during the MJO development, and these features were well simulated by SP. The cause of the difference in the MSE anomaly was analysed using the budget equation of the columnintegrated MSE. Among the terms in the budget equation, the largest deviation from the observation and SP for TK was found for the vertical advection term. Further analyses revealed that vertical advection of humidity in low-levels before the peak of organized convection differed between the cases. This means that environmental low-level advection caused by shallow convection was lacking in TK and was successfully simulated in SP. The difference found for SP is consistent with previous results and the characteristics of the scheme; this indicates that the scheme is able to adequately simulate shallow convection even though deep convection exists, an advantage originating from its initial design (Baba, 2019(Baba, , 2020cBaba and Giorgetta, 2020). Tiedtke scheme in ICON-A failed to simulate MJO because of overestimated entrainment of shallow convection (Baba and Giorgetta, 2020), and the spectral scheme could suppress the overestimation. In this study, the spectral scheme provided similar improvement such as the scheme adequately simulated properties of shallow convection, but the effect appeared as correctly simulated pre-existing and coexisting shallow convection during the MJO development.
Finally, statistical features of simulated MJO that may affect seasonal prediction skill were examined. The temporal variations of the MJO index implied that PC-based MJO strength was partly captured by the present model, but additional initialization for the atmospheric components was considered necessary for better agreement with the observation. An OLR anomaly-based analysis revealed that TK partly captured some MJO events while SP captured most of the MJO events though the scheme overestimated the number of the events. Moreover, the OLR anomaly simulated by SP appeared at the correct locations, and the frequency for the OLR anomaly versus longitude was better than TK when compared with the observation.
In conclusion, SP outperformed TK in terms of intraseasonal variability from various perspectives, with comparable climatological error and seasonal variability for long-range simulations. Although SP showed slightly worse climatological (mean state) error, but its fidelity for the anomalous fields were better than that of TK, as seen in Baba and Giorgetta (2020). In particular, the weak eastward propagation, weak strength and incorrect locations of MJO events simulated by TK were improved by SP. This reveals that the known MJO-mean state tradeoff problem (i.e., improving the MJO sometimes results in degrading climatology, Kim et al., 2011Kim et al., , 2018 is not present when employing the spectral scheme (other successful case can be seen in Bechtold et al., 2008), also indicating the merit of introducing the spectral scheme in SINTEX-F2.
Further improvement in the seasonal prediction skill, which is particularly affected by the intraseasonal variability, is therefore expected when using the spectral scheme. However, improved MJO simulation performance does not necessarily result in improved prediction performance as noted by Klingaman et al. (2015), who concluded that this might be due to the nonlinear interaction between the model error in the fidelity of MJO and the initial value error. Nonetheless, while this does not guarantee successful improvement in MJO prediction for the present system, the spectral scheme is expected to outperform the original scheme, since seasonal prediction is conducted over a longer period (2-5 months, Doi et al., 2016) than the MJO prediction (2-4 weeks, Kim et al., 2018). Furthermore, long-range statistical variability induced by the MJO is important for predicting the evolution of interannual variability (e.g., Tang and Yu, 2008;Rao et al., 2009;Zhang, 2013;Baba, 2020b), which is closely related to the fidelity of the MJO in climate simulations. Also, extended range prediction might become possible through the teleconnections by improved MJO simulation (e.g., Vitart and Molteni, 2010). Although the spectral scheme outperformed the original scheme, additional initialization (assimilation) techniques for the atmosphere may be necessary for further improvement of seasonal predictions. This is because the present simple SST nudging scheme indirectly initializes the atmosphere and is considered insufficient to predict the phase of MJO as seen in the comparison of the MJO index. Therefore, the initialization of the atmosphere using an additional assimilation technique (e.g., nudging scheme or three-dimensional variational data assimilation, 3DVAR) is considered necessary. Evaluation of the change of seasonal prediction skill by the spectral scheme when using different initialization techniques will be conducted in a future study.