Changes in El Niño characteristics and air–sea feedback mechanisms under progressive global warming

In this study, we investigate the potential changes of El Niño characteristics, including intensity, frequency and CP/EP El Niño ratio, under progressive global warming based on the 140-year CMIP6 model simulation outputs with the 1pctCO2 experiment. Major air-sea feedback mechanisms attributing to the changes are also examined. The CMIP6 ensemble means project a slight enhancement of El Niño intensity by about 2% and a modest increase of El Niño frequency by about 4% from the first to the second 70-year periods. It is found that these small changes result from the opposite response to global warming between CP and EP El Niño, i.e., the intensity of EP El Niño is projected to weaken by nearly 4.6% while the intensity of CP El Niño is projected to increase by about 4.5%. Since CP El Niño occurs more frequently than EP El Niño in CMIP6 simulations, this leads to a slight enhancement of the total El Niño intensity if these two types of El Niño were not separated. A similar situation occurs in projecting the future change of El Niño frequency, i.e., the frequency of EP El Niño is projected to decrease by about 1.4% while the frequency of CP El Niño is projected to increase by about 2%, thereby leading to a modest increase of the total El Niño frequency. By comparing the variance explained by key air-sea feedback mechanism between the two 70-year periods, we also note that the increased CP/EP ratio can be explained by the enhanced role played by the SF (seasonal footprinting) mechanism in a warmer atmosphere. Our study also points out that, as long as a climate model can correctly produce the intensity (variance) of major air-sea feedback mechanisms, the relationship between changes in El Niño characteristics and changes in feedback mechanisms can be physically robust. The potential changes of El Niño characteristics are investigated based on 140-year CMIP6 simulation outputs with the 1pctCO2 experiment. The ensemble means project a slight enhancement of El Niño intensity (2%) and a modest increase of El Niño frequency (4%) from the first to the second 70-year periods. CP El Niño becomes stronger and more frequent and CP/EP ratio increases due to the enhanced SF mechanism in a warmer atmosphere. The potential changes of El Niño characteristics are investigated based on 140-year CMIP6 simulation outputs with the 1pctCO2 experiment. The ensemble means project a slight enhancement of El Niño intensity (2%) and a modest increase of El Niño frequency (4%) from the first to the second 70-year periods. CP El Niño becomes stronger and more frequent and CP/EP ratio increases due to the enhanced SF mechanism in a warmer atmosphere.


Introduction
The El Niño-Southern Oscillation (ENSO) has long been recognized as the dominant natural climate variability arising from strong atmosphere-ocean interaction in the tropical Pacific basin (Bjerknes 1969;Philander 1985;Trenberth 1997).Several theories have been proposed to account for the excitation and maintenance of ENSO (Neelin et al. 2009;Wang 2018;Fang and Xie 2020).Bjerknes (1969) first hypothesized that interaction between atmosphere and ocean is essential to the establishment of ENSO.He recognized that the positive feedback between trade winds and zonal sea surface temperature (SST) gradients across the tropical Pacific would lead to the creation of ENSO.Although the Bjerknes feedback theory is accountable for the establishment of El Niño or La Niña phase, it cannot explain why the anomalous El Niño or La Niña conditions would end, how long they would last and why they had their peak during the boreal winter (i.e., El Niño or La Niña phase locked to the seasonal cycle).
Based on the Zebiak and Cane's (1987) coupled model framework, Jin (1997aJin ( , 1997b) ) proposed the rechargedischarge oscillator theory for ENSO.Specifically, the warm phase of ENSO is accompanied by westerly wind anomalies in the central Pacific and warm SST anomalies in the eastern Pacific.The resulting wind stress curl anomalies then lead to a divergence of the Sverdrup transport and a discharge of the tropical oceanic heat content.The discharge of oceanic heat content, in turn, generates an uplifting of the oceanic thermocline in favor of the transition from warm phase to cold phase of ENSO.Conversely, the same process but with an opposite sign will lead to the transition from cold phase to warm phase of ENSO.In short, the equatorward or poleward transport of oceanic heat content below the Ekman layer plays the key role in the recharge-discharge oscillator theory.
From another point of view, Vimont et al. (2003) proposed the seasonal footprinting theory to describe how the subtropical Pacific SST anomalies spread into the equator through the subtropical ocean-atmosphere coupling over several seasons to trigger an El Niño event.The seasonal footprinting signal generally starts from positive SST anomalies in the subtropical northeast Pacific, with a spatial pattern close to the Pacific Meridional Mode (Chiang and Vimont 2004), Victoria Mode (Ding et al. 2015), North Pacific Gyre Circulation (Di Lorenzo et al. 2015), or the North Pacific Meridional Mode (Jia et al. 2021).These warm SST anomalies would lead to southwesterly wind anomalies in the subtropical northeast Pacific that are opposing to the mean northeasterly trade winds, leading to the weakening of trade winds and surface evaporation.The reduced surface evaporation, in return, warms the SST on the equatorward side of the SST pattern through the windevaporation-SST feedback (Xie and Philander 1994), eventually triggering an El Niño event in the equatorial Pacific (Yu and Kim 2011;Anderson and Perez 2015).In contrast, negative seasonal footprinting signals (i.e., negative SST anomalies in the subtropical northeast Pacific) would trigger a La Niña event.
Over the past few decades, increasing attention has been paid to exploring the diversity of El Niño.For instance, Fu et al. (1986) noticed that there are two main El Niño patterns.The first one shows warm SST anomalies in the east and central Pacific areas and slightly below normal in the west.The second one exhibits nearly uniformly warm SST anomalies in the entire Pacific area.Xu and Chan (2001) later found that the onset time of the two types of El Niño behave differently.The eastern and central warming El Niño typically started in autumn; while the uniformly warming El Niño generally started in summer.They thus referred the two distinct patterns as the spring-onset type El Niño and summer-onset type El Niño, respectively.At about the same time, Larkin and Harrison (2005) pointed out that many El Niño events' warm SST anomalies occur near the central Pacific dateline and they named this type of El Niño as the dateline El Niño.The dateline El Niño was found to exhibit different temperature and rainfall variations from the traditional eastern Pacific El Niño.
While studying El Niño diversity has been a popular issue over the past few decades, a well-posed objective definition for measuring the diversity of El Niño did not exist until recently.Ashok et al. (2007) utilized the empirical orthogonal function (EOF) to identify two different types of El Niño, namely the typical and atypical El Niño, from the first and the second EOF modes, respectively.Their study showed that the atypical El Niño (also referred to as the El Niño Modoki) might influence the temperature and precipitation over many parts of the globe that are generally opposite to those of the typical El Niño.Kao and Yu (2009) later used the regression-EOF analysis to provide inductive integration of the two types of El Niño.Their study pointed out that the eastern Pacific El Niño (hereafter EP El Niño) is closer to the traditional El Niño which is caused mainly by large-scale air-sea coupling that involves changes in Walker circulation and thermocline across the equatorial Pacific.On the other hand, the central Pacific El Niño (hereafter CP El Niño) is caused by local air-sea interaction that shows changes of SST, surface winds and oceanic surface layer in the central Pacific.Since CP El Niño began to develop from the shallow oceanic surface layer, local atmospheric forcing is considered to play a more important role in exciting this type of El Niño.
Applying the multi-variate EOF (MEOF) analysis technique on the SST, surface winds and sea surface height (SSH) data during the historical period from 1958 to 2014, Yu and Fang (2018) recently noted that major atmosphere-ocean feedback mechanisms associated with El Niño can be captured by the first three MEOF modes.Their study further concluded that the Bjerknes feedback is a positive feedback process describing the enhancement of El Niño intensity.The seasonal footprinting mechanism involves extratropical forcing which tends to trigger an El Niño (or La Niña) in the equatorial central Pacific, thereby producing a more episodic and irregular El Niño evolution.The recharge-discharge oscillator mechanism, on the other hand, tends to force El Niño and La Niña to follow each other, accordingly generating a more cyclic El Niño evolution.
In the last 10 years or so, simulation outputs from the Coupled Model Intercomparison Project (CMIP) phase 3 (CMIP3) and phase 5 (CMIP5) have been widely used to project the potential changes of El Niño characteristics under global warming.For instance, Yeh and Kug (2009) analyzed simulation outputs of 11 CGCMs (coupled general circulation models) archived in the CMIP3 database and found that the occurrence ratio of CP/EP El Niño is projected to increase as much as five times under global warming due to a flattening of the thermocline in the equatorial Pacific.Based on ensemble-mean databases from CMIP3 and CMIP5, Cai et al. (2014) noted a doubling in the occurrence frequency of extreme El Niño events in the future in response to greenhouse warming.Their study also pointed out that the increased frequency extreme El Niño events arises from a projected surface warming over the eastern equatorial Pacific, facilitating more occurrences of atmospheric convection in the eastern equatorial region.Based on observational data during 1958-2014, Yu and Fang (2018), nonetheless, found that the impact of seasonal footprinting mechanism would become more important as the atmosphere gets warmer, thus leading to more CP-El Niño events rather than EP-El Niño events as claimed by Cai et al. (2014).A recent review article of Cai et al. (2021) noted that the El Niño SST variability and extreme El Niño events are likely to increase under greenhouse warming, with a stronger inter-model consensus in the latest generation of climate models, although uncertainties (e.g., location of SST variability and time of emergence) remain.
Some other studies emphasized the important role of background mean state in driving the El Nino diversity.By analyzing outputs from a 500-year preindustrial simulation run, Choi et al. (2011Choi et al. ( , 2012) ) found that the decadal change in frequencies of EP and CP El Nino could be attributed to the interaction between El Nino and mean state.The higher occurrence regime of CP El Nino is associated with a mean state of stronger zonal SST gradients in the equatorial Pacific and stronger trade winds to the east of the dateline.Similar arguments were reported in Capotondi et al. (2015) and Wang et al. (2019), in which they linked the more frequent occurrence of extreme El Nino events to the background SST warming in the western Pacific and the increased zonal and vertical SST gradients in the equatorial central Pacific.
In the present study, we attempt to explore (i) how the El Niño characteristics, including intensity, frequency and CP/EP El Niño ratio, might change under global warming?(ii) what air-sea feedback mechanisms are possibly responsible for such changes, using the newest CMIP6 model simulations with the 1pctCO2 experiment.While the 1pctCO2 scenario may not accurately reflect the recent global warming caused by non-CO2 factors such as aerosols, deforestation and urbanization, it can be very helpful in understanding high-emission future scenarios where CO2 is the main driver of climate change.The remainder of this paper is structured as following.Section 2 describes the data and analysis methods used.Section 3 evaluates the performance of CMIP6 models in reproducing the observed SST and cloud radiative effect (CRE) anomalies over the tropical Pacific during the mature phase of historical El Niño events.Sections 4 and 5 present the potential changes of El Niño characteristics and air-sea feedback mechanisms, respectively, as projected by CMIP6 simulations.Major findings are concluded in Sect.6.

Data and methodology
Several observational and model simulation datasets are used in this study.A brief introduction of these datasets is described below.

Data sources (i) ERSST data
The Extended Reconstructed Sea Surface Temperature (ERSST) data, generated by the National Oceanic and Atmospheric Administration (NOAA), was adopted as a reference to evaluate the performance of CMIP6 models in simulating the SST anomaly pattern associated with historical ENSO events.The ERSST dataset provides monthly SST data over the global oceans from 1900 to 2000, with a horizontal resolution of 2° × 2° (Smith and Reynolds 2003).The newest fifth version of ERSST (ERSSTv5) is used (Huang et al. 2017).

(ii) CERES/EBAF data
The Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) Surface data provides computed monthly mean radiation fluxes at the Earth's surface.The CERES/EBAF product is intended for studies that measure the climate models' performance, estimate the Earth's annual global mean energy budget, and infer the meridional heat transports (Wielicki et al. 1996).The CERES/EBAF dataset consists of 1° × 1° global coverage of surface upward/downward longwave, upward/downward shortwave and net fluxes under both clear-sky and all-sky conditions for calculating the cloud radiation effect.The period of CERES/EBAF data is from March 2000 to present.

(iii) ERA5 data
The fifth generation ECMWF reanalysis atmospheric data (ERA5) combines the model simulation and observations into a globally complete dataset to replace its previous generation product-the ERA-Interim reanalysis (Hersbach et al. 2020).The horizontal resolution of ERA5 is 0.25° × 0.25° with key atmospheric parameters on 37 pressure levels.The period of ERA5 data is from 1979 to present.In this study, monthly SST, zonal (u) and meridional (v) winds are used to conduct the MEOF analysis.

(iv) SODA data
The Simple Ocean Data Assimilation (SODA) data is an oceanic reanalysis product.It was built on the Modular Ocean Model, version 5 (MOM5) from National Oceanic and Atmospheric Administration (NOAA)/Geophysical Fluid Dynamics Laboratory (GFDL) (Carton et al. 2018).The latest version SODA3.12.2 used JRA-55 as meteorological forcing to conduct data assimilation.The horizontal resolution is 0.5° × 0.5° with 50 vertical levels.The period of SODA3.12.2 is from 1980 to 2017.The monthly SSH and surface wind stress magnitude (TAU) data, including zonal and meridional components, is used to conduct the MEOF analysis.

(v) CMIP6 historical simulation data
The historical simulation outputs archived in CMIP6 program are used to evaluate the performance of climate models in simulating the observed El Niño SST warming and the associated cloud radiative effect (CRE) patterns over the tropical Pacific.The CMIP6 historical simulation was forced by time-varying external and internal conditions from observations.Both naturally forced changes (e.g.due to solar variability and volcanic aerosols) and changes due to human activities (e.g.CO2 concentration, aerosols and land use) are considered.This historical simulation may serve as an important benchmark for assessing the performance of climate models through evaluation against observations (Eyring et al. 2016).

(vi) 1pctCO2 experiment data
The 1pctCO2 experiment, which assumes an idealized increase in CO2 concentration by 1% per year, is chosen for future projection.We choose this experiment instead of the Shared Socioeconomic Pathways (SSP) scenario because this experiment has been performed in all phases of CMIPs since CMIP2 and has been served as a consistent and useful benchmark for analyzing climate sensitivity and feedback.Most importantly, it provides a long simulation period (> 140 years) that ensures us to select enough El Niño cases for analysis.
In this study, CMIP6 historical simulation outputs from January 1979 to December 2005 are adopted to compare against observations.The 140-year simulation output from the 1pctCO2 experiment is taken for analysis [see http:// www-pcmdi.llnl.gov/ (Taylor et al. 2012, andEyring et al. 2016) for more details].The monthly variables used include SST, air temperature, specific humidity, geopotential height, horizonal winds, precipitation along with surface flux variables such as sensible heat, latent heat, longwave and shortwave radiation.All variables were interpolated onto the same 2.5° × 2.5° horizontal resolution, and their linear trends were removed prior to analysis to focus purely on the natural interannual variability.

Selection and classification of El Niño
Two different approaches are adopted to define an El Niño event.For observations and historical simulations, SST anomalies averaged over the Niño 3.4 region (5° N-5° S, 170° E-120° W), also referred to as the Niño 3.4 index, during the boreal winter (December, January and February) is used.Following Trenberth (1997), an El Niño event is defined when the 3-month running average of Niño 3.4 index is greater than 0.5 °C for five consecutive months.However, one should note that the above criteria cannot be applied to future projections due to significant discrepancies in simulating amplitudes of El Niño SSTA across the CMIP6 models.Since the 0.5 °C threshold used in observations is equivalent to a standard deviation value of 0.58 for El Niño SSTA, we thus use this value as a threshold for selecting the model El Niño events in future simulations.To avoid the weighting bias caused by a few models with very strong El Niño SSTA, the simulated SSTA in each CMIP6 model is scaled by its own standard deviation value prior to conducting the ensemble mean.This process ensures all models contribute the same weight in the ensemble mean products.
Following Kao and Yu (2009) and Yu and Kim (2013), the regressed Empirical Orthogonal Function (REOF) analysis was utilized to separate the El Niño events into CP and EP types.In practice, the SST anomalies regressed with the Niño1 + 2 (0°-10° S, 80° W-90° W) index from the original SST anomalies were removed prior to conducting the EOF analysis to identify the CP El Niño events.Likewise, the SST anomalies regressed with the Niño4 (5° S-5° N, 160° E-150° W) index from the original SST anomalies were removed prior to conducting the EOF analysis to identify the EP El Niño events.The two principal components obtained from the REOF analysis represent the El Niño intensity and are defined as the CP and EP El Niño indices, respectively.As in observations, a CP (or EP) El Niño event in model simulations is defined when the 3-month running average of CP (or EP) index is greater than 0.58 standard deviation value for five consecutive months.Moreover, the longest period of an El Niño event is set at 2 years.If the warming condition continues longer than 2 years, it will be regarded as two different El Niño events.
To examine the potential changes of El Niño characteristics (including frequency, intensity and CP/EP El Niño ratio) and air-sea feedback mechanisms under global warming, we divide the 140-year simulation outputs into two periods-the first and the second 70-year periods-to ensure that we have enough model El Niño cases for composite analysis.Accordingly, changes of El Niño characteristics and air-sea feedback mechanisms are calculated as differences between the first and the second 70-year periods (latter-former).

Air-sea feedback mechanisms
As in Yu and Fang (2018), we conduct the MEOF analysis on SST, surface winds, SSH and wind stress data over the domain (30° N-20° S, 120° E-80° W) covering the historical period of 1980-2005 to identify major physical processes controlling the evolution of El Niño in CMIP6 simulations, including Bjerknes feedback (hereafter BF), seasonal footprinting (hereafter SF) and recharge-discharge (hereafter RD) mechanisms.As shown in Fig. 1. the first MEOF mode (MEOF1) (left panels), accounting for about 15.6% of the total variance, represents the BF mechanism, characterized by warm SST anomalies in tropical central-to-eastern Pacific and westerly wind anomalies to the west and an east-west sloping of thermocline (represented by the slope of SSH).The BF mechanism can determine the intensity of a developing El Niño or La Niña event.The second MEOF mode (MEOF2) (middle panels), accounting for 8.9% of the total variance, represents the SF mechanism, characterized by positive SST anomalies from the subtropical northeastern Pacific to the equatorial central Pacific and southwesterly wind anomalies against the mean trade winds.The SF mechanism is often linked to the onset of an El Niño, or La Niña, event in the central Pacific.The third MEOF mode (MEOF3), accounting for only 3.2% of the total variance, represents the RD mechanism, featured by an increase in SSH along the equatorial Pacific and a decrease of SSH over most of the off-equatorial Pacific belt at about 10°N.The RD mechanism may facilitate the transition from El Niño phase to La Niña phase, and vice versa from La Niña to El Niño, forming a more regular ENSO cycle (Yu and Fang 2018).

Model evaluation
Before investigating the changes of El Niño characteristics, it is essential to evaluate the performance of various CMIP6 models in simulating the anomalous SST and CRE patterns during the mature phase (December, January and February) of El Niño in observations.We add CRE for evaluation here in addition to SST because the former has been recognized as an important source of uncertainty in simulating the interannual climate variability in the earlier CMIP5 models (Ying and Huang 2016;Li et al. 2018Li et al. , 2020)).Because some climate models participating in the CMIP6 program do not provide a complete set of radiation flux variables for CRE calculation, only 30 CMIP6 models are evaluated.
Figure 2 shows the observed spatial pattern of SST anomaly (SSTA) composited over the mature phase of 11 historical El Niño events (1979/80, 1982/83, 1986/87, 1991/92, 1994/95, 1997/98, 2002/03, 2004/05, 2006 (1979/80, 1982/83, 1986/87, 1991/92, 1994/95, 1997/98, 2002/ regions of both hemispheres.The Taylor diagram (Fig. 2b) shows that while there is a significant difference in representing the spatial variability of SSTA (as revealed by a wide range of the normalized standard deviation value from 0.5 to 1.5), CMIP6 models are generally skillful in simulating the observed SSTA pattern with correlation coefficients all greater than 0.75.The ensemble mean shows a high correlation of 0.87.We also notice that most CMIP6 models (21 out of 30) tend to underestimate the spatial variability of SSTA (i.e., the normalized standard deviation value < 1), implying a too-weak SSTA contrast during the mature phase of El Niño in CMIP6 simulations compared to observations.Some recent studies (Li et al. 2018;Wang et al. 2021;Li et al. 2022a, b) showed that biases in CRE simulation commonly seen in the state-of-the-art climate models may introduce anomalous low-level outflows to weaken the easterly trade winds over the tropical Pacific, in turn, reducing the upper-ocean mixing and leading to warm SST and positive SSH biases in the central Pacific to degrade the performance of CP El Niño simulations.Li et al. (2022a, b) further showed that this problem remains in many climate models contributing to the CMIP6 program due to the inappropriate treatment of the frozen hydrometeor radiative properties.Since surface wind stress, SST and SSH are important variables for extracting the key air-sea feedback mechanisms, evaluating the CMIP6 models' performance in simulating the CRE anomaly pattern associated with El Niño is necessary.Because the CERES/EBAF satellite data are available since March 2000, the simulated CRE change by 30 CMIP6 models at the surface level during the mature phase of 2 El Niño events (2002/03 and 2004/05) is evaluated against the CERES/EBAF data during their commonly available period (from Mar.2000 to Dec. 2005).It is worth noting that the abovementioned 2 El Niño events belong to CP El Niño.Therefore, the Taylor diagram shown in Fig. 3 simply measures the CMIP6 models' performance in reproducing the CRE anomaly pattern associated with CP El Niño, although conclusions drawn can be plausibly applied to EP El Niño.Surprisingly, the spatial pattern of CRE anomaly at the Earth's surface (Fig. 3a) is very different from the SSTA pattern typical of CP El Niño (Fig. 8a), showing a strong cooling over the central Pacific, with the peak amplitude stronger than − 20 watt m −2 .This cooling pattern tends to expand southeastward to impact the subtropical South Pacific.By contrast, the western Pacific is covered by a widespread CRE warming, with the peak amplitude slightly greater than 10 watt m −2 .Since cloud cover change is the major factor controlling surface CRE anomaly, it implies that clouds are enhanced mainly over the central and southern Pacific during the mature phase of El Niño rather than over the equatorial central and eastern Pacific where SSTA is maximum.To elaborate, the change in cloud area fraction composited over the mature phase of the 2002/03 and 2004/05 El Niño events from the CERES/ EBAF satellite data is plotted.As shown in Fig. 4, cloud area fraction significantly enhances over the south-central Pacific and decreases over the western Pacific, generally consistent with the CRE change shown in Fig. 3a.
Due to a poor simulation of cloud cover change, most CMIP6 models fail to reproduce the observed CRE anomaly associated with CP El Niño, with more than half (16 out of 30) of the members even producing a negative correlation (Fig. 3b).Only two models (ACCESS-CM2 and CMCC-CM2-SR5) barely pass the qualifying criteria, with a correlation coefficient slightly greater than 0.5 and a normalized standard deviation value close to 1.The failure in reproducing the CRE change associated with CP El Niño implies that there is much room for the CMIP6 models to improve the simulation of clouds in response to climate forcing.

Changes in El Niño characteristics
In this section, we investigate the changes of El Niño characteristics, including intensity, frequency and CP/EP El Niño ratio, in CMIP6 simulations with the 1pctCO2 experiment.The change is defined as the difference between the first and the second 70-year periods, corresponding approximately to a doubling of CO2 concentration during the 70-year time span between the two periods.Because some CMIP6 models (16 out of 30) do not provide sea surface height (SSH) output, which is a key variable to identify the air-sea feedback mechanisms associated with El Niño diversity, only 14 models are analyzed and discussed in the following contents (see Table 1 for a summary of the 14 CMIP6 models selected).
Figure 5a shows the projected changes of El Niño intensity between the two periods.Although significant discrepancies exist in projecting the future El Niño intensity changes (e.g., 6 CMIP6 models project an enhanced trend and 8 CMIP6 models project a weakened trend), the ensemble mean suggests a slight enhancement of El Niño intensity by about 2% from the first to the second 70-year periods, which is quite small considering a doubling of CO2 concentration during the 70-year time span between the two periods.Moreover, the model diversity appears to enlarge from the first to the second 70-year periods, implying an increasing uncertainty of simulation results.Specifically, the error bars show a slight enhancement of model diversity in simulating the El Niño intensity change from the first (1.264± 0.143) to the second 70-year periods (1.289 ± 0.164) and a relatively larger enhancement of model diversity in simulating the El Niño frequency change from the first (12.28± 3.02) to the second 70-year periods (12.78 ± 4.69).The above model diversity could be reduced by separating the El Niño events into CP and EP types as well as by considering only the models that can correctly produce the intensity (variance) of major air-sea feedback mechanisms, which will be discussed later.
Since the locations of El Niño SSTA centers can be different across various CMIP6 models, it is speculated that such a diversity in El Niño intensity change can be attributed to the opposite response to global warming between CP and EP El Niño.To elaborate, we separate the model El Niño events into EP and CP types (Table 2).Although discrepancies in projecting the intensity changes remain, the ensemble mean results suggest that the intensity of EP El Niño decreases by nearly 4.6% (from 1.30 to 1.24) while the intensity of CP El Niño increases by about 4.5% (from 1.12 to 1.17) from the first to the second 70-years periods.Because CP El Niño Figure 5b shows the projected change of El Niño frequency between the two 70-year periods.Although discrepancies in projecting El Niño frequency changes remain (e.g., 7 CMIP6 models project an increasing trend and 7 CMIP6 models project a decreasing trend), the ensemble mean suggests that the occurrence frequency of El Niño may increase by about 4% from the first to the second 70-year periods, which is small considering a doubling of CO2 concentration between the two periods.Again, the reason for such a modest change in El Niño frequency is caused by the opposite response to global warming between EP and CP El Niño.As shown in Table 3, the occurrence frequency of EP El Niño is projected to decrease by about 1.4% (from 4.92 to 4.85) while the CP El Niño is projected to increase by about 2% (from 7.54 to 7.69), thereby leading to a modest increase of the total El Niño frequency shown in Fig. 5b.Also, most CMIP6 models simulate more CP El Niño than EP El Niño during the entire 140-year period, again concurring with observations.To check whether the inter-model relationship supports the results of multi-model ensemble mean, we calculate the inter-model correlation coefficient between differences in SF variance and differences in CP/EP ratio for the 12 CMIP6 models (excluding AWI-CM-1-1-MR and NorCPM1 due to failure in detecting CP or EP El Niño).Surprisingly, the correlation between them is very low (r = − 0.09), showing little connection between SF variance and CP/EP ratio.However, when we consider only the models that can produce a strong enough variance of the SF mode (> 7%, as by results in Fig. 1), which include ACCESS-CM2, CESM2, CESM2-WACCM, CMCC-ESM2, MIROC6, NorESM2-LM and NorESM2-MM, the correlation coefficient largely enhances from − 0.09 to 0.61.Likewise, when we consider the intermodel correlation coefficient between differences in BF variance and differences in CP (or EP) El Niño intensity for the models that can produce a strong enough variance of the BF mode (> 15%, as implied by results in Fig. 1), which include ACCESS-CM2, CanESM5, CESM2, CESM2-WACCM, CMCC-ESM2, MIROC6, MRI-ESM2-0 and NorESM2-LM, the correlation coefficient associated with the intensity of CP El Niño is notably positive, with a value of 0.51, and the correlation coefficient associated with the intensity of EP El Niño is also positive, but with a significantly lower value of 0.36.The above findings clearly suggest that, as long as a climate model can correctly produce the intensity (variance) of major air-sea feedback mechanisms, the relationship between SF variance and CP/EP ratio is evident.The relationship between BF variance and El Niño intensity also exists, although BF tends to impose a stronger effect on the intensity of CP El Niño than the intensity of EP El Niño, generally consistent with the results of multi-model ensemble mean shown in Table 2.
In summary, CP and EP El Niño events respond to global warming in a very different manner, i.e., CP El Niño becomes stronger and more frequent but EP El Niño becomes weaker and less frequent in a warmer atmosphere.This opposite response accordingly leads to a modest change of El Niño intensity and frequency shown in Fig. 5 when these two types of El Niño were not separated.It is worth noting that the above findings are consistent with some previous studies that virtually no consensus could be achieved on the enhancement or weakening of El Niño under global warming (e.g., Zelle et al. 2005;Merryfield 2006;Yeh and Kirtman 2007;Latif and Keenlyside 2009;Collins et al. 2010;DiNezio et al. 2012;Stevenson 2012;Power et al. 2013;Kim et al. 2014;Cai et al. 2015).

Changes in air-sea feedback mechanisms
Because BF controls the intensity of El Niño and SF modulates the frequency of El Niño, especially the CP El Niño, we investigate how their explained variances might change under global warming.The discussion of RD is omitted here because its variance and projected change are the smallest among the three feedback mechanisms and, more importantly, RD controls the regularity from El Niño to La Niña rather than the intensity and frequency of El Niño (Yu and Fang 2018).
Figure 6 compares the differences in explained variance with the BF, SF and RD modes between the two 70-year periods.As in observations, most CMIP6 models suggest that BF is the leading mode associated with El Niño, followed sequentially by SF and RD modes, except for a few exceptions such as CanESM5 (BJ, RD, SF) and NorCPM1 (BJ, RD, SF) in the first 70-year period and NorCPM1 (BJ, RD, SF) in the second 70-year period that SF and RD modes switch their order.For the BF mode (Fig. 6a), most (9 out of 14) CMIP6 models project an upward trend of its variance under global warming and the ensemblemean result suggests an increase rate of nearly 3% from the first to the second 70-year periods, concurring with a slight enhancement of El Niño intensity from the ensemble-mean point of view shown in Fig. 5a.One example for supporting this argument is the AWI-CM-1-1-MR model (marked by *) that projects the largest increase of BF variance from 10.5 to 15% and the greatest enhancement of El Niño intensity by about 22.58% (see Figs. 5a  and 6a).
For the SF mode (Fig. 6b), the majority (12 out of 14) of CMIP6 models project an upward trend of its variance under global warming and the ensemble-mean results suggests a notable increase rate over 13%.Since SF has been recognized as an important mechanism triggering a CP El Niño, it implies that there will be more CP El Niño under global warming.Table 3 compares the differences in the numbers of EP and CP El Niño events along with the CP/EP ratios between the first and the second 70-year periods.Given the diversity in simulating the numbers of EP and CP El Niño events across various CMIP6 models, the ensemble mean suggests that the CP/EP ratio may increase by about 10% (from 1.53 to 1.68) from the first to the second 70-year periods, which is quite consistent with the enhanced role played by the SF mode under global warming shown in Fig. 6b.

Concluding remarks
In the present study, we investigate the potential changes of El Niño characteristics (including intensity, frequency and CP/EP El Niño ratio) under progressive global warming using 140-year CMIP6 simulation outputs with the Comparing the simulation differences between the first and the second 70-year periods exhibits significant discrepancies in projecting the future changes of El Niño characteristics, including intensity and frequency, under global warming.Given the diversity of CMIP6 projections, nonetheless, the ensemble mean suggests a slight enhancement of El Niño intensity (about 2%) and a modest increase of El Niño frequency (about 4%) from the first to the second periods.Notably, these changes are quite small considering approximately a doubling of the CO2 concentration during the 70-year time span between the two 70-year periods.
Our study further shows that the small changes in intensity and frequency should come from the opposite response to global warming between CP and EP El Niño.After dividing the simulated El Niño events into CP and EP types, we note that the intensity of EP El Niño weakens by nearly 4.6% (from 1.3 to 1.24) while the intensity of CP El Niño enhances by about 4.5% (from 1.12 to 1.17) from the first to the second 70-year periods.Because CP El Niño occurs more frequently than EP El Niño in CMIP6 simulations, this results in a slight enhancement (about 2%) of the total El Niño intensity.A similar situation occurs in projecting the future changes of El Niño frequency, i.e., the frequency of EP El Niño decreases by about 1.4% (from 4.92 to 4.85) but the frequency of CP El Niño increases by about 2% (from 7.54 to 7.69) from the first to the second 70-year periods, thereby only leading to a modest increase (about 4%) of the El Niño frequency.
To understand the physical causes behind the opposite response to global warming between CP and EP El Niño, the projected changes of BF and SF variances are examined.From the ensemble mean point of view, the percentage variance explained by the BF mode increases by about 3% from the first to the second 70-year periods.Since BF controls the intensity of El Niño, the slight enhancement of El Niño intensity under global warming seems to be related to the increased BF variance.Moreover, the CMIP6 models project a notable increase (about 13%) of the variance explained by the SF mode from the first to the second 70-year periods.Since SF is known for triggering a CP El Niño, the increased CP/EP ratio (about 10%) in CMIP6 ensemble-mean simulations can be explained by the notably enhanced role played by the SF mechanism in a warmer atmosphere.Our study also points out that, as long as a climate model can correctly produce the intensity (variance) of major air-sea feedback mechanisms, the relationship between changes in El Nino characteristic

Fig. 3
Fig. 3 Same as Fig. 2, but showing the performance of CRE anomaly simulations during the mature phase of 2 historical El Niño events (2002/03 and 2004/05) by individual CMIP6 models against the CERES/EBAF data.Only two El Niño events match with the CERES satellite service period (2001-2005).a The composited CRE anomaly pattern and b the Taylor diagram measuring the performance of individual CMIP6 models against the CERES/EBAF data

Fig. 4
Fig. 4 Same as Fig. 3a, but showing the spatial distribution pattern of cloud area fraction from CERES/EBAF data

Fig. 5
Fig. 5 Changes of El Niño (a) intensity (measured by the normalized Niño3.4 index) and (b) frequency between the first and the second 70-year periods in CMIP6 simulations.The percentage numbers in red (blue) indicate that the El Niño intensity, or frequency, increase (decrease) from the first to the second 70-year periods.In each diagram, the left and right histograms represent the ensemble mean values for the first and the second 70-year periods, respectively.The error bars, with the length denoting plus-and-minus one standard deviation, measure the diversity of CMIP6 models

Fig. 6
Fig. 6 The explained variances (in %) associated with (a) BF, (b) SF and (c) RD modes in various CMIP6 models, as well as the multi-model mean (MME), during the first (blue bar) and the second (red bar) 70-year periods.The ensemble-mean change rate between the two periods (laterformer) is displayed at the top-left corner of each panel

Fig. 7
Fig. 7 The correlation maps of the principal component time series (PC) derived from the SF mode of MEOF with the SSTA (shadings) and surface wind (vectors) anomalies over the tropical Pacific (20°S-30°N, 120°E-80°W) from lag 0 to lag 10 months derived from the 140-year CMIP6 ensemble mean simulations.Only wind vectors significant at the 90% level are displayed.The reference scale for winds at the top-right corner is 0.25 m s −1

Table 1
A list of 14 CMIP6 models used in this study

Table 2
The normalized intensity of EP and CP El Niño events Failure in satisfying the criteria of EP and/or CP El Niño is denoted by "Null".Numbers with boldface font indicate the ensemble mean values

Table 3
The number of EP and CP El Niño and the associated CP/EP ratio, respectively, in the first and second (in parentheses) 70-year period simulated by CMIP6 models Failure in satisfying the criteria of EP and/or CP El Niño is denoted by "Null".Numbers with boldface font denote the ensemble mean values