ENSO statistics, teleconnections, and atmosphere–ocean coupling in the Taiwan Earth System Model version 1

This study provides an overview of the fundamental statistics and features of the El Niño–Southern Oscillation (ENSO) in the historical simulations of the Taiwan Earth System Model version 1 (TaiESM1). Compared with observations, TaiESM1 can reproduce the fundamental features of observed ENSO signals, including seasonal phasing, thermocline coupling with winds, and atmospheric teleconnection during El Niño events. However, its ENSO response is approximately 2 times stronger than observed in the spectrum, resulting in powerful teleconnection signals. The composite of El Niño events shows a strong westerly anomaly extending fast to the eastern Pacific in the initial stage in March, April, and May, initiating a warm sea surface temperature anomaly (SSTA) there. This warm SSTA is maintained through September, October, and November (SON) and gradually diminishes after peaking in December. Analysis of wind stress–SST and heat flux–SST coupling indicates that biased positive SST–shortwave feedback significantly contributes to the strong warm anomaly over the eastern Pacific, especially in SON. Our analysis demonstrates TaiESM1's capability to simulate ENSO – a significant tropical climate variation on interannual scales with strong global impacts – and provides insights into mechanisms in TaiESM1 related to ENSO biases, laying the foundation for future model development to reduce uncertainties in TaiESM1 and climate models in general.

1 Introduction

The El Niño–Southern Oscillation (ENSO) is the primary mode of interannual and decadal climate variability in the tropics (Glantz, 2001; McPhaden et al., 2006). It also affects climate variations in subtropical and midlatitude regions across both hemispheres through teleconnection of Rossby waves (Diaz et al., 2001; Yeh et al., 2018). Therefore, its prediction is an essential part of global climate prediction on these scales (Latif et al., 1998). Additionally, ENSO is a crucial metric for climate model evaluation, especially for atmosphere–ocean coupling and associated physical feedbacks (Planton et al., 2021). Many studies have reported that Coupled Model Intercomparison Project (CMIP) models have successfully represented the basic features of observed ENSO, such as a recognizable ENSO life cycle and sea surface temperature (SST) pattern over the tropical central and eastern Pacific (Guilyardi et al., 2009, 2020; Lloyd et al., 2009, 2011; Bellenger et al., 2014). However, many model biases are also found in CMIP6 models (Beobide-Arsuaga et al., 2021; Capotondi et al., 2020a; Chen and Jin, 2021), resulting in 30 %–50 % uncertainties in future ENSO projection (Beobide-Arsuaga et al., 2021).

The paradigm based on theory and observations depicts ENSO as a closely coupled oscillation system between the ocean and the atmosphere (Jin, 1997; Latif et al., 1998; McPhaden et al., 1998, 2006; Neelin et al., 1998; Philander, 1989; Wang, 2018; Wang and Picaut, 2013). An El Niño event begins from a westerly wind initiated over the western equatorial Pacific, typically in the spring of the first year (i.e., March, April, and May in year 0; MAM⁰). The westerly wind drives more warm water towards the east and gradually warms the central Pacific. Around summer, eastward-propagating oceanic Kelvin waves are triggered over the central Pacific, reducing the upwelling at the tropical eastern Pacific, deepening the thermocline, and warming the sea surface temperature (SST). Consequently, the zonal SST gradient and the easterly wind in the tropical Pacific are reduced. Such a retreat of the easterly wind further reduces the SST gradient, causing the so-called Bjerknes feedback between the easterly wind and the SST gradient (Bjerknes, 1969; Cane, 2005). Through this feedback, the warm SST increases and reaches a maximum around the following winter (i.e., December of year 0 and January and February of year 1; DJF⁺¹). Furthermore, following the warm SST anomaly (SSTA), the center of deep convection activity shifts toward the central Pacific, increasing latent heat flux and reducing shortwave heat flux into the ocean surface through deep cloud cover. Such seasonal phase locking of an El Niño event is a crucial characteristic of the observed ENSO.

In contrast to the coupling nature of the atmosphere and ocean found in ENSO observations, the atmospheric feedbacks are found to dominate the modeled ENSO frequency and amplitude in CMIP models (Guilyardi et al., 2009; Lloyd et al., 2009, 2011; Bellenger et al., 2014; Beobide-Arsuaga et al., 2021). It has been noticed that the CMIP models tend to simulate a weaker Bjerknes feedback, namely, a weaker SST warming and westerly wind coupling. Furthermore, the heat flux–SST feedbacks are overemphasized in simulated ENSO dynamics, especially for the shortwave heat flux–SST feedback. Such overemphasized heat flux–SST feedback compensates for the weaker warming from the Bjerknes feedback, producing a seemingly realistic ENSO warming in CMIP models (Bayr et al., 2019). The biased Bjerknes and heat flux feedbacks were later found to be related to the biases of the seasonal variations of Walker circulation (Bayr et al., 2019, 2018). Such complexity resulting from intertwined atmospheric–ocean feedbacks makes it challenging for model developers to improve ENSO simulations without fully understanding how these mechanisms are represented in the coupled models.

Taiwan Earth System Model version 1 (TaiESM1; Lee et al., 2020) is the Earth system model developed at the Research Center for Environmental Changes (RCEC), Academia Sinica. It participates in the CMIP6 intercomparison activity and has been used to study major climate variabilities and regional climate features (Chen and Jin, 2021; Park et al., 2020). While its overall performance of climate mean states and major variations has been evaluated and documented in Wang et al. (2021), in this study, we conducted a more comprehensive investigation of ENSO's fundamental features and statistics in historical TaiESM1 simulations. We noticed that the ENSO amplitude increased significantly compared with the observations, with intense and prolonged warm SST from May to December. Especially over the eastern Pacific, the early onset of warming and sustained warming in September, October, and November (SON⁰) before peaks in December are the two primary prominent biases. We further analyzed the physical processes associated with ENSO's strong warming biases and found it is due to the biased positive feedback between SST and shortwave surface fluxes over the eastern equatorial Pacific. Such feedback is primarily attributed to the prevailing low clouds overlying the cold tongue region, with large cold biases in SON⁰. Our results provide the baseline for ENSO performance in TaiESM1 and suggest that the current ENSO biases in TaiESM1 are intertwined with biases of the mean state and seasonal variation of the tropical climate system. The remainder of this study includes the following. Section 2 describes the TaiESM1 and observational dataset used for model evaluation and the methodology for analyzing ENSO. Section 3 documents the basic characteristics of ENSO simulated in TaiESM1. More analysis focuses on seasonal variation of El Niño events in TaiESM1 and associated biases in Sect. 4. The study is concluded with a summary and discussion in Sect. 5.

2 Data, models, and methodology

Based on the Community Earth System Model version 1.2.2 (CESM1.2.2; Hurrell et al., 2013) developed by the National Center for Atmospheric Research and sponsored by the National Science Foundation and the Department of Energy in the United States, TaiESM1 includes several physical schemes developed in-house in RCEC. These designs include convective triggering (Y. C. Wang et al., 2015), radiation parameterization of three-dimensional topography (Lee et al., 2013), an aerosol scheme (Chen et al., 2013), and a cloud fraction scheme based on probability density function (Shiu et al., 2021). The ocean component is the same as CESM1 using the Parallel Ocean Program version 2 (POP2; Smith et al., 2010). An overall evaluation of climate variability in TaiESM1 shows that the simulated ENSO features stronger SST warming and atmospheric teleconnection compared with the base model CESM1 (Wang et al., 2021). In this study, we have conducted a more in-depth analysis of ENSO's fundamental features and statistics and identified the physical processes of ENSO biases within TaiESM1. The historical simulation of TaiESM1 from 1850 to 2014, driven by the forcing designed by CMIP6, is analyzed. The historical run is initiated from the pre-industrial control run of TaiESM1. It utilizes an atmospheric model with a horizontal resolution of 0.9^∘ latitude × 1.25^∘ longitude and 30 vertical layers. The community land model employed shares the same resolution as the atmospheric model. Additionally, the POP2 ocean model has a resolution of approximately 1.125^∘ in longitude and 0.47^∘ in latitude.

We evaluated the model's performance using the atmospheric variables from the Collaborative Reanalysis Technical Environment Multireanalysis Ensemble version 2 (MRE2; Potter et al., 2018). The MRE2 is a product of the ensemble average of seven reanalysis products, including CFSR (Saha et al., 2010), ERA-Interim (Dee et al., 2011), MERRA (Rienecker et al., 2011), MERRA-2 (Gelaro et al., 2017), JRA-25 (Onogi et al., 2007), JRA-55 (Kobayashi et al., 2015), and 20CRv2c (Compo et al., 2011). Studies have found that the ensemble average can reduce the errors of individual reanalysis for selected atmospheric variables (Potter et al., 2018). For those variables not provided in MRE2, such as cloud cover, we used the ECMWF reanalysis version 5 (ERA5), the most up-to-date reanalysis produced by ECMWF (Hersbach et al., 2020). While previous studies have identified differences in air–sea feedbacks among reanalysis datasets, the ensemble mean of multiple reanalysis datasets can be used as the best estimate by reducing random errors through averaging (Kumar and Hu, 2012). We also obtained precipitation data from the Global Precipitation Climatology Project (GPCP V2; Adler et al., 2003; Huffman et al., 2009), SST data from the Extended Reconstructed SST version 5 (ERSSTv5; Huang et al., 2017), and subsurface ocean data, such as sea surface height (SSH) and potential subsurface temperature, from the Simple Ocean Data Assimilation version 3.3.2 (SODA v3.3.2; Carton et al., 2018). The observational datasets used for this study span 1980 to 2018, except for ERSST, which covers the period from 1900 to 2018.

In this study, we employ regression and composite analyses as the primary tools to investigate ENSO features. We represent El Niño using indices over crucial regions, including Niño 3 (5^∘ N–5^∘ S, 150–90^∘ W), Niño 3.4 (5^∘ N–5^∘ S, 170–120^∘ W), and Niño 4 (5^∘ N–5^∘ S, 160^∘ E–150^∘ W), in the regression analysis. For the observational Niño 3.4 index, we use a base period between 1900 and 2014 from ERSSTv5, following the Niño index calculation of the Climate Prediction Center, NOAA (NOAA Climate Prediction Center, 2000). We also use model data from 1900 to 2014 as the base period for TaiESM1's historic run. To avoid impacts of model bias on longer timescales, such as interdecadal variation, we utilize the full length of available simulation data to obtain the most robust statistics of ENSO feature simulated by TaiESM1.

We choose the composite method instead of the regression map to better identify teleconnection signals associated with the El Niño events. In our preliminary analysis, the regressed maps of El Niño events show patterns similar to the composite events in the tropics, but with much weaker signals in the midlatitudes (not shown). Furthermore, as the ENSO events simulated in TaiESM1 show very symmetric alternations between El Niño and La Niña events, our composite based on El Niño events is followed by La Niña events. To build the El Niño composite, we choose strong events with the Niño 3.4 index larger than 1 standard deviation (i.e., 1.22 ^∘C) over the simulation period. The ERSSTv5 dataset includes eight El Niño events in 1982, 1986, 1987, 1991, 1994, 1997, 2002, and 2009. In comparison, TaiESM1 simulated 21 events throughout the entire historical simulation. As the El Niño events simulated by TaiESM1 exhibit a very strong amplitude, most of the composite fields of these events passed the significance test at the 95 % confidence level (not shown). Therefore, we will not denote regions passing significant tests in the composite fields in the following analysis.

In our analysis of TaiESM1's ability to simulate ENSO diversity, we examined its capability to distinguish between eastern Pacific (EP) and central Pacific (CP) El Niño events using the Niño 3–Niño 4 approach (Kug et al., 2009). During the historical period of TaiESM1, we identified 23 CP events and 17 EP events. This higher frequency of CP events in TaiESM1 is consistent with previous findings in CMIP models (Capotondi et al., 2020b; Chen et al., 2017; McPhaden et al., 2011). It is worth noting that the total number of El Niño events exceeds 21 when using the El Niño 3.4 index, as some EP events transition into CP events during the mature phase. In contrast, during the observation period from 1980 to 2014, we found a total of four EP events and four CP events. In the subsequent sections, we will also discuss the composites of EP and CP events in further detail.

3 Basic statistics of ENSO in TaiESM1

3.1 Niño 3.4 SST variability

Figure 1 shows the mean SST state (white contour) and monthly standard deviation (color shading) in ERSST and TaiESM1 over the tropical Pacific. The monthly standard deviation denotes the deviation of monthly SSTs from long-term monthly mean SST. TaiESM1 has a much more substantial equatorial Pacific SST variability, with an elongated region of high SST variation extending further into the warm pool region. Such a westward extension is collocated with the equatorial cold tongue, indicated by the 27 ^∘C isotherm thick white contour laying over the tropical eastern Pacific in the climatological SST mean field. TaiESM1 simulated a westward extension of the cold tongue compared with ERSSTv5. TaiESM1 also overestimated SST variation over other tropical oceans, including the Indian Ocean and warm pool.

Figure 1 Mean SST (white contours) and SST monthly standard deviation (color shading) for (a) detrended ERSSTv5 over 1958–2014 and (b) 1900–2014 of the TaiESM1 historical simulations. The contour interval is 3 ^∘C, and the thick white contour indicates the 27 ^∘C isotherm.

Figure 2 shows the Niño 3.4 SST index based on the ERSSTv5 and TaiESM1 historical simulations. TaiESM1 has more oscillatory ENSO signals with alternating cold and warm phases during 3–5 years compared with the observations. The Niño 3.4 index's standard deviation of ERSSTv5 is 0.84 ^∘C, whereas that of TaiESM1 is 1.22 ^∘C. The simulated ENSO amplitude decreased after 1980 when the global temperature increased (Fig. 2a, b). When comparing SST and surface wind fields between two 30-year periods, specifically 1950–1980 and 1984–2014, during which TaiESM1 exhibits distinct ENSO variability, Fig. S1 demonstrates a notable shift in the background state towards a La Niña-like state. This shift is characterized by an increased zonal temperature gradient over the tropical Pacific and strengthened trade winds during the period of 1984–2014 in comparison to the period of 1950–1980. Previous studies investigating the ENSO response to changes in the observed mean state have indicated that such an increase in zonal wind stress can lead to a weakening of feedback mechanisms associated with El Niño (Fedorov et al., 2020; Zhao and Fedorov, 2020). This aligns with the observed decrease in ENSO variability in TaiESM1 during the period of 1984–2014.

The power spectrum of Niño 3.4 confirms what is found in the time series of Niño 3.4, namely a much larger amplitude of ENSO in TaiESM1 (Fig. 2c). The amplitude of major peak between 3 and 4 years is around 250 ^∘C² per month, while that of observed peak is around 75 ^∘C² per month in ERSSTv5. Similarly, around the secondary peak with a period of 5 to 6 years in observations, TaiESM1 also shows two spectral peaks at 6 years and 8 years with stronger amplitude, respectively. Such model bias in representing ENSO-related spectral peaks has long been noticed in CMIP models and is still one of most challenging questions for climate models (Jha et al., 2014). Figure 2d shows the seasonal cycle of ENSO SST variance in ERSSTv5 and TaiESM1. Compared with the observations, the peak months simulated by TaiESM1 occurred in boreal winter, with a 1-month delay from the observed and a larger amplitude that was 1.5 times the observed value.

Figure 2 Normalized time series of Niño 3.4 index in (a) ERSSTv5 (1900–2018) and (b) the TaiESM1 historical run (1900–2014). The standard deviation of each dataset is noted on the upper-right side of the panels. The corresponding spectrums are shown in (c), with a black line for ERSSTv5 and a blue line for TaiESM1. (d) The seasonal cycle of SST variance in ERSSTv5 (black) and TaiESM1 (blue).

3.2 Atmosphere–ocean coupling of ENSO

As a coupled oscillation system, the coupling of atmosphere and ocean plays an important role in ENSO dynamics. To see how such coupling is simulated in TaiESM1, Fig. 3 shows the regressed rainfall, wind stress, and SST to the Niño 3.4 index in the observations (Fig. 3a) and TaiESM1 (Fig. 3b). In Fig. 3a, the observation shows the west–east displacement of wind stress and warm SST over the tropical Pacific. The strong wind stress at 160^∘ E is collocated with rainfall due to (mostly meridional) moisture convergence at the west and north edges of the warm SSTA, which was located in the central-eastern equatorial Pacific and was approximately 1^∘ higher than the western equatorial Pacific. It is important to note that the major sea surface temperature anomaly (SSTA) did not occur in the eastern equatorial Pacific, where interannual variance was the highest. Instead, the SSTA occurred to the west of the region with the maximum variance. In TaiESM1 (as shown in Fig. 3b), the warm SST center is located at approximately 120^∘ W, which is further east than in the observational data. Moreover, the magnitude of the SST anomaly is approximately 50 % greater in TaiESM1 than in the observations. Compared with ERSST, the warm SSTA in TaiESM1 is meridionally narrower and more zonally elongated to the western Pacific around 155^∘ E, causing stronger zonal and meridional SST gradients. An increase in the meridional SST gradient induces stronger meridional wind and moisture convergence over the equatorial Pacific. Zonally, the stronger westerly wind extends wider from 155^∘ E to 120^∘ W near the eastern edge of the island of New Guinea. As a result, TaiESM1 produces strong wind stress and more deep convection (shown by rainfall; color shading in Fig. 3) to the north and west sides of the warm SSTA. Figure S2 shows the regressed magnitude of wind stress onto the Niño 3.4 index and marks longitude center with dashed lines. While TaiESM1 reproduces the longitudinal center of wind stress response at 140^∘ E as in observations, the response magnitude of wind stress to SST increase is weaker in TaiESM1 than in the observations.

Figure 3 Regression map of precipitation (mm d⁻¹, color shading), wind stress (s⁻¹, vectors), and sea surface temperature (SST; ^∘C, contours) on the normalized Niño 3.4 index for (a) the Global Precipitation Climatology Project (GPCP), Simple Ocean Data Assimilation (SODA), Extended Reconstructed SST (ERSST), and (b) TaiESM1 historical simulation.

Figure 4 shows the regressed SSH to the wind stress averaged over the Niño 3 region (5^∘ N–5^∘ S, 150^∘ E–90^∘ E) to show the thermocline response to the strengthening of equatorial wind stress over the western Pacific. In the observation, a west–east dipole of thermocline response is found in Fig. 4a, showing thermocline deepening over the eastern Pacific (marked as black square in Fig.4) and shallowing over the western subtropical Pacific. Compared with the observations, TaiESM1 has captured this west–east dipole of SSH response to equatorial wind stress, but with a much stronger magnitude over the eastern equatorial Pacific (Fig. 4b). Such a strong response indicates that the SSH in TaiESM1 is more responsive than the observed to the wind stress and can easily lead to an El Niño state through the Bjerknes feedback by reducing the zonal SST gradient when the wind stress anomaly in the central equatorial Pacific initiates.

Figure 4 Regression map of sea surface height (SSH, cm N⁻¹ m², color shading) upon the normalized wind stress averaged over the Niño-4 region of 5^∘ S–5^∘ N and 160^∘ E–160^∘ W (a) from Simple Ocean Data Assimilation (SODA) and (b) the historical run of TaiESM1. The black square represents the Niño 3 region.

In addition to the two components of atmosphere–ocean coupling related to the wind stress, we also examine the heat flux–SST coupling related to ENSO. Figure 5a and b show the shortwave radiation regressed to the Niño 3.4 index in MRE2 and TaiESM1. In the observations, the reduction in shortwave fluxes prevails in the tropics with an increased warm SSTA in the Niño 3.4 region because of emerging deep convection reflecting more shortwave radiation back (Fig. 5a). A zonal gradient of shortwave fluxes is shown from −160 W m⁻² K⁻¹ over the western Pacific (i.e., 170^∘ E) to −80 W m⁻² K⁻¹ over the eastern Pacific (i.e., 120^∘ W). Overall, TaiESM1 reproduced the negative feedback patterns in the deep tropics (3^∘ S–3^∘ N), but with a much stronger shortwave reduction of −200 W m⁻² K⁻¹ over the western Pacific in response to the warm SSTA in the Niño 3.4 region (Fig. 5b). Such a pattern is consistent with a stronger rainfall response of TaiESM1 in Fig. 3b, suggesting that the stronger deep convection response reflects more shortwave fluxes over the western Pacific. In contrast, we observed an increase in downwelling shortwave flux over the subsidence regions adjacent to the Intertropical Convergence Zone (ITCZ) at both 10^∘ N and 10^∘ S over the eastern Pacific. One notable feature in TaiESM1 is that the shortwave fluxes can increase up to 60 W m⁻² K⁻¹ over the tropical eastern Pacific (i.e., 120 to 100^∘ W), in contrast to a decrease in observations (Fig. 5). Such a difference suggests that there is a biased cloud radiative response over the eastern Pacific when an El Niño event occurs, which may induce biased heat flux–SST coupling in TaiESM1. We will further examine and discuss this bias in Sect. 4.

Figure 5 Regression map of tropical surface downwelling shortwave radiation (RSDS; W m⁻²; color) upon the normalized Niño 3.4 index for (a) the MRE2 ensemble (1980–2017) and (b) TaiESM1 historical simulation (1900–2014).

3.3 Composite of El Niño structure and teleconnection

To evaluate the structure of El Niño events in TaiESM1, we compile strong El Niño events with a Niño 3.4 index larger than 1 standard deviation of the entire time series (i.e., larger than 1.22 ^∘C). Under this definition, there are 21 Niño events in the TaiESM1 historical run and 9 events in the MRE2 ensemble from 1980 to 2015.

Figures 6 and 7 present the seasonal variation of El Niño events in the tropics and its teleconnection pattern in the midlatitudes in the MRE2 ensemble and TaiESM1 by showing the 2 m surface temperature (color shading; Fig. 6), sea level pressure (SLP; contours; Fig. 6), precipitation (color shading; Fig. 7), and 300 hPa stream function (contours; Fig. 7). Overall, TaiESM1 reproduced the observed spatial structures and teleconnection patterns associated with El Niño; however, consistent with the over-simulated El Niño signals, TaiESM1 produces much stronger tropical SST warming and teleconnection in extratropical regions than the observations in all four seasons. As early as June, July, and August in the first year (JJA⁰), TaiESM1 already simulates an SST anomaly of 2 ^∘C over the eastern Pacific (Fig. 6a, e) with a clear rainfall response in the central Pacific (Fig. 7a, e). A zonal dipole of surface temperature and rainfall between the eastern and western Pacific forms earlier than in the observations (color shading in Figs. 6a, e and 7a, e). In SON⁰, the warm SST anomaly grows even stronger and expands over the entire tropical Pacific in TaiESM1. As a result, a very clear teleconnection similar to that of DJF⁺¹ can already be found in the Northern Hemisphere, including horseshoe-shaped cooling in the western Pacific, and over Eurasia and the United States (Fig. 6b, f). In the meantime, TaiESM1 captures the responses in the Southern Hemisphere in DJF⁺¹, including the warm and dry response over northern South America, the opposite responses over southern South America, and the north–south dipole between Eurasia and South Asia (Fig. 6c, g). As for the rainfall response, TaiESM1 realistically simulates the shift of deep convection from the western Pacific to the central Pacific as the warm SST occurs in DJF⁺¹; however, the westward shift of tropical SSTA causes the surface temperature response pattern to also shift westward, resulting in enhanced stronger cooling in East and Southeast Asia. The cooling further extends into the Indian Ocean, causing an Indian Ocean Dipole-like response as depicted in Fig. 6g. In contrast to the weaker SSTA and surface temperature impacts in the observations (Fig. 6d) in MAM⁺¹, a strong teleconnection pattern in the surface temperature over the extratropical regions is sustained into MAM⁺¹ in TaiESM1 (Fig. 6h). Furthermore, the over-response of the rain band over the Indian Ocean could be due to the mean rainfall biases simulated in TaiESM1, as noticed by Wang et al. (2021). In terms of the atmospheric circulation anomaly, TaiESM1 successfully captures the southward Rossby wave propagation from the central Pacific to the southeastern Pacific (contours in Figs. 6a, e and 7a, e) in JJA⁰. From SON⁰ and into DJF⁺¹, the teleconnection in TaiESM1 intensifies as the warm SST over the equatorial Pacific develops into the mature stage of El Niño (Figs. 6b, f; 7b, f). TaiESM1 reproduces the Rossby wave train response emitted from the equatorial Pacific into North America and the resulting dipole of surface temperature over North America during DJF⁺¹ (Fig. 6c, g). In line with the stronger temperature and rainfall responses observed during MAM+1, TaiESM1 exhibits El Niño-related circulation anomalies across the tropics (as shown in Figs. 6d, h and 7d, h).

Figure 6 The surface temperature (color shading) and SLP (contours; contour interval is 1 hPa; contours smaller than zero are dashed) of the El Niño composite in (a, e) JJA⁰, (b, f) SON⁰, (c, g) DJF⁺¹, and (d, h) MAM⁺¹ based on the MRE2 ensemble (left column) and TaiESM1 historical run (right column).

Figure 7 As in Fig. 6 but showing precipitation (mm d⁻¹; color shading) and 300 hPa stream function (contours; contour interval of 2 × 10⁶ m² s⁻¹ anomalies; contour values smaller than zero are dashed).

Figure 8 shows the seasonal evolution of ocean subsurface potential temperature averaged over the selected El Niño events. The composites show four seasonal means from JJA⁰ when an El Niño event was identified to MAM⁺¹ in the following year. The green line shows the location of the 20 ^∘C subsurface isotherm (Z20) during El Niño events, and the gray line shows the climatological Z20. During JJA⁰ in the observations, Z20 deepens in the eastern Pacific and shallows in the western Pacific (Fig. 8a). While TaiESM1 realistically simulated the climatological Z20 depth, it overestimated the response of subsurface temperatures and simulated a flatter Z20 profile in JJA⁰ (Fig. 8e). Such an overestimation bias in TaiESM1 continues from the beginning of the ENSO evolution through SON⁰ and DJF⁺¹ (Fig. 8b, c, f, and g). Accompanied by the warm bias over the eastern Pacific, the cold bias developed in SON⁰ when the cool water started to form in the tropical western Pacific (Fig. 8b–d, f–h). Such a zonal dipole of subsurface temperature bias in TaiESM1 manifests a basin-wide response of ocean circulations during El Niño events. Moreover, both warming and cooling from the surface to 100 m depth were much more pronounced in the model than in the observation, especially over the eastern Pacific. An unrealistic warming in the central-eastern equatorial Pacific is also notable, reflecting the unrealistic westward extension of positive SSTA. This bias led to an anomalously strong SST gradient up to 2 ^∘C between 180 and 150^∘ W, consistent with the widespread strong westerly anomalies. In the meantime, the composite of subsurface zonal currents corresponding to El Niño events simulated in TaiESM1 is shown in Fig. S3. A strong westerly current anomaly suggests that zonal advection may also play a role in driving strong El Niño (Fig. S3). More analysis is needed to determine which model components are more responsible for these biases.

Figure 8 Equatorial cross-section (5^∘ S–5^∘ N) of the El Niño composite of the potential temperature anomaly (color shading) in (a, e) JJA⁰, (b, f) SON⁰, (c, g) DJF⁺¹, and (d, h) MAM⁺¹ based on SODA v3.3.2 (left column) and the TaiESM1 historical run (right column). The gray line shows the climatological 20 ^∘C isotherm (Z20), and the green dashed line shows the Z20 in the Niño state.

3.4 ENSO diversity and teleconnection

Based on the EP and CP events identified by the Niño 3–Niño 4 approach, we make longitudinal profiles of SSTA for both EP events and CP events simulated in TaiESM1 compared to observations (Fig. S4a). Our analysis reveals that TaiESM1 generally exhibits warmer SSTA over the tropical Pacific, particularly east of the 150^∘ E line. Notably, during EP events, the model shows SSTA that can be as high as 1 ^∘C over the eastern Pacific region. To further analyze the impacts of EP and CP events, we constructed composites of surface temperature and SLP based on four EP events and four CP events during the observation period (1980–2014), as well as CP and EP events during the historical TaiESM1 period (1900–2014).

Regarding the EP composite, TaiESM1 successfully captures the observed features over the tropical Pacific, as illustrated in Fig. S4b and c. However, the CP events identified in TaiESM1 exhibit elongated warm SSTA in the tropical region but with weaker teleconnections to the midlatitudes in the Northern Hemisphere (Fig. S4d and e). Additionally, the warming over North America is less pronounced and retreats towards the polar region in TaiESM1, whereas the observed cold surface temperature anomaly is replaced by a warm anomaly. These discrepancies suggest that biases in the model's mean state contribute to the model's biases in ENSO diversity and teleconnection patterns observed in TaiESM1, a common issue seen in other climate models as well (Ham and Kug, 2012).

4 Linking the warm SST bias during El Niño events with simulated seasonal mean states

In this section, we analyze the seasonal life cycle of the strong ENSO signals simulated in TaiESM1 to understand the intense tropical warming anomaly of El Niño events and link this bias with seasonal mean states. Based on the El Niño events defined in the previous section, we construct the seasonal cycle of El Niño events by plotting a Hovmöller diagram averaged over the Equator between 3^∘ S and 3^∘ N.

Figure 9 shows the Hovmöller diagram of SST anomalies and 1000 hPa zonal winds along the Equator (3^∘ S–3^∘ N) based on the strong El Niño event selected from the ERSSTv5 and TaiESM1 historical run. In the observation in May, the warm SSTA occurs over the dateline and the westerly wind anomaly starts to propagate to 135^∘ W. In JJA⁰, the warm SST slowly develops at 180–135^∘ W, progressively propagates to the eastern equatorial Pacific with westerly anomalies, and reaches maximum amplitude in November, December, and January (Fig. 9a). In contrast, the warm SST when initiated in May intensifies almost simultaneously in the basin east of 150^∘ W in TaiESM1 (Fig. 9b). The warm anomaly quickly reaches 2 ^∘C over the entire eastern equatorial Pacific through JJA⁰ and into September. Such warming seems to be coupled, with one branch quickly extending to the eastern equatorial Pacific after the initiation in May, whereas the westerly anomaly's major branch is well coupled with SSTA over the western Pacific most of the time. The warming continued developing even when the westerly anomalies in the eastern equatorial Pacific weakened in JJA⁰ and reached a maximum in January as in observations. This early development of a warm SSTA in the eastern equatorial Pacific in May and continuous warming in JJA⁰ and SON⁰ are two primary model biases that may contribute to strong El Niño events in TaiESM1.

Figure 9 Hovmöller diagram of composite SST anomalies (color) and zonal wind at 1000 hPa (contour) along the Equator (3^∘ S–3^∘ N) for El Niño based on (a) observations (ERSST–MRE2 ensemble from 1980 to 2015) and (b) the TaiESM1 historical run (1900–2014).

To understand the heat flux–SST coupling, we examined similar composites of net surface heat flux (Fig. S5a, b) and latent heat fluxes (Fig. S5c, d), which are the two heat fluxes prominent in observational El Niño events. Notably, in the observations, the seasonal variation of surface net heat fluxes in MRE2 is primarily controlled by latent heat fluxes, especially from May to November, as other heat fluxes play a minor role (Fig. S5c, d). However, the spatial patterns of net surface fluxes of TaiESM1 are primarily dominated by shortwave surface fluxes (Fig. S5d) and amplified by latent heat fluxes (Fig. 2b). Figure 10 shows the same Hovmöller diagram but for downwelling surface shortwave flux (color shading) and surface temperature (contours) composited over the El Niño events. Evolutions of shortwave radiation and SST during the life cycle of El Niño events in observations and TaiESM1 are shown in Fig. 10a and b, respectively. In the observations, upward shortwave fluxes of about 10 to 20 W m⁻² (i.e., shown as negative in Fig. 10) are seen during the entire El Niño period, which is well collocated with the warm SSTA, due to the increased shortwave reflection of deep convection triggered by warmer SST. Such reflection of shortwave fluxes by deep clouds induces negative feedback (i.e., higher SST, large reflected shortwave radiation) between shortwave radiation flux and SST in the tropics (Fig. 10a). The negative feedback intensified in September, reached the maximum in December, and continued into the following February (Fig. 10a). In contrast to the observations, TaiESM1 produced unrealistically strong negative feedback over the western equatorial Pacific and Intertropical Convergence Zone (ITCZ) regions (Fig. 10b), as well as very strong positive shortwave radiation anomalies near the equatorial eastern Pacific after May in year 0 (Fig. 10a). The increase in shortwave influx starts from the eastern Pacific in May and gradually extends to 125^∘ W in November, a feature that was not seen in observations. The westward progression of shortwave radiation increase contributes to erroneously strong SST warming, and its westward extension is seen in TaiESM1 from May to Jan in year 1 (Fig. 10b) instead of the observed simultaneous warming in the eastern Pacific.

Figure 10 A composite of El Niño events for RSDS (W m⁻²; color) and surface temperature (^∘C; contour) in (a) MRE2 and (b) TaiESM1. Both variables are averaged over the Equator (3^∘ S–3^∘ N).

Figure 11 Climatological seasonal cycle of the surface temperature (shaded area; ^∘C), low clouds (gray contours; %), and 500 hPa vertical velocity (green contours; hPa s⁻¹) for the tropical Pacific (3^∘ S–3^∘ N) in the (a) MRE2 ensemble and (b) TaiESM1, as well as (c) their differences. (d) The seasonal cycle of surface temperature at 100^∘ W in MRE2 (black) and TaiESM1 (blue).

Figure 12 Longitudinal height cross-section of the El Niño composite of September, October, and November in year 0 (SON⁰) along the Equator (5^∘ S–5^∘ N) based on (a) ERA5 and (b) TaiESM1. The cloud fraction (%) of the El Niño composite is shown in color shading, and the zonal and vertical winds (units: hPa s⁻¹) are shown as green vectors.

To understand the cause of the strong SSTA development from May to November identified in Fig. 9 and its relationship with shortwave radiation fluxes in Fig. 10, we analyzed the seasonal cycle of low clouds (gray contours), vertical velocity (green contours), and SST (color shading) in the observations and TaiESM1 in Fig. 11a and b. Figure 11c shows the differences between TaiESM1 and MRE2. Figure 11d shows the climatological seasonal SST cycle at 100^∘ W to show the SST differences where cloud biases are prominent. On the seasonal timescale in observations, the SSTA is closely coupled with anomalies of vertical velocity and low-level clouds. Cold surface temperature is commonly collocated with excessive low-level clouds and subsidence (Fig. 11a). While TaiESM1 also exhibits a clear seasonal variation like the observations in the eastern tropical Pacific (Fig. 11b), the simulated seasonal cycle is rather asymmetric with a colder bias during MAM⁰ and SON⁰ and warm biases during JJA⁰ and DJF⁺¹ (Fig. 11c). Compared with the observations, TaiESM1 warmed approximately 1 month later during MAM⁰ and cooled deeper over the eastern Pacific in SON⁰ (Fig. 11c). The warmer MAM⁰ sea surface in TaiESM1 provides a smaller zonal SST gradient and may lead to the earlier onset of the eastward-propagating westerly anomaly found in MAM⁰ through the Bjerknes feedback (Fig. 9b). However, this cold surface temperature bias during SON⁰ provides a cold lower boundary for low stratus clouds to develop, leading to an environment for the strong positive feedback between shortwave fluxes and SST during El Niño events. Such impacts of the cold tongue bias are also found in CESM1 and CESM2, which share the same ocean model (i.e., POP2) as TaiESM1 (Y. C. Wang et al., 2015; Wei et al., 2021).

For a further look into the bias in SON⁰, Fig. 12 shows the changes in cloud fraction and zonal circulations in SON⁰ during the composited El Niño events. TaiESM1 successfully reproduced the eastward shift of the upward branch related to deep convection in response to the warm SST during El Niño events (Fig. 12a, b). However, TaiESM1 produced a much stronger response in circulations and cloud cover than the observations (Fig. 12b). Especially over the eastern Pacific, a stronger upward motion anomaly occurs and about 10 %–20 % of low clouds are reduced in TaiESM1, indicating a dramatic reduction of the low-cloud regime (Fig. 11c). As a result, more shortwave heat influxes are allowed into the ocean surface in SON⁰ when El Niño events occur (Fig. 10b). Consistent with the seasonal variation of SSTA in Fig. 9b, the increased shortwave fluxes thus help to keep the warm SSTA in SON⁰ over the eastern Pacific after the warm SSTA initiates in MAM⁰ in TaiESM1 (Fig. 10b).

5 Summary

This study documented ENSO's fundamental statistics and features in TaiESM1, a CMIP6 participant. Compared with the observational dataset, TaiESM1 has captured many prominent observed ENSO features, including a 3–5-year spectrum peak, seasonal phasing, the evolution of warm SST, deepening of the subsurface layer during an El Niño event, and teleconnection patterns in midlatitudes. However, the simulated El Niño signals in TaiESM1 are much stronger and more prolonged than the observed ENSO signals. Such strong signals are shown as intense warm SSTA over the tropical Pacific within the spatial structure of the composite of the El Niño event. In the meantime, in response to the tropical warm SSTA, the teleconnection wave activities associated with El Niño events are much stronger in TaiESM1, significantly affecting temperature and rainfall in high-latitude regions, although with a similar spatial pattern to the observations.

To understand the cause of the warm El Niño biases, we investigated the seasonal cycle of strong El Niño events in TaiESM1. Compared with the observed delay of SST propagation to the central Pacific in May after the wind anomaly's initiation, the simulated SST warming quickly propagates toward the eastern Pacific along with the westerly wind anomaly in TaiESM1. The SST warming continued through November and reached a maximum in December, showing a stronger and more prolonged warm period than the observation. Moreover, in contrast to the negative feedback found in observations, our analysis shows that the strong El Niño warm anomalies in TaiESM1 is due to the strong SST–heat flux positive feedbacks, especially over the eastern Pacific in SON⁰. Further analysis shows that this biased feedback is due to the response of the spurious low-cloud regime off the west coast of South America simulated in TaiESM1. This biased cloud regime results from the seasonal variation of the cold tongue over the eastern Pacific, consistent with previous studies using the CESM family and CMIP models (Ham and Kug, 2012; Wei et al., 2021). During El Niño events, the stratus clouds over the eastern Pacific gradually diminished due to the warmer SST, allowing solar radiation to warm the ocean surface. This result leads to a positive feedback of downward solar radiation and SST over the eastern Pacific, an opposite sign of the observed relationship. This biased relationship is typical in the CMIP5 and CMIP6 models (Bayr et al., 2019; Beobide-Arsuaga et al., 2021).

In summary, TaiESM1 can reproduce many fundamental features of ENSO. However, it still possesses several biases shared by other CMIP6 models, including the lack of randomness of El Niño events, El Niño magnitudes that are too strong, and early SST warming in the early stage of El Niño events. The strong El Niño strength, our analysis found, mainly results from the biased atmosphere–SST coupling accompanied by biases of the mean state and seasonal cycle of the cold tongue and Walker circulation, which is consistent with previous studies with CMIP5 and/or CMIP6 models (Bayr et al., 2018; Chen et al., 2021). However, to resolve this bias (and others), more detailed analysis, including process-related metrics, and more model experiments, such as atmosphere-only and ocean-only experiments, are required to determine the cause and effects of the observed biases. For TaiESM1, we plan to implement ocean-only experiments with the ocean component POP2, allowing us to quantify the ocean's response to biased winds and radiation fluxes. We will also conduct AMIP-type simulations to investigate the development of westerly wind anomalies under biased SST conditions. Combined with process-oriented diagnosis, these model experiments will help us to determine and better comprehend the causes and effects of these observed biases. Also echoed with other studies of ENSO evaluations of CMIP6, our analysis suggests that a more process-based development strategy focusing on atmosphere–ocean coupling rather than a feature-based evaluation of ENSO is needed to reduce the uncertainty of ENSO simulations and future ENSO projection in climate models.