The dynamical seasonal prediction system is based on an fully coupled global ocean–atmosphere circulation model (CGCM) called SINTEX-F ver. 2 developed under the EU–Japan collaborative framework (Masson et al., 2012; Sasaki et al., 2013), which has a dynamical sea-ice model, together with a higher-resolution relative to its previous version (Luo et al., 2003, 2005; Masson et al., 2005). This system adopts a relatively simple initialization scheme based only on the nudging of the SST data (Doi et al., 2016) and a three-dimensional variational ocean data assimilation (3DVAR) method by taking three-dimensional observed ocean temperature and salinity data into account (Doi et al., 2017). In consideration of the uncertainties of both initial conditions and model physics, the original system had 12 members for the predictions initiated on the first day of each month. The 12-member system has recently been upgraded to a bigger ensemble system of 108 members using the Lagged Average Forecasting (LAF) method (Doi et al., 2019, 2020b). Based on this new system, we conducted the reforecast experiments from the nine initialized dates (1st−9th) in every month of 1983–2020. The prediction anomalies were determined by removing the model climatology at each lead-time and ensemble member using the reforecast outputs over the period 1983–2015.

To evaluate the prediction results, we used the NCEP/NCAR reanalysis data (Kalnay et al., 1996) and the NOAA Interpolated Outgoing Longwave Radiation (OLR), for atmospheric variables. The monthly climatologies of these datasets were calculated by averaging the monthly data from 1983 to 2015, and then the anomalies were derived through deviations from those climatologies. The linear trends were removed for all of the analyses.

In this study, the correlation coefficient and root-mean-square error (RMSE) between the anomaly of prediction and the verification value based on the observation/reanalysis datasets were adopted as the deterministic prediction skill scores. The p-value was calculated based on a paired t-test. We also assessed the potential predictability by dividing the predicted variability into signal (S) and unpredictable noise (N) components. S is the common feature among the ensemble members, and could be the predictable component of variability mainly due to tropical SST forcing, which was estimated as the ensemble mean of the 108 members. N was estimated as the ensemble spread that provided a measure of uncertainties among the predictions arising from the random effects of chaotic atmosphere-ocean variability. By calculating the SNR, we measured robustness of the prediction of the signal as well as the potential predictability. When the ensemble mean is small relative to the ensemble spread, the prediction is potentially difficult due to the large uncertainty.

Anomalies among the ensemble members (defined as deviations from the ensemble mean) and their co-variability in the ensemble space may provide useful insights into possible precursors and teleconnection patterns related to a climate event (Ma et al., 2017; Ogata et al., 2019; Doi et al., 2020a,b,c). Here, we calculated an inter-ensemble correlation among the 108-member ensemble predictions in target months. The Pearson correlation measures how two continuous timeseries co-vary over time and indicate the linear relationship as a number between – 1 and 1. In contrast, inter-ensemble correlation measures how two predictions co-vary over ensemble members by adding the ensemble phase space to the conventional time dimension. Here, we calculated it as

Here nt is the number of target years, ne is the ensemble size (108), X(t, e) is predictions of the target index [area-averaged anomalous temperature, e.g., monthly surface air temperature regionally averaged over most part of Japan (129°-141°E, 30°N−37°N) for Figures 3, 5; Niño3.4 index (SST anomaly averaged in 170°W−120°W, 5°S−5°N) for Figure 4] and a function of time (t) and ensemble member (e), X¯ is the average of X(t, e) over target years and all ensemble members, Y(x, y, t, e) is predictions of the target variables (anomalous temperature, OLR, geopotential height (GH) at 200 hPa, and velocity potential at 200 hPa) and a function of two-dimensional horizontal (x, y), time (t) and ensemble member (e), Y(x,y)¯ is the average of Y(x, y, t, e) over target years and all ensemble members, respectively.

Figure 1 shows seasonality of skills of 1-month lead prediction (e.g., prediction of December issued on 1st−9th of November) of monthly surface air temperature regionally averaged over most part of Japan (129°-141°E, 30°N−37°N). Based on the correlation skill, temperature in September is the most predictable (the skill: ~0.7), then, that in May and October are the second (the skill: ~0.5) and the third most predictable (the skill: ~0.4), respectively (Figure 1A). September is the only month in which the model prediction skill beats the persistence. The correlation skills for the other months are not statistically significant beyond the 95% confidence levels on the paired t-test.

The reanalysis data shows the seasonality of the standard deviations of the averaged air temperature; large in wintertime and small in summertime (Figure 1B). It is not captured well by the ensemble mean prediction. Because atmospheric intrinsic variability included in the observation is largely removed in the ensemble mean, this suggests that the model is not skillful at predicting the amplitude. Interestingly, the ensemble spread shows a strong seasonality, large in wintertime and small in summertime, similar to what is seen in the standard deviation of the reanalysis data. The large standard deviation in the reanalysis data and the large ensemble spread in the predictions for wintertime can be partly due to strong atmospheric internal variability such as blocking events (e.g., Takaya and Nakamura, 2005). Note that a relatively large ensemble spread is also seen in June (Figure 1B), when the Baiu season onsets (e.g., Tomita et al., 2011). For other summer months, the ensemble spread is relatively small, in particular for August (~0.5°C), which is almost half of that in wintertime. In August, the Bonin High is formed around Japan, which is relatively stationary, as a separate system from the Tibetan and North Pacific anticyclones (Neyama, 1968; Enomoto et al., 2003). The Bonin High is also stable in some years, keeping aside the interannual variations, associated with limited propagations of stationary Rossby waves along the Asian jet (Enomoto, 2004). Those features may explain the small ensemble spread in August. Figure 1B also shows the relationship between the RMSE and the ensemble spread. The differences among the ensemble members are larger than the difference between the ensemble mean prediction and the observed value in September. In other words, the uncertainty of the prediction is overestimated (Fortin et al., 2014). In wintertime prediction, the opposite occurs, where the model underestimates the sources of error. The problem is common to most dynamical seasonal prediction systems based on GCMs (Scaife and Smith, 2018). Those errors in the prediction of the amplitude can be corrected by a postprocessing step that would amplify the predictable signal based on the correlation skills (Eade et al., 2014).

The SNRs are low in all calendar months (0.3–0.4) except for those in August and September (< 0.4) (Figure 1C), when the correlation skill and the RMSE are also relatively high. This suggests that the SNR can partly indicate the seasonality of the prediction skill.

Interestingly, we found the correlation skill of December becomes high and beats the persistence when we selected years in which the SNR is larger than the 1983–2020 mean of the SNR (Figure 1A). Therefore, we will focus on December as well as the most predictable September hereafter. Note that the SINTEX-F prediction system is not skillful in the 2-month lead prediction for all months.

At a first glance, the model seems to be skillful in predicting most of the significant events from the timeseries for September (Figure 2A). The correlation skill is 0.72. Figure 2A shows that the SNR has substantial interannual variability (0.1–1.4), although the ensemble spread is almost constant. In high SNR years when SNR is higher than its 1983–2020 mean, the correlation skill is higher (~0.9) compared to the values considered with all the years (~0.7) (Figure 2B). Even for the low-SNR years, when SNR is lower than its 1983–2020 mean, the correlation skill is still about 0.4, which is comparable to that of the persistent skill.

We selected warm years when the anomaly in the reanalysis data was above +0.5°C (1998, 1999, 2007, 2010, and 2019) and cold years when the anomaly in the reanalysis data was below – 0.5°C (1987, 1992, 1993, 1996, 1997, 2009, 2014, and 2015), then analyzed the inter-ensemble correlations. For the warm years, ensemble members that predicted a La Niña Modoki-like state (Ashok et al., 2007) mostly tend to predict an intensely warmer-than-normal condition over Japan (Figure 3A). They also tend to predict an enhanced convection (Figure 3C) and a positive SST anomaly (Figure 3A) over the Philippine Sea around 20°N and a positive GH200 anomaly over Japan (Figure 3E). Note that for Figure 3C, negative correlation (blue color) means that ensemble members that predicted an intense warmer-than-normal condition over Japan mostly tend to predict enhanced convection. The ensemble co-variability may suggest the teleconnection from the Philippine Sea into Japan via the enhanced local Hadley circulation associated with the La Niña-Modoki-like state. The suppressed (enhanced) convective activity over the western/central equatorial Pacific (tropical Indian Ocean) (Figure 3C), and the positive (negative) velocity potential anomaly at 200 hPa over the western Pacific and the tropical Indian Ocean (the eastern Pacific) (Figure 3G) suggest an enhancement of the Walker circulation in the Indo-Pacific basin, which also may contribute to the positive temperature and GH200 anomaly over Japan. The inter-ensemble correlations suggest that the La Niña Modoki could work as a potential source of the 1-month predictability of temperature over Japan. This may be partly consistent with the La Niña Modoki developing in autumn (Zhang et al., 2011). Note that the negative GH anomalies at 200 hpa associated with the cold SST anomaly are also seen in the central and western tropical Pacific (Figures 3A,E), which are similar to a typical La Niña Modoki state. However, they are not simply linked with the velocity potential at 200 hPa (Figure 3G), which looks similar to the canonical La Niña-like state. The ensemble co-variability may suggest a pattern resembling a mixture of a La Niña Modoki-type and a canonical-type (Karnauskas, 2013).

Interestingly, we could not find the ensemble co-variability between prediction of the tropical Pacific condition and prediction of atmospheric variables over Japan for the cold years (Figures 3B,D,F,H). The cold years include the super El Niño years: 1997 and 2015, then we also focused on those 2 years (Figure 4). As similar to Figures 3B,D,F,H, we could not find the ensemble co-variability between prediction of the super El Niño and prediction of atmospheric variables over Japan, although it suggested a weakening of the Walker circulation in the tropical Pacific (Figures 4G,H). Since the El Niño-Modoki and its impacts on East Asian summer are different from those of ENSO (Weng et al., 2007), further studies are necessary for deeper understanding of the uncertainty in their predictions.

In stark contrast to September, the model seems to be not skillful in predicting most of significant events in December at a first glance (Figure 2C). However, when we selected 15 years in which the SNR is higher than the 1983–2020 mean of the SNR (1984, 1987, 1988, 1989, 1990, 1994, 1997, 2000, 2004, 2006, 2008, 2009, 2010, 2017, 2019), the correlation skill becomes 0.67, while that in the other 23 years it is – 0.15 (Figure 2D). The difference is statistically significant beyond the 99% confidence levels based on the t-test.

For the warm years, ENSO-Modoki contributions are also found as was seen in September, but with the opposite sign; warmer-than-normal conditions over Japan are linked with the El Niño-Modoki like condition in the tropical Pacific (Figure 5A). Ensemble members that predicted an El Niño Modoki-like state mostly tend to predict a suppressed convection (Figure 5C), a negative SST anomaly (Figure 5A) and a negative velocity potential anomaly (Figure 5G) over the Philippine Sea and the Maritime Continent, and a positive GH200 anomaly over Japan (Figure 5E). It may suggest the weakened local Hadley circulation in the western Pacific associated with the El Niño Modoki. The ensemble co-variability are also consistent with a possible teleconnection pattern from suppressed convection over the western Pacific associated with El Niño/El Niño Modoki into warm winters in Japan via an emanation of the Rossby wave train in the upper troposphere as discussed by previous works (Sakai and Kawamura, 2009; Ueda et al., 2015; Kuramochi et al., 2021).

In addition, ensemble members that predicted warmer conditions in the western tropical Indian Ocean tend to predict warm temperature around Japan (Figure 5A). They also tend to predict an enhanced convection in the western tropical Indian Ocean (Figure 5C), a negative GH anomaly over southern China and a positive GH anomaly in the upper troposphere over Japan (Figure 5E). The ensemble co-variability suggests a possible teleconnection from the western tropical Indian Ocean via stationary Rossby waves along the subtropical jet. Although this is partly seen by the wave activity fluxes at 200 hPa (Takaya and Nakamura, 2001; Morioka et al., 2014) for the warm years (Figure 5E), we could not find the significant ensemble correlations between prediction of the temperature around Japan and prediction of the wave activity flux along the subtropical jet. Doi et al. (2020b) demonstrated the propagation pathways of stationary Rossby waves from the western tropical Indian Ocean to Japan for the record-breaking warm winter case over Japan in 2019/2020 by the SINTEX-F, when the positive SST anomaly in the western tropical Indian Ocean was observed after the super Indian Ocean Dipole (IOD) (Doi et al., 2020a). As a result, the southern penetration of cold air masses into Japan from high latitudes was weakened, causing warmer winter.

Interestingly, we could not find the ensemble co-variability between prediction of the tropical Pacific and Indian Ocean condition and prediction of atmospheric variables over Japan for the cold years (Figures 5B,D,F,H). Further studies are necessary for deeper understanding of such asymmetries.