Introduction

The El Niño–Southern Oscillation (ENSO) is the most prominent interannual variability mode in the global climate system. The preference of the ENSO peak phase in boreal winter is known as the ENSO phase-locking phenomenon and is one of the fundamental features of ENSO1,2,3,4. The solar-forced seasonal cycle in the tropical climate system modulates both the growth and phase transition, and this seasonal modulation is responsible for making El Niño and La Niña peak preferentially during the boreal wintertime5 (Supplementary Fig. 1a). Several studies demonstrated that the mechanism of ENSO’s winter phase-locking is dominated by the seasonal modulation of the ENSO’s sea surface temperature (SST) growth rate5,6,7,8,9,10. The instability is strongest in the boreal summer to autumn and weakest in the boreal spring, which is controlled by the zonal advective and thermodynamic feedbacks, leading to the maximum ENSO sea surface temperature anomalies (SSTA) tending to occur at the end of the calendar year11. However, previous studies of the ENSO winter phase-locking mainly focused on the SSTA in the Niño3 or Niño3.4 region. The SSTA phase-locking phenomenon in the far eastern equatorial Pacific (FEP) and coastal regions is entirely different from that in the central-eastern Pacific. Although El Niño was named because fishermen in the late 19th century noticed a counter-current appear off the coast of Peru around Christmas time, the SSTA maximum at FEP or Peru coast usually occurs during the boreal summer from April to July (Supplementary Fig. 1b, c).

Rasmusson and Carpenter found that the warming SSTA observed near the Peru coast typically occurred prior to the development of warming SSTA in the central equatorial Pacific, particularly before 19801. This observation suggests that the coastal SSTA can serve as a valuable precursor for predicting subsequent basin-scale ENSO warming. However, this relationship did not hold during the strong El Niño event of 1982/8312, and it appears to be less common in subsequent events. In general, these warming SSTA events in the FEP and South American coast regions could be followed or preceded by warming SSTA in the central equatorial Pacific during El Niño events1,12,13,14, as shown in Supplementary Table 1. However, the warming SSTA in the FEP and Peru coast regions events may also occur independently of basin-scale SSTA warming15. To distinguish the warming SSTA events occurring in the FEP and South American coast regions from conventional El Niño events characterized by warming in the central-eastern Pacific, these coastal warm events are commonly referred to as coastal El Niño.

It is now widely recognized that the conventional El Niño event, characterized by warming in the central-eastern Pacific, is primarily caused by westerly wind anomalies in the central Pacific through the adjustments of equatorial waves. However, there are still discrepancies in the literature about what are the main driving mechanisms of the coastal El Niño phenomenon. According to observational analysis, Hu et al. suggested that coastal El Niño events differ significantly in evolution, mechanism, and timing16. These coastal El Niño events can arise from various factors, including the equatorially centered intertropical convergence zone (ITCZ) during extreme basin-scale El Niño events decaying phase, downwelling Kelvin waves driven by westerly wind bursts in the western Pacific, and westerly surface wind anomalies in the eastern equatorial Pacific. Based on in situ observations, Takahashi and Martínez argued that the strong northerly winds across the equator and southward shift of ITCZ are possible to initiate coastal warming such as the coastal El Niño event in 192515. Peng et al. also emphasized the important role of Kelvin waves generated by westerly winds in the equatorial Pacific and coastal wind anomalies in the 2017 coastal El Niño event through observations and ocean model experiments17. While previous studies extensively investigated the evolution of coastal El Niño events through dynamical and numerical analyses13,14,15,16,17,18,19,20, the knowledge and research in understanding aspects of SSTA in the far eastern Pacific and coastal regions are relatively limited compared to the extensive and systematic studies of SSTA in the central-eastern Pacific region. There is a need for further research and understanding of the SSTA dynamics in these specific regions to enhance our overall comprehension of coastal El Niño events.

Additionally, although previous studies thoroughly investigated the characteristics and mechanisms of ENSO phase-locking in Niño3 or Niño3.4 regions, the summer phase-locking of FEP SSTA or coastal El Niño is still not well described. A general understanding of the nature and dynamics of the SSTA phase-locking at FEP and coastal regions in the observation and massive simulations of a large number of expensive climate models has still yet to emerge. Therefore, our study aims to fill this gap by utilizing observational data and comparing them with simulation results from the CMIP5 and CMIP6 models, and by adopting a conceptual stochastic model framework to investigate the characteristics and dynamics of summer phase-locking in the far eastern equatorial Pacific.

Results

Features of summer phase-locking in the far eastern Pacific

The observed SSTA has an apparent feature of phase-locking in the central-eastern Pacific where the ENSO-related SST variability is strong. As shown in Fig. 1a, both the peak phase histogram and seasonal variance of SSTA show significant boreal winter phase-locking located in the Niño3.4 region (5°S–5°N, 170°–120°W), which has been the main area of focus for previous phase-locking studies in the literature. However, the phase-locking phenomenon in the FEP (5°S–5°N, 90°–80°W) and South American west coast regions (Niño 1+2; 10°S–0°, 90°–80°W) is completely different from that in the central-eastern Pacific. In these regions, there is a preference for summer phase-locking shown both by seasonal variance (red bars in Fig. 2b and Supplementary Fig. 2a) and peak phase histogram of SSTA (red bars in Fig. 2c and Supplementary Fig. 2b), which means that the preferred peak of SSTA tends to occur in the boreal summertime from April to July. Although the preference of phase-locking is strong in both the Niño3.4 and FEP regions, their preferred peak months are during opposite phases (Figs. 1a and 2c). It is worth noting that the histogram of FEP SSTA has a second peak in the wintertime, but its probability is much lower than in summer (Fig. 2c). The most significant air-sea interaction region for ENSO is in the central-eastern Pacific, but the ENSO SSTA in that region can directly influence the SST variability throughout the rest of the Pacific, resulting in a small probability of winter peaks in the FEP or other Pacific regions.

Fig. 1: Characteristics of the equatorial Pacific SSTA in observations.
figure 1

a Peak phase histogram (bars) and seasonal variance (curve) of Niño3.4 SSTA. b Regression (red curve) and correlation coefficients (blue curve) of equatorial SSTA (averaged between 5°S–5°N) against the Niño3.4 index for different longitudes. c SSTA standard deviations of ENSO forcing (red curve) and residual SST components (blue curve) for different longitudes. In (ac), the vertical lines or shading denote the range of the minimum and maximum values between different observed datasets.

Fig. 2: Characteristics of the equatorial far eastern Pacific SSTA in observations.
figure 2

a Ensemble average of time series of FEP SSTA (gray bars) and its ENSO forcing (red curve) and residual SST components (blue curve). The peaks of FEP SSTA event are shown as circles. b Seasonal variance and c peak phase histogram of FEP SSTA (red bars) and its residual SST component (blue bars). In (ac), the vertical lines or shading denote the range of the minimum and maximum values between different observed datasets.

Figure 1b shows the regression and correlation of equatorial SSTA (averaged between 5°S–5°N) with Niño3.4 SSTA, which are used to examine the influence of ENSO’s central-eastern Pacific SSTA on the rest of the equatorial Pacific. The regression and correlation coefficients are large and close to one around the Niño3.4 region but decrease in the eastward and westward directions. The coefficients drop to zero near 160°W and become negative values in the western Pacific, indicating a negative relationship between ENSO-related SSTA in the central-eastern and western Pacific. In the far eastern Pacific region (90°–80°W), the correlation and regression coefficients decrease to 0.7 and 0.6–0.8, respectively, suggesting that the ENSO forcing is weakened here and the residual SSTA variability would be as important as the ENSO forcing. To further investigate the relative roles of ENSO forcing and residual SSTA variability, the equatorial SSTA at different longitudes are separated into an ENSO-related component and a residual SSTA component by a linear regression model:

$$T\left(x,t\right)=E\left(x\right){T}_{N34}\left(t\right)+{T}_{R}\left(x,t\right)$$
(1)
$${T}_{R}\left(x,t\right)=T\left(x,t\right)-E\left(x\right){T}_{N34}\left(t\right)$$
(2)

where E is the regression coefficient of equatorial SSTA (T; averaged between 5°S-5°N) at different longitude (x) against the Niño3.4 index (\({T}_{N34}\); 5°S–5°N, 170°–120°W). The residual SSTA (\({T}_{R}\)) is defined as the difference between T and \(E\,\cdot \,{T}_{N34}\) as described in Eq. (2).

Figure 1c shows the amplitudes of ENSO forcing (\(E\,\cdot\, {T}_{N34}\)) and residual SSTA (\({T}_{R}\)) components of equatorial Pacific SSTA, as the first and second terms in Eq. (1), respectively. In the central-eastern Pacific around 170°–110°W, the ENSO forcing component (red curve) is significantly larger than the residual SSTA component (blue curve), indicating that ENSO-related SSTA is dominant in this region and that the residual SSTA process is minor. In contrast, in the FEP region (around 90°–80°W), the residual SSTA component is as large as the ENSO forcing component, suggesting that the process associated with residual SSTA is important there.

The time series of FEP SSTA and its ENSO forcing and residual SSTA components are shown in Fig. 2a. Although the FEP SSTA is composed of the \({T}_{R}\) (blue curve) and \(E\,\cdot\, {T}_{N34}\) (red curve), the peak of FEP SSTA events is usually in the same phase with \({T}_{R}\) instead of \(E\,\cdot\, {T}_{N34}\). This feature is also reflected in the seasonal variance and phase histogram of FEP SSTA and its residual SSTA in Fig. 2b, c. After removing the direct ENSO forcing from total FEP SSTA, only the residual SSTA remains. The variance and peak phase preference in the boreal winter are significantly reduced but are still evident in summer (blue bars in Fig. 2b, c). This result suggests that the summer variance and peak preference of the FEP SSTA are contributed by the residual SSTA component, while the ENSO forcing provides only a small probability of a winter peak.

The strong correlation (R = 0.98) between the observed ensemble average FEP and Niño 1 + 2 SSTA, as well as their similar phase-locking characteristics, suggests the presence of common underlying background fields or physical processes that contribute to the SSTA phase-locking across the South American coast, Nino 1 + 2, and FEP regions. Our analysis exclusively focuses on the FEP SSTA but the results are consistent with those obtained from Niño 1 + 2 SSTA analysis.

Dynamics of summer phase-locking in the far eastern Pacific

In the previous section, we found that the summer SSTA phase-locking in the FEP region is determined by the residual SSTA process. To further investigate the dynamics of the residual FEP SSTA, we adopt a linear stochastic-dynamical model (SDM) for the residual FEP SSTA as follows:

$$\frac{d{T}_{R}}{{dt}}=\alpha \left(t\right){T}_{R}\left(t\right)+\beta \left(t\right){T}_{N34}\left(t\right)+f\left(t\right){h}_{w}\left(t\right)+\sigma \left\langle \xi \left(t\right)\right\rangle$$
(3)

As described in “Methods,” the terms on the right-hand side of Eq. (3) represent the seasonal damping caused by the residual SSTA with a damping rate α, the influence of ENSO-related SSTA with coefficients β, and thermocline anomaly with coefficients f. The direct ENSO forcing for the residual FEP SSTA (\({T}_{R}\)) was removed by the linear regression model, as shown in Eq. (2). However, there may still be an indirect ENSO forcing that might influence \({T}_{R}\). For example, the Niño3.4 SSTA-induced wind change can generate oceanic waves that may affect FEP SSTA with a delay. For this reason, our SDM still includes the ENSO forcing, where the \(\beta {T}_{N34}\) and \(f{h}_{w}\) terms reflect the ENSO’s impact on \({T}_{R}\) with consideration of the response timescale of \({T}_{R}\). More information regarding the suitability of SDM for residual FEP SSTA can be found in Supplementary Note 1.

After providing the seasonally varying parameters with artificial noise (random noise) as well as the observed Niño3.4 index and western Pacific thermocline anomalies, the SDM-simulated residual FEP SSTA (\({T}_{R}\)) is estimated by integrating Eq. (3). Total FEP SSTA (T) can be obtained by adding back the Niño3.4 index as in Eq. (1). As shown in Fig. 3a, b, the SDM simulation can perfectly reproduce the seasonal variance and phase histogram of total FEP SSTA (T) and residual FEP SSTA (\({T}_{R}\)) for observations. This indicates that the linear stochastic-dynamic model explains the annual mean of amplitude, seasonal variation of amplitude, and phase-locking performance of FEP SSTA very well.

Fig. 3: Features of FEP SSTA phase-locking for SDM simulations fitted to observations.
figure 3

a Seasonal variance and b peak phase histogram of FEP SSTA (red bars) and its residual SST component (blue bars). The peak phase histogram of residual FEP SSTA (bars) for c SDM-α simulations, d SDM-β simulations, and e SDM-f simulations. In (ce), the red curves with dots denote the annual mean-removed α, β, and f, and red circles on the right y-axis indicate the annual mean of α, β, and f. In (ae), the vertical lines or shading denote the range of the minimum and maximum values between different observed datasets.

To further discuss the individual contribution of each SDM parameter to the residual FEP SST phase-locking, we adopt three additional SDM simulations: (1) only considering the seasonal modulation of SSTA damping rate α (SDM-α; set β and f as annual mean); (2) only considering the seasonal modulation of ENSO-related SSTA coefficient β (SDM-β; set α and f as annual mean); and (3) only considering the seasonal modulation of ENSO-related thermocline coefficient f (SDM-f; set α and β as annual mean). When only the seasonality of the SSTA damping rate is considered (SDM-α), the phase histogram clearly shows a summer preferred peak month with a strong annual cycle of α, which turns from positive to negative in boreal summer (Fig. 3c). The feature of the SDM-α histogram is similar to that of the SDM histogram (Fig. 3b, blue bars), suggesting that the summer phase-locking of residual FEP SSTA is mainly controlled by the seasonal modulation of the damping process α. In contrast, for the SDM-β and SDM-f simulations that consider only the seasonal variation of the coefficients β and f, respectively, there is only very weak phase-locking (Fig. 3d, e).

It is noteworthy that although the annual cycle amplitude of β (Fig. 3d) is comparable to that of α (Fig. 3c), the contribution of β to phase-locking is weak because the absolute value of annual mean of β is smaller than that of α. The annual mean of α is about −3.5 year−1 but the annual mean of β is only −0.1 year−1 (red circles on the right y-axis), resulting in much weaker contribution of the ENSO-related SSTA term \(\beta \,\cdot\, {T}_{N34}\) to the residual FEP SSTA phase-locking. The compensation between the seasonal evolution of Niño3.4 index (winter maximum) and β (spring maximum) can also reduce the contribution of \(\beta \,\cdot\, {T}_{N34}\). However, additional sensitivity experiments (Supplementary Note 2) point out that the contribution of \(\beta\, \cdot\, {T}_{N34}\) to phase-locking remains weak even when the peak phase preference of Niño3.4 index is removed or shifted to boreal summer (Supplementary Fig. 12). The contribution of the ENSO-related SSTA process can be increased only when the absolute value of annual mean of β increases. For the ENSO-related thermocline term \(f\,\cdot \,{h}_{w}\) (Fig. 3e), its minor contribution is because both the amplitude of seasonal variation of f and absolute value of annual mean of f are small compared to that of α.

Biases of FEP phase-locking in climate models

In this section, we examine the performance of FEP phase-locking in the CMIP5 and CMIP6 models and try to identify the reasons for the failure of climate model simulations of FEP phase-locking. First, the simulated peak phase histogram and seasonal variance of Niño3.4 SSTA show a very weak boreal winter phase-locking for the CMIP ensemble and about 72% of the total models have the winter-preferred peak of the ENSO SSTA (Fig. 4a). This result suggests that, as previous studies have shown, although most models perform well in simulating the ENSO peak calendar month, most climate models fail to simulate the strength of the ENSO peak preference9,11,21,22,23,24,25,26. Figure 4b shows the influence of ENSO’s central-eastern Pacific SSTA on the equatorial Pacific SST variability. The basic features of the ensemble mean are similar to the observations (Fig. 1b). The largest regression and correlation coefficients are located around the Niño3.4 region and decrease in its eastward and westward directions, however, there are still large spreads between the models. Moreover, the coefficient drops to zero near 150°W, much further west than observed, reflecting the excessive westward extension of the ENSO SSTA in climate models26,27. In the FEP region, the regression coefficients of SSTA and Niño3.4 index are weaker than the observed values in most climate models (Fig. 4b and Supplementary Fig. 3a), indicating that ENSO has a less efficient effect on the FEP SSTA. The amplitude of the ENSO forcing component, however, is still as large as that of the residual SSTA component for FEP SSTA (see Fig. 4c and Supplementary Fig. 3b).

Fig. 4: Characteristics of the equatorial Pacific SSTA in climate models.
figure 4

a Peak phase histogram (bars) and seasonal variance (curve) of the Niño3.4 SSTA. The bottom panel shows the percent of models for different preferred peak month and the triangle indicates observed peak month. b Regression (red curve) and correlation coefficients (blue curve) of equatorial SSTA (averaged between 5°S–5°N) against the Niño3.4 index for different longitudes. c SSTA standard deviations of ENSO forcing (red curve) and residual SST components (blue curve) for different longitudes. In (ac), the vertical lines and shading denote the range of minimum and maximum values between different models, and in (a) the vertical lines and shading correspond to histogram and seasonal variance, respectively.

The performance of FEP phase-locking in the climate models is shown in Fig. 5a. In CMIP simulations, only about half of the simulations (60% of total models) have the preferred peak month of FEP SSTA in boreal summer from April to July, and the phase-locking strength of the CMIP ensemble is smaller than the observed (red bars in Fig. 5a), indicating that most climate models fail to simulate the performance of SSTA phase-locking in FEP region. After the ENSO forcing \(E\,\cdot\, {T}_{N34}\) is removed from the total FEP SSTA, the residual FEP SSTA phase-locking strength of the CMIP ensemble is still much smaller than observed (blue bars in Fig. 5a), but most models have the peak month in boreal summer (82% of total models). This demonstrates that the poor performance of ENSO’s impact and phase-locking leads to bias in the peak month of the total FEP SSTA phase-locking in the climate models.

Fig. 5: Characteristics of the equatorial far eastern Pacific SSTA in climate models.
figure 5

a Peak phase histogram of FEP SSTA (red bars) and its residual SST component (blue bars). b Same as (a) but for SDM simulations fitted to CMIP models. In (a, b), bottom panel shows the percent of models for different preferred peak month and the triangle indicates observed peak month. The vertical lines and shading denote the range of minimum and maximum values between different models.

To further investigate the reason for the weak-bias of phase-locking strength of FEP SSTA, the SDM is also applied to the simulated residual FEP SSTA. The SDM simulations can also reproduce the phase histogram of total and residual FEP SSTA in climate models (Fig. 5b and Supplementary Figs. 47), confirming that the reproduction of residual SSTA by SDM simulations can be used for observations and models, and it is a universal method. The individual contribution of each SDM parameter to the residual FEP SST phase-locking for CMIP models is shown in Fig. 6. As with the observations, the summer phase-locking of residual FEP SST is mainly dominated by the seasonal modulation of SSTA damping rate α (SDM-α) with very small contributions from the ENSO-related SSTA term β (SDM-β) and thermocline term f (SDM-f). However, as shown in Figs. 5a and 7a, the strength of the phase-locking preference of the simulated residual FEP SSTA is much weaker compared to the observed. This bias in the phase-locking strength is mainly due to the small contribution of the SSTA damping term (SDM-α) in climate models (Figs. 6a and 7a), which is not related to the ENSO-related SSTA (SDM-β) and thermocline terms (SDM-f).

Fig. 6: Features of residual FEP SSTA phase-locking for SDM simulations fitted to climate models.
figure 6

The peak phase histogram of the residual FEP SSTA (bars) and annual mean-removed α, β, and f (red curves with dots) for a SDM-α simulations, b SDM-β simulations, and c SDM-f simulations. In (ac), bottom panel shows the percent of models for different preferred peak month and the triangle indicates observed peak month. The vertical lines and shading denote the range of minimum and maximum values between different model simulations.

Fig. 7: Distribution of the residual FEP SSTA phase-locking strength in climate models.
figure 7

The boxplots (blue bars) represent the assessment of the climate models. a The phase-locking strength of residual FEP SSTA for CMIP models and different SDM simulations. b Annual cycle amplitude and c annual mean of α, β, and f for climate models. The annual cycle amplitude is defined as half of the difference between the maximum and minimum values. d The dependence of the annual mean of α on the phase-locking strength of SDM simulations in climate model sensitivity experiment (gray dots) and its boxplot for the distribution of the SDM-simulated phase-locking strength. In (ad), the red circles indicate the mean of observed values and red vertical lines denote the range of minimum and maximum values between different observed datasets. For boxplots, the outlier blue dots are values that are more than 1.5 IQR for models.

A previous study found that in climate models the weak-biased ENSO’s SSTA phase-locking in the central-eastern Pacific is due to the small annual cycle amplitude of the SST growth rate11. However, for the FEP SSTA, the annual cycle amplitudes of the SDM parameters, especially the SSTA damping rate α, are comparable to the observed values (Fig. 7b). This indicates that the annual cycle amplitude of α is not responsible for the weak-bias of residual FEP SSTA and SDM-α phase-locking. The weak phase-locking of SDM-α is primarily attributed to the much smaller annual mean of the SSTA damping rate α in the climate models (Fig. 7c). In a linear stochastic-dynamical model, with the same annual cycle amplitude of damping rate, a more stable system (i.e., a more negative annual mean of damping rate) has a weaker phase-locking strength28. The relationship between the phase-locking strength of residual FEP SSTA and the annual mean and annual cycle amplitude of SDM parameters is illustrated in Supplementary Fig. 8. It is essential to note that the phase-locking behavior is influenced by multiple factors, which leads to the smaller correlation between phase-locking strength and SDM parameters. However, the relationship between phase-locking strength and the annual mean of α is relatively significant (R = 0.38), indicating the dominant role of the mean SSTA damping rate in controlling the spread of FEP SSTA phase-locking in climate models.

To further examine the sensitivity of the annual mean of damping rate α to the phase-locking strength, the provided coefficients in the climate model’s SDM simulations are replaced by the SDM parameters of observations, with the exception of the annual mean of α. That is, all SDM parameters, including the annual cycle of α, are set to observed values except for the annual mean of α. Figure 7d shows the dependence of the phase-locking strength on the annual mean of α in this sensitivity experiment. The increase in the preference strength of residual FEP SSTA phase-locking is a consequence of the increased mean value of α (Fig. 7d), which is consistent with previous results. The distribution of the SDM-simulated phase-locking strength (boxplot in Fig. 7d) is similar to that of the climate models (first boxplot in Fig. 7a), with most values between 0.2 and 0.3, confirming that the weak-biased strength of total and residual FEP SSTA phase-locking in climate models is attributed to an overestimation of the SSTA damping feedback (i.e., overly negative SSTA damping rate α).

Discussion

In this study, the mechanism of boreal summer SSTA phase-locking in the far eastern equatorial Pacific (FEP) and its biases in the climate models are investigated based on the linear regression analysis and linear stochastic-dynamical model (SDM). In observations, the preferred peak of SSTA events in the FEP and South American coastal regions tends to occur in the boreal summer from April to July, with fewer winter peak events. By separating the direct ENSO forcing (Niño3.4 index) from the FEP SSTA, we found that the summer peak preference is attributed by the residual SSTA component, while the ENSO forcing provides only a small probability of winter peak. Further analysis for examining the individual contribution according to the SDM simulations shows that the summer phase-locking of FEP SSTA is mainly controlled by the seasonal modulation of the SSTA damping rate α, with the strong annual cycle turning from positive to negative in boreal summer. In contrast, the contribution of ENSO-related SSTA and thermocline processes are very small.

In the climate models, the dynamics of FEP SSTA phase-locking are similar to the observations, which are dominated by the residual SSTA damping process. However, most climate models still fail to simulate a realistic summer SSTA phase-locking in the FEP region, where 40% of models do not have a peak month in the boreal summer and the preference strength of phase-locking is much smaller than observed in the climate models. Erroneous ENSO’s impact and phase-locking behaviors lead to the bias in the peak month of total FEP SSTA phase-locking in the climate models. The weak-bias in the phase-locking strength is mainly due to the smaller contribution of residual SSTA damping feedback, which is caused by the too negative annual mean of damping rate α.

In our results, we also found that the lead-lag relationship between the central-eastern Pacific SSTA and far eastern Pacific SSTA is more pronounced after decomposing the SSTA into ENSO forcing and residual SSTA components. The residual FEP SSTA led the Niño3.4 SSTA by 6 months before the mid-1970s, however, after the mid-1970s the Niño3.4 SSTA led the residual FEP SSTA instead (Fig. 2a). Although this decadal change was extensively studied by early researchers and attributed to changes in background status from the mid- to late-1970s29, the relative contributions of ENSO forcing and residual FEP SSTA in different decades remain unclear and require further study.

We provide a perspective for understanding the characteristics of FEP SSTA phase-locking and reveal the biases in current climate models. This study also suggests that when evaluating the seasonality of SST variability in the Pacific Ocean in climate models, it is important to study not only the winter SSTA phase-locking in the Niño3.4 region, but also to consider the summer SSTA phase-locking performance in the far eastern Pacific and coastal regions.

Methods

Observations and CMIP models

The observational monthly subsurface ocean temperature is obtained from four global ocean reanalysis products: Ocean Reanalysis System 5 (ORAS5) from 1958 to 201830; Simple Ocean Data Assimilation Version 2.2.4 (SODA224) from 1950 to 201031 and Version 3.3.1 (SODA331) from 1980 to 201932; and Global Ocean Data Assimilation System (GODAS) from 1980 to 201933. The observational results presented in this study are based on the ensemble average of the analysis of four datasets.

Monthly output from the historical simulations of Coupled Model Intercomparison Project Phase 5 (CMIP5)34 and Phase 6 (CMIP6)35 are used to evaluate the performance of the SSTA phase-locking in climate models. For the 45 CMIP5 and 40 CMIP6 models (Supplementary Table 2), only the first ensemble member of each model (r1i1p1 for CMIP5 and r1i1p1f1 for CMIP6) is analyzed and the period 1850–2005 is chosen to match the period of each model.

The interannual anomalies are calculated by subtracting the climatology from 1850–2005 for CMIP models and the dataset period for observations. The linear detrending is performed for calculating the anomalies.

Definition and metrics of SSTA phase-locking

An event of SSTA is defined as a 3-month running averaged time series of SSTA exceeding ±1.0 standard deviation. The peak of a positive (negative) SSTA event is selected as the maximum (minimum) peak of the 3-month running averaged SSTA time series within a 10-month time window to avoid double peaks in a single event.

The peak phase histogram is a probability histogram normalized to the sum of the values of 1 for all 12 calendar months. To quantitatively evaluate the characteristics of SSTA phase-locking from phase histogram, two metrics of phase-locking proposed by Chen and Jin9 are used here. The preferred peak month (\({\varphi }_{p}\)) is defined as the calendar month with the largest value in the 3-month averaged SSTA peak phase histogram. The metric for preference strength (\({\varphi }_{s}\)) of phase-locking is estimated as follows:

$${\varphi }_{s}=\frac{4}{3}\left({\varphi }_{\max }-\frac{1}{4}\right)$$
(4)

where \({\varphi }_{\max }\) is the sum of the 3-month values of the histogram centered on the peak month \({\varphi }_{p}\). \({\varphi }_{s}=1\) means complete phase-locking within 3 months and, conversely, \({\varphi }_{s}=0\) means there is no preferred peak month within 12 calendar months.

Linear stochastic-dynamical model

To examine the dynamics of the residual SSTA process after removing the ENSO forcing, we adopt a linear stochastic-dynamical model (SDM), which is driven by stochastic forcing. Our SDM considers seasonally varying damping process and ENSO-related processes, which is an extension of the original stochastic climate model36,37 that can be used to describe slowly varying climate variability. For the residual SSTA, the SDM can be written in the following form:

$$\frac{d{T}_{R}}{{dt}}=\alpha \left(t\right){T}_{R}\left(t\right)+\beta \left(t\right){T}_{N34}\left(t\right)+f\left(t\right){h}_{w}\left(t\right)+\sigma \left\langle \xi \left(t\right)\right\rangle$$
(5)
$$\frac{d\xi }{{dt}}=-m\xi +w(t)$$
(6)

where \({T}_{R}\), \({T}_{N34}\), and \({h}_{w}\) are the monthly area-averaged residual SSTA, Niño3.4 SSTA (5°S–5°N, 170°–120°W), and western equatorial Pacific thermocline anomalies (depth of the 20 °C isotherm over 5°S–5°N, 120°E–160°E), respectively. The first term on the right side of Eq. (5) reflects the seasonal damping by the residual SSTA with a damping rate α. The second and third terms on the right side of Eq. (5) represent the influence of ENSO-related SSTA and thermocline anomaly with coefficients β and f, respectively. The last term in Eq. (5), stochastic forcing, has an amplitude of σ and a normalized Gaussian distributed red noise \(\left\langle \xi \right\rangle\) (the angle brackets indicate the normalization). For the red noise Eq. (6), w is the white noise term and \(1/m\) denotes the decay timescale. More information about the derivation and applicability of SDM to the residual FEP SSTA can be found in Supplementary Note 1.

The seasonally varied SDM parameters for the observations and simulations are estimated by using the seasonal multivariate linear regression (Supplementary Method 1). The noise amplitude σ and decay rate m are determined as the standard deviation and exponential decay rate of the residual of SDM. The seasonality of σ is not considered here due to the insignificant effect on the seasonal modulation of SSTA. The SDM simulations are driven by the artificial stochastic forcing with applied interpolated pentad (5 days) ENSO-related SSTA (\({T}_{N34}\)) and thermocline anomalies (\({h}_{w}\)) from observations and CMIP models. All SDM simulations have 200 ensemble members and are integrated using a modified Euler’s method with a pentad time step for the corresponding period of each dataset. The monthly average dataset simulated by SDM are used to analyze in this study.