Phytoplankton and chlorophyll (Chla) can affect the absorption of shortwave radiation and thus the vertical distribution of heat in the upper ocean (Lewis et al., 1983). Some previous works already showed that seasonal cycle of Chla is important for modifying and reducing the bias of annual mean and seasonal cycle of shortwave absorption and, especially in the tropics, where the thermocline depth is shallow, the Chla concentration is large, and the shortwave radiation is strong (Nakamoto et al., 2001; Sweeney et al., 2005; Löptien et al., 2009).

Some climate models, including the SINTEX-F2 CGCM (Masson et al., 2012), use the satellite-observed climatology of chlorophyll to compute the shortwave absorption in the upper ocean. However, interannual variations in Chla are important as some previous works demonstrated that interannual variations of Chla could dampen El Niño/Southern Oscillation (ENSO) variability by 10–40% in CGCM experiments and hybrid coupled physics-biogeochemistry models (Timmermann and Jin, 2002; Jochum et al., 2010; Kang et al., 2017; Tian et al., 2021). Those studies also presented possible processes; during El Niño, equatorial upwelling and therefore phytoplankton nutrient supply is reduced and the mixed layer depth is increased, which can reduce Chla concentrations and associated biological heating (Figure 1). Reduced Chla also could deepen the mixed layer depth via reduced surface warming and subsurface cooling. Tian et al. (2021) also showed that interannual Chla anomalies could modulate vertical redistribution of penetrative shortwave radiation between the mixed layer and subsurface layers and affect ocean stratification and vertical mixing, which weaken the thermocline and Ekman feedback related to ENSO. Park et al. (2018) showed that the Chla-induced thermodynamic response is quicker compared to the wind-induced dynamical response, and thus can work as a precursor of ENSO using a long-term global Earth system model simulation incorporating marine biogeochemical processes. A further study of such physical-biogeochemical feedbacks via Chla may newly find room for improvement in processes and predictions associated with ENSO (Lee et al., 2014; Zhang et al., 2018; Tian et al., 2020).

To the best of our knowledge, an exploration of possible impacts of the interannual variations of Chla on seasonal predictions by a dynamical prediction system has not yet been presented, thus it is the focus of this study. We evaluated such a response by using the SINTEX-F2 prediction system, which has demonstrated its outstanding performance of predicting ENSO, Indian Ocean Dipole, and other tropical climate phenomena (see Section Methods).

We begin by exploring historical time series of the Niño3.4 Index, which is defined as SST anomalies averaged over the central-western tropical Pacific (170°W−120°W, 5°S−5°N). At a first glance, the original F2 system seems to be skillful in predicting most of significant events at least one season ahead, though the onset stage and the amplitude of the events are not well predicted (Figure 2A). Contrary to our expectation, we could not find any significant differences between the F2 and the F2chT. This may be due to the fact that ENSO is not a dominant factor for the interannual variations of Chla in the Nino3.4 region. The Nino3.4 region is one of major high-nutrient/low-chlorophyll (HNLC) regions in the world (Coale et al., 1996), where the abundance of phytoplankton is low owing to the unavailability of iron (Behrenfeld et al., 2006). In contrast, a relatively high concentration of iron is observed in the western equatorial Pacific (Ryan et al., 2006). Therefore, the regressed spatial patterns of interannual Chla anomalies on the Nino3.4 index in observations showed that the maximum of the negative Chla anomalies associated with El Niño is found mainly over the eastern edge of the Western Pacific Warm Pool, which corresponds to the separation between warm, low-salinity waters of the warm Pool and cold, high-salinity upwelled waters of the Pacific cold tongue in the central-eastern equatorial Pacific (Zhang et al., 2018). Actually, based on a horizontal map of the differences in the prediction skills of SST in June–August between the twin experiments (Figure 3), we found the statistically significant skill improvement in the eastern edge of the Western Pacific Warm Pool (140°E−160°E, equator −10°N) above the 95% confidence levels on a paired t-test. In addition, we can find the skill degradation in the northern central tropical Pacific region. The SINTEX-F2 model has an underestimation bias of ENSO amplitude (Figure 2A). The bias may be enhanced via the damping effect of interannual variations of Chla. Although we should be careful of this skill degradation, the drop in skill is not so large relative to the skill improvements in other regions. Therefore, we have not discussed much about this minor skill degradation and focused more on the skill improvement part in this study.

Comparing the observed and predicted time series of the SST averaged in the target region, we found that the difference between the F2 and F2chT was striking in June 2015, which was statistically significant above the 90% confidence levels on a paired t-test (Figures 2B, 4A,B). The observation shows development of negative SST anomaly in June 2015, which was associated with occurrence of a super El Niño (Chen et al., 2017). Different from the observed anomalies, the ensemble mean prediction by the F2 system shows development of positive SST anomaly from January 1st, 2015, which reached the maximum in June 2015. That unrealistic development of the positive SST anomaly was toned down by the F2chT system.

The F2chT still could not capture well the development of the negative SST anomaly as seen in the observation. The observed cooling process in the target region is mainly due to the air-sea coupling dynamical development processes of the 2015 super El Niño (Chen et al., 2017; Ineson et al., 2018). However, the positive SST anomaly in the Nino3.4 region, an important aspect of the occurrence of the 2015 super El Niño, was not predicted well (Figure 2A). The deficiency might be due to the so-called spring prediction barrier of ENSO events (Latif et al., 1998; Ren et al., 2016; Behera et al., 2021). Actually, the SINTEX-F2 model captured the occurrence of the 2015 super El Niño after May 2015, although the amplitude was underestimated in the ensemble mean (Doi et al., 2019a). A strong westerly wind burst activity in 2015 may be important for successful prediction of the 2015 tropical Pacific condition (Ineson et al., 2018). However, the SINTEX-F2 as well as many climate models in the world have significant biases in their representation of the westerly wind burst activity (Tan et al., 2020). We need further efforts to improve the intraseasonal variability in the model (Baba, 2021), which is beyond the scope of this study.

The differences between the two re-forecast experiments disappeared when the model was initialized on July 1^st, 2015 (Figure 3B). The SINTETX-F2 ensemble mean has an amplitude underestimation bias of the ENSO development for both experiments (Figure 3). We need further analysis from a viewpoint of probabilistic prediction as well as deterministic prediction for the rare and extreme case dynamics by conducting large ensemble members' re-forecast experiments (Doi et al., 2019a).

We could also see slight improvement in the target region for the latter half of 2009, which is related to an El Niño event. However, for other ENSO events (El Niño; 2002/03, 2006/07, La Niña; 1999/2000, 2007/08, 2010/11 2020/2021), we could not find any significant differences. The improvement of the skill over the target region shown in Figure 3C is mostly due to the 2015 super El Niño event.

The time series of the mixed-layer heat budget analysis shows that the contribution from the surface heat flux was dominant to the difference between the F2chT and the F2 during March 15, 2015–April 30, 2015 (Figure 5C), when the balance between the warming by the meridional advection and the entrainment and the cooling by the vertical diffusion played an important role on the heat budget in the target region by the F2 and F2chT (Figures 5A,B). Then, the cooling by the meridional advection and the vertical diffusion was dominant in May 2015 for the difference between the F2chT and the F2. Although the interpretation of the mixed-layer heat budget was too difficult to exactly specify the reason for the differences in the simulations, we focused on the early stage during March 15, 2015–April 30, 2015. A horizontal distribution of the mixed-layer heat budget averaged in March 15, 2015–April 30, 2015 shows contributions of several key processes; the zonal and meridional advections, vertical diffusion, as well as surface heat flux (Figure 6). For the surface heat flux contribution in April 2015 over the target region, the latent heat loss associated with the enhanced trade wind was dominant (Figures 7A–D). The enhanced trade wind seems to be related to the colder SST in the east of the target region (Figure 7E), where the mixed layer depth was also increased (Figure 7F). The deeper mixed layer depth at that region could also be due to the reduced Chla (Figure 7G) via reduced surface warming and subsurface cooling. This is partly consistent with the damping process of El Niño via Chla discussed by the previous work (Jochum et al., 2010). The deeper mixed layer depth may propagate westward and extend to the target region as Rossby waves. Moreover, the difference in the SST prediction between the F2 and the F2chT could contribute to a difference in the rainfall prediction (Figure 4C). Therefore, considering the interannual variation of Chla is important for predicting not only the oceanic conditions but also atmospheric conditions via air-sea coupling in the tropical Pacific.

The F2chT predictions show larger ensemble spread relative to the F2 prediction (Figure 4A). Interestingly, 9 out of 12 members of the F2chT show cooling tendency, while 3 ensemble members show warming tendency. It might be interesting to compare the difference between the warmest 3 ensemble means of the F2chT and the F2 ensemble mean with that of the other 9 ensemble means of the F2chT and the F2 ensemble mean (Figure 8). Although the other 9 ensemble means are similar to the F2chT ensemble mean, the warmest 3 ensemble means are almost a mirror image of the F2chT ensemble mean; the time series of the mixed-layer heat budget analysis shows that the contribution from the surface heat flux was dominant to the difference between the warmest 3 ensemble means and the other 9 ensemble means during February 1, 2015–March 31, 2015 (Figures 8C,D). The latent heat loss, associated with the weaker trade winds, was the dominant contributor to the net surface heat loss during February–March 2015 over the target region (Figures 8E,G,I). The mixed layer depth was also decreased (Figure 8K). Those larger uncertainties among the ensemble members of the F2chT may be due to the shallower mixed-layer depth during January–February 2015 (Figure 9A), which could be due to the increased Chla (Figure 9B) via enhanced surface warming and subsurface cooling.

In order to present a skill assessment from a view point of impacts of interannual variations of Chla on seasonal predictions, we conducted twin re-forecast experiments with the SINTEX-F2 dynamical seasonal prediction system: one system used the observed climatology of Chla to compute the shortwave absorption in the upper ocean, while the other used the observed Chla with year-to-year (interannual) variations. Although the interannual Chla impacts on predictions of SST in the central-eastern equatorial Pacific were limited, improvements in predictions of SST were found over the eastern edge of the Western Pacific Warm Pool, in particular for the 2015 super El Niño year. The results showed that considering the interannual variation of Chla in a dynamical seasonal prediction system is potentially important for improvement of seasonal climate predictions.

Several mixed-layer processes are found to play a role in the improvements in the SST predictions. Those seem partly related to the damping process of El Niño via Chla. However, the exact contributions of those processes are complicated to be fully understood at this stage.

Several factors could affect the results. For example, although this study used the satellite-observational data of surface Chla concentration, the maximum Chla concentration is actually observed at the subsurface ocean in the tropical Pacific (Le Borgne et al., 2002; Lee et al., 2014; Yasunaka et al., 2021). Developing a sustainable ocean observing system with Biogeochemical (BGC) Argo (Bittig et al., 2019) will shed light on 3-dimentional variations of Chla, which could not be observed only by the satellites. Along with the development of the observational system, a new scheme to compute the shortwave absorption in the upper ocean based on the vertical profile of Chla should be developed. For example, Manizza (2005) uses the entire vertical profile of Chla to compute the induced biological heating at each vertical level. Such a parameterization might enhance the impact of Chla in the tropical Pacific (Park et al., 2014). In addition, a seasonal prediction system based on a global Earth system model to resolve physical-biogeochemical feedbacks via Chla may open a new door to discovery of predictability of variations in the marine system associated with the ENSO, including ocean warming, acidification, deoxygenation, biological production, and biodiversity (e.g. Park et al., 2019). On the other hand, we should be careful of possible skill degradation by those approaches. For example, Lim et al. (2018) demonstrated that an overestimation bias of mean Chla state in the tropics in their earth system model enhanced cold SST bias in the tropical Pacific. Further studies in such research streams are necessary. The tropical Pacific may serve as an optimal test-bed for those studies, because the interannual dynamical variations associated with the ENSO are already predictable with high accuracy relative to other basins.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

TD performed the seasonal prediction experiments, analyzed the observation data, and model prediction outputs. Both authors contributed to designing the research, interpreting results, and writing the manuscript.

This research was supported by JSPS KAKENHI Grant 20K04074.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.