A study of the effects of westerly wind bursts on ENSO based on CESM

Numerous works have indicated that westerly wind bursts (WWBs) have a significant contribution to the development of El Niño events. However, the simulation of WWBs commonly suffers from large biases in the current generation of coupled general circulation models (CGCMs), limiting our ability to predict El Niño events. In this study, we introduce a WWBs parameterization scheme into the global coupled Community Earth System Model (CESM) to improve the representation of WWBs and to study the impacts of WWBs on El Niño-Southern Oscillation (ENSO) characteristics. It is found that CESM with the WWBs parameterization scheme can generate more realistic characteristics of WWBs, in particular their location and seasonal variation of occurrence. With the parameterized WWBs, the skewness of the Niño 3 index is increased, in better agreement with observation. Eastern Pacific El Niño and central Pacific El Niño events could be successfully reproduced in the model run with WWBs parameterization. Further diagnoses show that the enhanced horizontal advection in the central Pacific and vertical advection in the eastern Pacific, both of which are triggered by WWBs, are crucial factors responsible for the improvements in ENSO simulation. Clearly, WWBs have important effects on ENSO asymmetry and ENSO diversity.


Introduction
Westerly wind bursts (WWBs) are high frequency wind anomalies that occur over the tropical western and central Pacific Ocean. They typically have a duration of 6 days, zonal scale of 20° longitude and amplitude of 7 m/s (Harrison and Vecchi 1997;Lengaigne et al. 2003;Seiki and Takayabu 2007). The occurrence of WWBs can induce surface eastward current anomalies, trigger downwelling tropical Kelvin waves, and induce a warming effect in the equatorial central and eastern Pacific Ocean (McPhaden et al. 1988;Lengaigne et al. 2003). Observations and modeling show that WWB activities strongly affect the diversity, asymmetry, and evolution of El Niño-Southern Oscillation (ENSO) events (McPhaden 1999;Lian et al. 2014;Fedorov et al. 2014;Hu et al. 2014;. It has also been found that WWBs play a key role in limiting the forecasting accuracy for El Niño events (McPhaden et al. 2015). The addition of observed WWBs in numerical models can indeed produce better forecasts (Perigaud and Cassou 2000;Vitart et al. 2003).
Forced with observed sea surface temperature (SST), some atmospheric general circulation models (GCMs) can reproduce the basic features of WWBs, including their occurrence and variability (Vecchi et al. 2006;Miyama and Hasegawa 2014). However, it has been a challenge for coupled models to characterize WWBs realistically. Currently, almost all coupled GCMs show large biases in the simulation of WWBs, especially with regard to occurrence frequency and amplitude Lian et al. 2018). For example, Community Climate System Model version 4 (CCSM4) severely underestimates both the number of WWBs and their strength in boreal spring over the tropical Pacific region (Lian et al. 2018).
Many efforts have been made to improve the representation of WWBs in climate models for better understanding, modeling, and predicting ENSO. This effort is typically implemented by parameterizing WWBs in climate models. Although some climate models see WWBs as stochastic processes that do not depend on the ENSO cycle (Kessler et al. 1995;Moore and Kleeman 1999), many studies based on observations have found that the occurrence of WWBs is highly related to the ocean state (Harrison and Vecchi 1997;Vecchi and Harrison 2000;Batstone and Hendon 2005;Eisenman et al. 2005;Tziperman and Yu 2007). Indeed, WWBs preferentially occur prior to and during the development of El Niño events and tend to migrate eastward with the expansion of the Western Pacific warm pool (McPhaden 1999;Yu et al. 2003;Huang et al. 2019). A widely accepted hypothesis is that the occurrence of WWBs is semi-stochastic, that is, WWBs should be partially modulated by the background SST and partially determined by stochastic processes (Gebbie et al. 2007). The WWBs parameterization schemes currently used in climate models fall within this framework (Gebbie et al. 2007;Gebbie and Tziperman 2008;Lian et al. 2014;Thual et al. 2016). For example, in Gebbie et al. (2007), the trigger of an individual WWBs event is stochastic, but the probability of occurrence depends on the extent of the warm pool.
Since the 1980s, various kinds of numerical models with different degrees of complexity, from simple models (Hirst 1986), to intermediate coupled models (Zebiak and Cane 1987;Chen et al. 2004), to hybrid coupled models (Barnett et al. 1993;Tang and Hsieh 2002), and to fully coupled ocean-atmosphere model (Barnston et al. 2012;Luo and Lee 2016), have been employed for ENSO simulation and prediction. A literature review reveals that the ENSO models employed with WWBs parameterization schemes are basically noise-free models. Some of these are simple coupled models or intermediate coupled models (Perigaud and Cassou 2000;Lian et al. 2014).  and  used fully coupled models to explore the roles of WWBs in ENSO but the models they used lack atmospheric variability over the WWBs domain, allowing to directly add WWBs into models by parameterization schemes. It was reported that the bias toward cold events is reduced by SST-dependent WWBs parameterization, but the skewness of ENSO is still negative and the models actually produce more extreme cold events.
There has been little work to introduce WWBs into a fully coupled model that generates its own WWBs for two key reasons. The first is that such CGCMs produce spurious WWBs (Lian et al. 2018). It has been a major challenge to filter misrepresented WWBs before introducing a parameterization scheme without hurting the model's dynamical and thermos-dynamical balance. The second is that all CGCMs have large model biases in the simulation of tropical climatology, such as the extent of the warm pool and cold tongue in the Pacific Ocean, which are critical in all WWB parameterization schemes. With the rapid increase in available computing resources, fully coupled models have become the main ENSO operational prediction systems. Therefore, improving WWBs representation in global coupled models is a crucial issue for improving the simulation and forecast capabilities of ENSO.
In this study, we use the global coupled Community Earth System Model (CESM), released by the National Center for Atmospheric Research (NCAR), as an example for parameterizing WWBs in a fully coupled model that has shown spurious WWBs. We emphasize the improved modeling of ENSO features, including its asymmetry and diversity, by introducing WWBs scheme in CESM. This introduction is a necessary step for further employing ENSO prediction.
The remainder of the present paper is organized as follows. Section 2 describes the model, the WWBs parameterization scheme, and the method of filtering CESM's built-in WWBs. Section 3 presents the comparison of observed WWBs, built-in WWBs in CESM, and parameterized WWBs. Section 4 investigates the ENSO simulation in CESM in response to built-in WWBs and parameterized WWBs. Section 5 provides a conclusion of these results.

Model and method
CESM is a coupled climate model composed of seven modules: atmosphere, land, river runoff, ocean, sea ice, land ice, and ocean wave, as well as a coupler (Vertenstein and Craig 2012). It has been widely used in climate study, including ENSO research (Deser et al. 2012). In this study, we use the CESM version 1.2.2. The atmospheric component used here is CAM4, and the oceanic component is POP2. The horizontal resolution is set as f09_g16, which refers to atmospheric gridding at 0.9° latitude × 1.25° longitude, and a Greenland pole 1° ocean grid. We conduct two sets of experiments in this study. The first one is the control run without any change in CESM. The second one is the forced run in which we introduce a WWBs parameterization scheme into CESM. In both experiments, the CESM is run 70 years under present-day forcing, initialized from a state at which the model reaches the climatological equilibrium under present-day climate forcing.
The WWBs parameterization scheme is based on Gebbie et al. (2007). The probability of occurrence is defined as: where x pool is the warm pool edge that is defined as the longitude of 28.5 °C isotherm, and P 0 is a constant value of WWBs occurrence probability. A WWBs event is triggered only if a random number between 0 and 1 is less than P(t). In both temporal and spatial variation, WWBs are assumed to have a Gaussian shape. The zonal wind stress can be defined as: (1) where A is the magnitude, T is the duration, L x and L y are the zonal and meridional width, T 0 is the time of maximum wind, and x 0 and y 0 are the center longitude and latitude. The definition of WWBs event used in this study is the same as that in Lian et al. (2018), namely, the surface westerly wind anomalies averaged over the latitudinal band between 5°N and 5°S should be above 5 m/s, sustain for 2 or more days, and show a zonal extension of at least 10° in longitude. Based on this criteria, we identify each WWBs event and retrieve its corresponding feature parameters defined in Eqs. (1) and (2) during the period from 1948 to 2016, using the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/ NCAR) reanalysis daily surface zonal wind dataset and the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) dataset (Rayner et al. 2003). The horizontal resolutions of NCEP/NCAR reanalysis and HadISST dataset are 2.5° longitude × 2.5° latitude and 1° longitude × 1° latitude, respectively. Averaging each parameter over all WWBs events, we obtain climatological values for these parameters, which are used in Eqs. (1) and (2) except the center longitude of WWBs of x 0 , namely, P 0 = 0.06/day, A = 0.15 N/ m 2 , T = 6 days, L x = 22°, L y = 6°, T 0 = t init + 12.5 days ( t init is the triggering time), and y 0 = 0. These values are very similar to those used in previous research (Harrison and Vecchi 1997;Eisenman et al. 2005;Gebbie et al. 2007). For x 0 , its observed climatological value is x pool − 25°. However, we found that CESM has a systematic bias in the position of the warm pool, which dominates x 0 . The warm pool edge is 15° more eastward in the model than in observation (see Fig. 2). Thus, we modify x 0 as x 0 = x pool − 40°. This modification is also verified by sensitivity experiments. When x 0 is set as x 0 = x pool − 25°, WWBs occur excessively over the central and eastern Pacific but too little over the western Pacific.
As WWBs are misrepresented in CESM, we have to filter spurious WWBs before introducing the WWBs parameterization scheme into the model. This is implemented by reconstructing the wind stress field. As in Gebbie et al. (2007) and Tziperman (2008, 2009), when WWBs is triggered, the total zonal wind stress can be decomposed as: where ̄ * is the seasonal climatology of the non-WWB field, and � * is the non-WWB field anomaly. One widely used method to reconstruct � * is to seek its statistical relation to SST anomaly (SSTA) by singular value decomposition (SVD ;Neelin 1990;Chang et al. 2001;Zhang et al. 2005;Zheng and Zhu 2010). Here we use the 30-year daily (3) x = WWB +̄ * + � * , mean surface zonal wind stress and daily SST output from the CESM control experiment to construct � * by the SVD approach. That is, the surface zonal wind stress anomalies and SSTA are first calculated relative to seasonal climatology. To remove the WWBs contained in the wind stress anomalies, we identify each WWB during the period of the control run based on the above definition. Wind stress anomalies greater than 0.06 N/m 2 , which is approximately equivalent to the 5 m/s of the WWB definition, are subtracted from the zonal wind stress anomalies during the identified WWB period. Then, the monthly mean is applied to the daily wind stress anomalies and daily SSTA. In short, the procedure of adding parameterized WWBs as below: (1) at each time step, the condition of the occurrence of WWBs is examined. If it is satisfied, the zonal wind stress in the ocean model is replaced by Eq. 3, which uses oceanic model SST as inputs and SVD technique, and then the CESM integrates forward. If the condition is not triggered, the original CESM integrate forward. (2) Repeat (1) until the completion of the experiment period. Since we only focus on the WWBs occurring in the tropical Pacific region, the reconstructed wind stress anomalies are restricted to between 20°N and 20°S and between 120°E and 70°W. Thus, the filtering and parameterized WWBs are both done in the framework of coupling run. It should be noted that only deterministic part of WWBs is filtered, if the model SST variation is regarded to be completely deterministic. However, one may argue that the modeled SST, which determines WWBs here, could be random in nature due to the nonlinearity or randomness of coupled GCMs. As such, the filter also removes some stochastic parts of WWBs.
In order to verify the simulations from the two experiments against observations, monthly mean potential temperature and salinity, and zonal wind stress from Simple Ocean Data Assimilation (SODA, version 3.3.1) dataset (Carton et al. 2018) are used here. The resolution is 1/2° × 1/2°× 50 lev. Potential temperature and salinity fields are used to calculate the depth of 20 °C isotherm (D20) that represents the thermocline depth. Figure 1a shows the climatologically longitudinal distributions of WWBs occurrence over the tropical Pacific region for the observational data, the CESM built-in experiment (also called control run sometime hereafter) and the parameterized experiments. The definition of WWBs is the same for the observational data and the CESM control experiment. As can be seen, WWBs occurrence in observation approximately follows a normal distribution, centered in the central equatorial Pacific around 160° E-180° E. This distribution is well simulated in both experiments. The  To further explore the spatial and seasonal variations in WWBs occurrence, we calculate the occurrence rate in the Western Pacific (WP, 120° E-160° E) and Central Pacific (CP, 160° E-160°W) as a function of calendar month. The WP and CP are classified based on the criteria used in Lian et al. (2018), which span most of the area of WWBs occurrence. Compared with observation, CESM shows more WWBs occurring in these two regions, especially in the WP region. The occurrence rate of WWBs in the two regions is more realistic under the WWBs parameterization scheme than under the CESM control experiment. Over the WP and CP regions, there are 1.87 and 2.36 WWBs occurring per year in observation. In the built-in control experiment, these numbers are overestimated: 2.88 and 2.80 per year, respectively. While in the forced run with the WWBs parameterization scheme applied, the average numbers of WWBs are 1.59 and 2.75 per year, respectively, which are closer to the observed averages.

WWBs comparison
The monthly distribution of WWBs occurrence is shown in Fig. 1b-d. As can be seen, over the entire equatorial Pacific, observed WWBs peak in boreal winter and reach their minimum occurrence in boreal summer (Harrison and Vecchi 1997). This is considerably different in the control run, in which the occurrence of WWBs is underestimated in boreal winter and overestimated in spring and autumn, resulting in a spurious seasonal variation. Using the parameterization scheme, the occurrence of WWBs is more realistically simulated in boreal autumn and spring. However, the occurrence in boreal summer is not significantly improved by the parameterization scheme, which results in a weak seasonal variation in the parameterized WWBs.
In the WP region, the control run shows excessive WWBs in boreal winter and spring. The parameterized WWBs in this region have a seasonal variation of occurrence similar to that in observation. In the CP region, the seasonal variations of observed WWBs, built-in WWBs, and parameterized WWBs are similar to those in the entire Pacific region (Fig. 1b). In short, the parameterization scheme significantly improves the modeled occurrence of WWBs in the WP and CP regions, including their longitudinal distribution and their monthly frequency. The improvement is more significant during the times when WWBs preferentially occur: boreal autumn, winter, and spring.

Effects of WWBs on ENSO
The oceanic response to parameterized WWBs is the focus of this study. Before discussing the effects of WWBs on ENSO characteristics, we investigate how WWBs affect the tropical mean state and seasonal cycle. As shown in Fig. 2, the simulated SST in the CESM is warmed than observed record throughout the tropical Pacific basin. That is why we modify the equation for the location of WWBs as mentioned in Sect. 2. The pattern and amplitude of the simulated tropical surface zonal wind stress are similar to those in observation. However, the wind stress is slightly underestimated by about 0.1 dyne/cm 2 over the western Pacific, and overestimated by about 0.15 dyne/cm 2 over the central Pacific around 10° N. The differences of SST and zonal wind stress between the control run and forced run are quite small. With the parameterized WWBs, the zonal wind stress is increased by about 0.15 dyne/cm 2 over the western and central Pacific Ocean. The SST is slightly cooled over the western Pacific but warmed by about 0.2 °C over the central eastern Pacific.
Like the annual mean state, the tropical seasonal cycle is also weakly affected by parameterized WWBs. With   Fig. 3 a Seasonal cycle of Niño 3.4 SST (unit: °C) and climatological annual mean thermocline depth (D20, depth of the 20 °C isotherm) from observation (black), the control run (red) and the forced run (blue) the parameterization scheme, Niño 3.4 SST is slightly increased by about 0.2 °C (Fig. 3a), with more warming during boreal autumn and winter. However, the CESM has a semiannual seasonal cycle in both control and forced run, which is different form the observation. As shown in Fig. 3b, the climatological D20 is shallower over the western equatorial Pacific but deeper over the eastern Pacific in both control and forced run than in observation, especially for the forced run. The weaker thermocline slope in the forced run can be a reminiscent of Bjeckness feedback, suggesting that parameterized WWBs are stronger than the spurious counterpart in the control run over the equatorial Pacific.
In order to investigate the impacts of WWBs on the statistical features of ENSO events in detail, this section is divided into three parts. In Sect. 4.1, we examine the effects of WWBs on ENSO evolution, especially on the evolution of El Niño events. The effects of parameterized WWBs on ENSO asymmetry are explored through the analysis of the Niño 3 index in Sect. 4.2. In the last part, we discuss the ENSO diversity in response to WWBs and the possible mechanisms behind it. Figure 4 shows the evolutionary trajectories of the surface zonal current anomalies, thermocline depth anomalies, and equatorial Pacific SST anomalies averaged from 5° N to 5° S in the forced run. As can been seen, El Niño events, especially extreme El Niño events, are always companied by strong westerly zonal wind stress, as noted in some previous works (Kerr 1999;Vecchi and Harrison 2000;McPhaden 2004). During the development of El Niño events, the westerly wind stress tends to migrate eastward from the western to the central Pacific, leading to a strong eastward surface zonal current anomaly and potentially resulting in warming in the central and eastern Pacific Ocean.

Effects of WWBs on ENSO evolution
We conduct the composite analysis of all El Niño and La Niña events based on observation and the two model experiments. Here, El Niño and La Niña events are defined Black lines denote regions with zonal wind stress greater than 0 dyne/ cm 2 . Red bars in c indicate El Niño events by the SSTA greater than 0.5 °C and less than − 0.5 °C, respectively, as averaged over the Niño 3.4 region (5°N-5°S, 170°W-120°W). Figure 5 shows the composite SSTA in observation and in the two experiments, as well their differences. We can find that the amplitude of SSTA over the central and eastern Pacific is underestimated for El Niño events but overestimated for La Nina events in the control run. In the forced run, the amplitude of El Niño strength over the Niño 3.4 and Niño 3 regions is markedly increased by about 0.4 °C compared with that in the control run, whereas the amplitude of La Niña event is only slightly changed in these regions. The results are similar when the composite analysis is performed using the Niño 3 region to define ENSO events (not shown). The stronger El Niño amplitude in the eastern Pacific suggests that parameterized WWBs could improve the simulation of EP El Niño, in particular, the extreme El Niño events in the CESM.
Further investigation emphasizes the effects of WWBs on El Niño phase. Shown in Fig. 6 is the evolution of surface zonal wind stress anomalies, SST anomalies, and D20 anomalies during an El Niño cycle. This figure, obtained using the composite analysis based on observation and the 70-year coupling run, tracks these anomalies through the development, mature, and decay phases. As can be seen, during the development phase of El Niño, the westerly wind anomalies prevailing over the western and central Pacific are underestimated in the control run compared with observed counterparts. In the forced run with parameterized WWBs, the amplitude of the westerly wind anomalies is much enhanced.
The westerly wind stress anomalies drive warm waters eastward, accompanied by the expansion of the warm pool. The strongest wind stress anomalies occur during the period from boreal summer to winter, resulting in the maximum warming of eastern Pacific Ocean sea surface water in winter. As expected, accompanied with stronger wind stress anomalies in the forced run, a deeper thermocline exists in the eastern Pacific region which is important to the development of El Niño, especially for the eastern Pacific El Niño . The evolution and the amplitude of the thermocline depth anomalies by parameterized WWBs are much closer to those in observation. Consequently, the amplitude of warming anomalies is strengthened in the forced run and in better agreement with observed result. In the composite, the eastern Pacific SST anomalies peak in boreal winter, which is consistent with the phase locking of ENSO as seen in observational data, although the maximum warming still does not extend to the coast of South America.

Effects of WWBs on ENSO asymmetry
The histograms of the Niño 3 index (3 month running mean of SSTA in the Niño 3 region) are shown in Fig. 7. Observational data shows that the SSTA warming anomaly is stronger than its cooling counterpart, with a skewness value of 0.89, indicating a strongly asymmetric distribution of the Niño 3 index. The control run fails to represent this asymmetry well. As can be seen in Fig. 7b, a normal distribution approximately holds for the Niño 3 index, with a skewness value of 0.11. In other words, CESM is unable to well characterize the ENSO asymmetry. By contrast, the forced run has a larger skewness of 0.44 and the Niño 3 index has a more asymmetric distribution as shown in Fig. 7c. These are significantly different from the results of CCSM3, an ancestor of the model used here, as reported in . They found that in the forced run with SST-dependent WWBs, although the bias toward cold events is reduced, the skewness of ENSO is still negative and CCSM3 actually produce more extreme cold events (sometimes even exceed − 3.0 °C) than warm events.
The improvement of ENSO asymmetry in our forced experiment is mostly because WWBs are typically triggered during the growth phase of warm events. The enhanced westerly wind anomalies can generate eastward downwelling Kelvin waves and anomalous eastward currents, thereby  (Bjerknes 1969). Namely, the intensity of ENSO warming phase is significantly amplified by WWBs activities. This explains why extreme El Niño events can often occur with the occurrence of strong WWBs. More details about the physical processes will be discussed in next subsection. By improving the representation of WWBs, the forced run can successfully reproduce the features of ENSO asymmetry.

Effects of WWBs on ENSO diversity
Many studies have shown that El Niño warming centers are typically distributed in the central Pacific and eastern Pacific Ocean, with corresponding events defined as central Pacific El Niño (CPEN) and eastern Pacific El Niño (EPEN) events (Ashok et al. 2007;Kao and Yu 2009;Capotondi et al. 2015). The mechanisms of the two kinds of El Niño events have been a hot issue in ENSO study, producing several classic hypotheses. One of these hypotheses, proposed by Lian et al. (2014) and , is that WWBs perturbations play an important role in ENSO diversity and trigger CPEN events through zonal heat advection. However, they only used the Zebiak-Cane model (Zebiak and Cane 1987), a highly simplified coupled model, to verify this hypothesis.  employed CCSM3 and CCSM4 with a coarse resolution to investigate the effect of WWBs on ENSO diversity. They found that SST-dependent WWBs are much more effective on EP El Niño than on CP El Niño. It is therefore of interest to further explore the impact of WWBs on ENSO diversity using an advanced coupled model with a finer resolution. Figure 8 shows two leading rotated empirical orthogonal function (REOF) modes of tropical SSTA from the control run, the forced run, and observation, respectively. Compared with traditional EOF, REOF can better characterize the spatial distribution of signal variance, and thus, it is widely applied in pattern recognition and classification (Kaiser 1958;Lian and Chen 2012). For the observational data, the first two modes clearly represent the EP and CP El Niño, respectively. In the control run, however, CESM can only produce the EP El Niño as shown in its first mode, but fails to reproduce the CP El Niño in the second or higher modes. When the parameterized WWBs are added in the forced run, both the EP and CP El Niño types can be well simulated in the leading two modes, which are similar to those in observation. Thus, the hypothesis of ENSO diversity mechanism, proposed by , is well supported not only by simple model but also by the fully coupled model. This agreement suggests the generality and robustness of this mechanistic hypothesis.
We can gain deeper insight by further exploring the physical processes underlying the role of WWBs in ENSO diversity. To achieve this, we want to identify the main contributions to the development of the two El Niño types. For this purpose, we use a heat budget analysis of the temperature anomaly of the ocean mixed layer (from surface to 50 m depth), as in Kug et al. (2009). The change in temperature anomaly with time can be decomposed as: where u, v, w, and T denote, respectively, the zonal current, meridional current, vertical velocity, and temperature averaged from 5° N to 5° S in the mixed layer. x, y, and z indicate the zonal, meridional, and vertical direction. Overbars and primes signify the monthly mean climatology and anomalies. Q ′ net denotes net surface flux anomaly which is the sum of latent heat, sensible heat, short wave radiation and long wave radiation fluxes anomalies. ρ, C p and H denote sea water density, specific heat and mixed layer depth, respectively. The terms −u � xT and −v � yT are anomalous zonal and meridional advection terms. The −ū x T � and −v y T � (4) terms denote mean zonal and meridional advection terms. The −w � zT term denotes Ekman feedback, and the −w z T � term denotes thermocline feedback. The − u � x T � , − v � y T � , and − w � z T � terms are nonlinear advection terms. R is the residual term.
As in Lian et al. (2014), EP El Niño events are defined as those whose averaged SSTA between 5° N and 5° S and between 110°W and 90°W is larger than 0.6 °C and lasts at least 6 months. Whereas, CP events are defined as those whose averaged SSTA between 5° N and 5° S and between 175°W and 155°W is continuously above 0.5 °C for at least 3 months. Over the 70-year period of the model experiment, there are 12 CPEN events and 14 EPEN events in the forced run.
The temporal evolution of each heat budget term during the development and decay periods (9 months before and after the peak SSTA) of each El Niño type is shown in Fig. 9. The temperature anomaly tendency reaches the maximum around 4 and 5 months before the peak warming for CPEN and EPEN, respectively. During the development phase of CPEN, the meridional and zonal advection terms are the dominant factors whereas the vertical advection seems negligible. A further diagnosis, as shown in Fig. 9c, reveals that zonal advection is dominated by the anomalous zonal advection term ( − u � xT ), which probably results from a stronger zonal flow driven by WWBs. Meridional advection is highly related to the mean meridional advection term ( −v y T � ), suggesting that poleward mean flow expands the positive SST anomalies around the equator. Vertical advection appears to remain very minor during the entire CPEN event cycle. These results confirm the importance of horizontal advection for the development of CPEN events as suggested in some previous works (Kug et al. 2009;Lian et al. 2014).
In contrast to CPEN, the dominant contributor to the development of EPEN is vertical advection, followed by meridional advection. Zonal advection shows only a very small contribution to EPEN. Among the vertical advection terms, the thermocline feedback term ( −w z T � ) gives a significant contribution to the positive SSTA tendency, as shown in Fig. 9d. Comparable significant contributions during the development phase are also made by the Ekman feedback term ( − w � zT ) and mean meridional advection term ( −v y T � ), indicating their importance in the growth of EPEN (see Fig. 9d). The net heat flux term has a negative effect on the development of El Niño event, especially on EP El Niño type. Further analysis found that the short wave radiation flux dominates the net heat flux term during the development of CPEN, while the latent heat flux dominates the net heat flux term during the development of EPEN (figure not shown here).
The difference in the contributable items between CPEN and EPEN is a sign of different sets of physical processes and climatological features in the two event types. For example, the thermocline is climatologically deep in the western and central Pacific but shallow in the eastern Pacific. As a result, the thermocline variation is weaker in the central Pacific but stronger in the eastern Pacific, resulting in a In order to show why CESM cannot simulate CPEN in the control run. Figure 10 explores the differences between the control run and forced run during the cycle of an El Niño event, with particular regard to the nine advection terms. In the central Pacific region, the parameterized WWBs in the forced run can produce eastward zonal current anomalies in the central Pacific, as shown in Fig. 4a, which increases the anomalous zonal advection − u � xT , as shown in Fig. 10a. The enhanced mean zonal and meridional advection terms ( −ū x T � and −v y T � in Fig. 10b, f) in the central Pacific might result from the forced run's stronger SST anomalies (shown in Fig. 6j). The vertical advection terms, either from thermocline feedback or Ekman feedback or both, are slightly changed in the central Pacific. While in the eastern Pacific region, because the WWBs deepen the thermocline and weaken the upwelling, the thermocline feedback and Ekman feedback terms are enhanced significantly which are typically the dominant factors to control the SSTA variability there. This is consistent with the conclusion on the importance of thermocline feedback to the development of EPEN . The difference of heat flux contribution between the control run and the forced run is negligible, compared with those of advection terms (figure not shown here).
These results indicate that, with parameterized WWBs, the central Pacific experiences enhanced horizontal advection (Fig. 10d, h), and the eastern Pacific experiences enhanced vertical advection (Fig. 10l); these are the main drivers of the occurrence of CPEN and extreme EPEN. In other words, lacking of sufficient horizontal advection in the central Pacific and vertical advection in the eastern Pacific in the control run is probably a main factor responsible for the absences of CPEN and the underestimate of the magnitude of EPEN.

Conclusion and discussion
Over the past decades, ENSO study has been an intensive research field in oceanic and atmospheric science research. With the great efforts of the researchers, significant progress of ENSO study has been made (Jin et al. 2008;Barnston et al. 2012). However, ENSO prediction still suffers uncertainties. For example, the latest El Niño Fig. 9 Oceanic mixed layer heat budget during CPEN event and EPEN event based on the forced run. Black, green, orange, red, blue and purple lines in the top panels denote the temperature anomaly tendency, residual contribution, heat flux term, zonal advection, meridional advection and vertical advection, respectively. Solid, dotted and dashed lines in the bottom panels indicate three different contributions to each advection term as shown in Eq. 4. Unit: °C/month event from 2014 to 2016 poses a severe challenge to the classic ENSO theory-based forecasting models . Almost all models failed and missed the strongest warming event when predictions were initialized in early 2015. A main reason is probably that almost all of these models failed to realistically capture the WWBs activities those occurred in spring 2015 (Chen et al. 2017).
In this study, we have attempted to improve the WWBs representation in a fully coupled model, which serves as a first step toward the final goal of improving ENSO prediction skill. For this purpose, we introduced a WWBs parameterization scheme into CESM based on the framework of Gebbie et al. (2007). We comprehensively evaluated the improvement of both the WWBs simulation and Eq. 4. The right row denotes the differences between the control and forced run in total zonal, meridional and vertical advection. Unit: °C/ month the corresponding features of ENSO by the parameterization scheme. This evaluation was implemented using two parallel experiments, one with CESM built-in WWBs (control run) and the other with parameterized WWBs (forced run). It was found that the CESM built-in WWBs show significant differences from observed counterparts, especially in the WWBs occurrence. Before adding the parameterized WWBs into CESM, we first removed the built-in WWBs produced in CESM by reconstructing the non-WWBs wind stress field using Singular value decomposition (SVD) technique.
With the observation-derived parameters, it was found that the WWBs parameterization scheme can produce more realistic characteristics of observed WWBs, particularly their location and seasonal variation of occurrence. The effects of WWBs on ENSO in these experiments are consistent with some previous works Fedorov et al. 2014;Hu et al. 2014;. With the parameterized WWBs, surface wind stress is intensified, enhancing the surface current and eastward downwelling Kelvin waves, thereby resulting in deeper thermocline and greater SST anomalies. Meanwhile, the magnitude of El Niño is greater than that in the control run, while the amplitude of La Niña is only slightly affected. Consequently, the skewness of the Niño 3 index in the forced run is increased and in better agreement with observation. Namely, WWBs could play an important role in generating the asymmetry between El Niño and La Niña events.
Analysis shows that CESM with built-in WWBs in the control run fails to reproduce CPEN type. By contrast, when the WWBs parameterization scheme is introduced into the model, the simulation of CP El Niño events could be improved, especially at the location of maximum warming. A further investigation on the heat budget shows that the cycle of CPEN and EPEN events is dominated by different feedback terms. The anomalous zonal advection and mean meridional advection are two major factors responsible for the development of CPEN, whereas the vertical (both thermocline and Ekman feedbacks) and mean meridional advection dominate the growth of EPEN. With the parameterized WWBs, zonal and meridional advection in the central Pacific and vertical advection in the eastern Pacific are significantly strengthened in CESM, and result in increased SST anomalies in the central and eastern Pacific regions. The intensified advection terms driven by WWBs could be responsible for the occurrence of CPEN and extreme El Niño events.
In short, CESM with parameterized WWBs has better simulation skills in terms of ENSO asymmetry and diversity. One may argue that, the forced run results in the mean states of SST and surface zonal wind further from observation over the central and eastern Pacific, as shown in Fig. 2. However, it producers a better simulation of ENSO. This better simulation is probably from a random correction to some spurious features of ENSO. It is most due to the fact that the WWBs parameterization can better capture Bjerknes feedback in model. Previous works have noted that most climate models underestimate the wind-SST feedbacks, especially zonal advective feedback and thermocline feedback, which are thought to play an important role in ENSO physics (Kim et al. 2014;Ferrett and Collins 2019). The feedback terms can be defined as a combination of mean states and the responses of atmosphere to SST change and ocean to wind change. The forced run with a flatter thermocline, i.e., shallower thermocline in the west and deeper in the east, the weaker trade wind over the equator, and warmer temperature over the eastern Pacific, however, could result in a stronger response of ocean to wind stress, thereby enhancing the feedback processes as Kim et al. (2014) has reported in their work. In the forced run, although the model biases of mean states are slightly larger, the evolution of surface wind anomalies, thermocline depth anomalies and SST anomalies are better simulated as shown in Fig. 6. As a result, accompanied with stronger wind anomalies and thermocline depth anomalies, the positive advection terms during El Nino development phase are enhanced and bring better simulation of ENSO characteristics (Fig. 10). In other words, one may not conclude that a better representation of mean states in climate models must lead to a better simulation of ENSO. Besides that, the parameterized WWBs can result in a longer life span of the warm event than observation. In addition, we explored only one WWBs parameterization scheme, which assume some parameters related to the basic characteristics of WWBs, including the magnitude, duration, and center latitude, are constant. Other schemes remove the assumption, for example, the one proposed by Tziperman (2008, 2009) describes the characteristics of WWBs as linear functions of SST but with flow-dependent characteristic parameters. Thus, it seems necessary to conduct further study using different WWBs parameterization schemes, which is under our way. Nevertheless, this work investigates the effect of WWBs on ENSO in the framework of a globally fully coupled model, offering a theoretical significance and practical importance for improving ENSO predictions.