Sensitivity of ENSO Simulation to the Convection Schemes in the NESM3 Climate System Model: Atmospheric Processes

The El Niño-Southern Oscillation (ENSO) is the most prominent climate system in the tropical Pacific. However, its simulation, including the amplitude, phase locking, and asymmetry of its two phases, is not well reproduced by the current climate system models. In this study, the sensitivity of the ENSO simulation to the convection schemes is discussed using the Nanjing University of Information Science and Technology Earth System version 3.0 (NESM3) model. Three convection schemes, including the default, the default coupled with the stochastic multicloud model (SMCM), and the default used in the Coupled Model Intercomparison Project Phase 6 (CMIP6), are implemented. The model results reveal that the low-level cloud cover and surface net shortwave radiation are best represented over the tropical Pacific in the model containing the SMCM. The simulations of the ENSO behavior’s response to changes in the convection scheme are not uniform. The model results reveal that the model containing the SMCM performs best in terms of simulating the seasonal cycle of the sea surface temperature anomaly along the equatorial Pacific, the phase locking, and the power spectrum of ENSO but with a modest ENSO amplitude. Compared to the model containing the default convection scheme, the coupling of the default scheme and the SMCM provides a good simulation of the ENSO’s asymmetry, while the model containing the CMIP6 convection scheme outperforms the others in terms of the simulation of the ENSO’s amplitude. Two atmospheric feedback processes were further discussed to investigate the factors controlling the ENSO’s amplitude. The analyses revealed that the strongest positive atmospheric Bjerknes feedback and the thermodynamic damping of the surface net heat flux occurred in the model containing the CMIP6 convection scheme, suggesting that the atmospheric Bjerknes feedback may overwhelm the heat flux damping feedback when determining the ENSO’s amplitude. The results of this study demonstrate that perfectly modeling and predicting the ENSO is not simple, and it is still a large challenge and issue for the entire model community in the future.


INTRODUCTION
The El Niño-Southern Oscillation (ENSO), which has been studied extensively over the last few decades, plays a central role in modulating the tropical and global climate and weather via its strong interannual variability. Moreover, the warm and cold phases of the ENSO have diverse impacts in different regions; that is, some regions suffer severe flooding while other suffer severe droughts (e.g., Chiodi and Harrison, 2015). Moreover, the ENSO also exerts tremendous influences on the activities of typhoons (e.g., Wang et al., 2013), the Madden-Julian Oscillation (MJO) (e.g., Wei and Ren, 2019), the Indian Ocean Dipole (e.g., Hong et al., 2008;Fan et al., 2017), the monsoon system (e.g., Torrence and Webster, 1999), and other climatic phenomena. Thus, accurately representing and predicting the ENSO using climate system models (CSMs) is crucial for socioeconomic development.
Simulating and predicting the ENSO is a significant challenge due to its complicated physics, which involve atmospheric and oceanic processes and their interactions. Although significant progress has been made in developing CSMs by improving the physical parameterizations and increasing the horizontal and vertical resolutions during the last few decades, common problems related to ENSO simulations, such as the spurious phase locking and excessively regular variations, are still encountered. As the results of the fifth phase of Coupled Model Intercomparison Project (CMIP5) reported, there is a large spread of ENSO simulations with respect to individual models, but there were no significant changes compared to the third phase of the CMIP (CMIP3) in terms of the multimodel mean state (e.g., Guilyardi et al., 2012). The simulated cloud radiation feedback during the ENSO cycle exhibits great uncertainty in coupled models (Chen et al., 2013). In the majority of the CMIP3 and CMIP5 models, the response of the shortwave radiation anomaly to the ENSO cycle and the two types of El Nino is severely underestimated (Chen et al., 2018;Chen et al., 2019a, b), and it exhibits a westward-shifted bias (Ge and Chen, 2020). Moreover, the model results from the sixth phase of the CMIP (CMIP6) still exhibit diversity with respect to the phase locking of the ENSO (McKenna et al., 2020).
Previous studies have addressed the fact that the atmospheric processes play vital roles in shaping the ENSO's properties (e.g., Guiyardi et al., 2004;Kim et al., 2008;Watanabe et al., 2010;Lloyd et al., 2011). Based on a coupled model, Lengaigne et al. (2006) concluded that the peak and termination of intense El Niño events are tightly linked to the seasonal evolution of solar insolation. In addition, ENSO-associated circulation can be largely altered by convective momentum transport (Kim et al., 2008). Moreover, the convective schemes in the atmospheric component of CSMs have been considered to be a deterministic factor in simulating the ENSO's behaviors (Wu et al., 2007;Neale et al., 2008;Lloyd et al., 2012). For example, when the deep convection is represented by the Kerry-Emanuel scheme in the coupled L'Institut Pierre-Simon Laplace (IPSL) model, the amplitude of the ENSO has been reproduced .
The Nanjing University of Information Science and Technology (NUIST) Earth System Model (NESM) was developed by the Earth System Model Center of NUIST. Since version 1.0 of the NESM (NESM1) was launched (Cao et al., 2015), it and its descendants have been successfully applied to the study of monsoon dynamics (e.g., Li et al., 2017;Li and Wang, 2018;Cao et al., 2019a and atmosphere-sea iceocean interactions in the Southern Ocean (Ma et al., 2020b). Moreover, version 3 of the NESM (NESM3) has been registered with the CMIP6 (Cao et al., 2018). The associated data produced by the NESM3 have been submitted to the CMIP6 and can be downloaded. This study incorporates three convective schemes into the atmosphere component of the NESM3 in order to study their effects on ENSO simulations.
The rest of this paper is arranged as follows.

NESM3 and Experimental Design
The NESM3 was deliberately constructed and developed to participate in the CMIP6. It consists of the ECHAM6.3.02 atmospheric component (Stevens et al., 2013), the NEMO3.4 ocean component (Madec, 2012), and the CICE4.1 sea ice component (Hunke and Lipscomb 2010). They are coupled together using the OASIS-MCT3.0 coupler (Valcke and Coquart, 2015). Additionally, the JSBACH land model (Raddatz et al., 2007) is included in ECHAM6.3.02. The standard resolution (SR) version of the NESM3, in which the horizontal resolutions of ECHAM56.3.2 (with 47 vertical layers) and JSBACH (with 5 vertical layer) are T63 spectral truncation (approximately 1.875°) and a nominal 1°is applied to NEMO3.4 (with 46 vertical layers) and CICE4.1, is used in the study. For detailed information related to the NESM3 and its improvements, refer to Cao et al. (2018).
Three convection schemes are implemented in the NESM3. First is the default convection scheme in ECHAM6.3.02 (referred to as NESM_CTRL); the second adopts the modifications of the convection scheme made by Yang and Wang (2018) (referred to as NESM_CMIP6); and the third couples a stochastic multicloud model (SMCM) with the default convection scheme in ECHAM6.3.02 (referred to as NESM_SMCM). The modified convection scheme used in NESM_CMIP6 has been demonstrated to significantly improve MJO simulations , and it is already used in the NESM3 in the CMIP6. Further information about the modified convection scheme in NESM_CMIP6 can be found in Yang and Wang (2018) and Yang et al. (2020). In addition, the incorporation of the SMCM into ECHAM6.3.02 improves MJO simulations (Peters et al., 2017;Ma et al., 2019a) and monsoon simulations Ma et al., 2020a;Ma and Jiang, 2020). Detailed discussions of the SMCM can be found in Peters et al. (2017). Additionally, brief descriptions of the three convection schemes are included in Supplementary Material Text S1. Moreover, no other modifications were made except to the convection scheme in the atmospheric component of the NESM3. For each convection scheme, one experiment was conducted, and 100year integrations were run with the external forcings for the 1990s. The last 50 years of monthly data were used for the analysis in this study.

Datasets and Methodology
Monthly sea surface temperature (SST) data from the National Oceanic and Atmospheric Administration Extended Reconstructed SST (NOAA ERSST) v4.0 ranging from 1956 to 2005 were used in this study (Huang et al., 2005). Monthly cloud cover data from 1984 to 2008 were obtained from the International Satellite Cloud Climatology Project (ISCCP) (Rossow et al., 1996). Monthly cloud radiation forcing data from 1984 to 2007 were obtained from the National Aeronautics and Space Administration Global Energy and Water Exchanges (NASA/GEWEX) surface radiation budget (SRB) project (Gupta et al., 2006). In addition, the 1984-2009 OAFlux monthly surface fluxes (Yu and Weller, 2007), the 1991-2010 zonal wind stress data derived from the National Center for Environmental Prediction and National Center for Atmospheric Research (NCEP NCAR) reanalysis (Kalnay et al., 1996), and the 1984-2009 ORAS4 (Ocean Reanalysis System 4) (Balmaseda et al., 2013) were used in this study. Moreover, five monthly SST datasets  from CMIP6 of the NESM3 historical simulation (i.e., r1i1p1f1, r2i1p1f1, r3i1p1f1, r4i1p1f1, and r5i1p1f1) were also used. It should be noted that, for convenience, the observation and reanalysis data are all referred to as observations in this study. Additionally, the modeled and observational atmospheric data were horizontally interpolated onto a 2.5°× 2.5°grid, whereas the ocean data were horizontally interpolated onto a 1°× 1°grid. The anomaly in each field was obtained by removing the long-term averaged seasonal cycle for each grid.
In order to explore how the convection scheme influences the modeled ENSO, two atmospheric feedback processes are discussed in addition to the basic evaluations (e.g., ENSO's amplitude, phase locking, and period). One is the atmospheric Bjerknes feedback, τ′ Niño4 µSSTA Niño3 , which is measured by the linear regression coefficient between the time series of the zonal wind stress anomaly, τ ' , averaged over the Niño-4 region (5°S-5°N, 160-210°E), and the time series of the SST anomaly (SSTA) averaged over the Niño-3 region (5°S-5°N, 210-270°E) (Lloyd et al., 2012;Bellenger et al., 2014). The atmospheric Bjerknes feedback calculated in this study represents the nonlocal response of the zonal wind stresses to a given equatorial Pacific SSTA, which is a positive feedback and also measures the coupling strength between the Niño-3 SSTA and the Niño-4 zonal wind stress anomaly (e.g., Zebiak and Cane, 1987;Chen et al., 2016;Chen et al., 2019c). The other is the so-called heat flux feedback, Q ' αSSTA, which is defined as the regression coefficient between the net surface heat flux anomaly and the SSTA in a specific region, and it measures the local negative feedback between the SSTA and the changes in the surface net heat flux (e.g., Jin et al., 2006). Two additional objective metrics, the pattern correlation coefficient (PCC) (e.g., Zhu and Li, 2017) and the normalized root mean square error (NRMSE) (e.g., Lee and Wang, 2014), were also calculated to gauge and quantify the model's performance. The PCC and NRMSE are calculated as follows: where covar(·, ·) and std(·) are the covariance and standard deviation, respectively, and X (x 1 , x 2 , /, x M ) and Y (y 1 , y 2 , /, y M ) are the observations and model simulation, respectively.

Responses of Cloud and Surface Radiation
There is no doubt that modifying the convection scheme results in changes in the cloud properties. Figure 1 shows the annual-mean low-level, middle-level, and high-level cloud cover in the observations and model simulations. With respect to the low-level cloud cover ( Figures 1A-D), overall, all of the simulations underestimate it over the tropical Pacific compared to the observations, especially over the western Pacific. It should be noted that the model simulations exhibit different biases. The center of the negative bias in NESM_CTRL is located in the western Pacific, while those of NESM_SMCM and NESM_CMIP6 are located in the central Pacific. Moreover, among all the simulation, NESM_CMIP6 most severely underestimates the low-level cloud cover. Based on the PCC and NRMSE calculated over the tropical Pacific (20°S-20°N, 100-280°E), NESM_SMCM has the largest PCC, which is 0.86 versus 0.66 for NESM_CTRL and 0.65 for NESM_CMIP6. Moreover, NESM_CTRL has the smallest NRMSE, which is 1.45 versus 2.02 for NESM_SMCM and 2.51 for NESM_CMIP6. Regarding the middle-level cloud cover, NESM_CTRL underestimates it over the tropical Pacific ( Figure 1F), while NESM_SMCM and NESM_CMIP6 overestimate it ( Figures 1G,H). A comparison of all of the simulations indicates that NESM_CMIP6 has the best pattern correlation coefficient (PCC 0.87) and NESM_SMCM has the smallest bias (NRMSE 1.33). In addition, similar distributions were found for the modeled high-level cloud cover, which are overestimated compared to the observations ( Figures 1I-L). Figure 2 shows the annual mean of the total cloud cover, surface net solar radiation, and longwave radiation. As in the observations, the minimum amount of total cloud cover occurred over the southern central Pacific (Figure 2A), while the minimum total cloud cover in all of the simulations occurred in the southeastern corner of the tropical Pacific (south of 10°S) ( Figures 2B-D). Additionally, three simulations overestimated the total cloud cover over the western Pacific, with the largest amount overestimation occurring for NESM_CMIP6. The PCC and NRMSE calculated over the tropical Pacific indicate that NESM_SMCM has slightly better performance; i.e., it has the largest pattern correlation (0.77) and the smallest bias (0.71). The discrepancies in the cloud cover no doubt lead to differences in the radiation distribution at the surface. Regarding the surface net solar radiation, first of all, the general spatial patterns were captured by all of the simulations but with amplitude and position biases. Compared to the observations ( Figure 2E), NESM_CTRL has a comparable strength but is located too far to the west ( Figure 2F), while NESM_SMCM and NESM_CMIP6 simulate the correct position but underestimate the strength ( Figures 2G,H). Moreover, the minimum solar radiation occurred over the western Pacific in NESM_CMIP6 due to the maximum amount of total cloud cover being located there, which potentially induces the weakest zonal gradient in the sea surface temperature between the equatorial eastern and western Pacific. Similar distributions were found for longwave radiation ( Figures  1I-L). The PCC and NRMSE demonstrate that NESM_SMCM has the best performance in terms of simulating the solar and longwave radiations, that is, 0.84 and 0.74 for solar radiation and 0.65 and 0.77 for longwave radiation, respectively.

Sensitivity of ENSO Simulations to Convection Schemes
The responses of the ENSO's basic features to the convection schemes are discussed in this section. Figure 3 shows the annual mean and standard deviation of the SST. In the observations, the 28°C isotherm of the equatorial SST extends westward to 170°W Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 596442 ( Figure 3A). Among the three simulations, ECHAM_CMIP6 captures the eastward extension of the 28°C isotherm ( Figure 3D), while the others mismodel this feature; that is, the 28°C isotherm of the equatorial SST extends eastward to 170°E in NESM_CTRL ( Figure 3B) and extends eastward to the date line in NESM_SMCM ( Figure 3C). Moreover, colder and narrower cold tongues occurred in all of the simulations, while warmer SSTs occurred along the coast of Peru, indicating that simulated trade winds are weaker than those observed. Moreover, a colder warm pool occurred in NESM_CTRL and NESM_CMIP6, whereas a positive SST bias in the warm pool occurred in NESM_SMCM. The coldest warm pool and the warmest SST in the equatorial eastern Pacific occurred in ECHAM_CMIP6, which is consistent with the distribution of the surface net solar radiation, implying that the weakest zonal gradient of the equatorial SST was modeled. As indicated by the standard deviation of the SST, all of the simulations underestimated the variability of the equatorial SST, suggesting that less ENSO events and weaker ENSO amplitudes were modeled compared to the observations ( Figures 3E-H). It should be noted that NESM_CMIP6 performed best in terms of simulating the ENSO's amplitude (PCC of 0.90), while NESM_CTRL had the smallest PCC (0.78). It has been pointed out that the seasonal cycle of the SST anomalies (SSTA) along the equator plays a vital role in ENSO evaluation (e.g., An and Choi, 2013). The observed seasonal cycle of the SSTA is reproduced by all of the simulations with a warm phase in the first half of the year and a cold phase in the second half of the year (Figures 4A-D). However, biases were also found in simulations, especially in the eastern equatorial Pacific. Larger positive and negative biases were found in NESM_CTRL and NESM_CMIP6, while they were alleviated in NESM_SMCM ( Figures 4E-G). The PCC and NRMSE in the longitudinal range of the Niño 3 indicate that NESM_SMCM performs best in terms of simulating the seasonal cycle of the equatorial SSTAs compared to Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 596442 6 NESM_CTRL and NESM_CMIP6, potentially indicating that the ENSO's evolution is best simulated by NESM_SMCM.
The observed variability of the ENSO in the Niño 3, which exhibits a strong phase locking of the seasonal cycle with a minimum during the March-May period and a maximum during the November-January period in terms of the standard deviation of the SSTA during each month, is not well resolved by the model simulations ( Figure 5A). The seasonal cycle of the ENSO's variability is completely reversed in NESM_CTRL and NESM_CMIP6, with a maximum SSTA standard deviation during the March-May period and a minimum during the November-January period. Moreover, the maximum and minimum SSTA standard deviations in NESM_SMCM occurred in May and November, respectively. A seasonality metric, which is defined as the ratio of the seasonal cycle of the SSTA standard deviation between the November-January and March-May periods, was calculated to further assess the seasonal amplitude and phase of the modeled ENSO (Bellenger . The seasonality metrics in NESM_CTRL, NESM_SMCM, and NESM_CMIP6 were 0.86, 1.08, and 0.86, respectively, which are smaller than those in the observations (1.57). Moreover, the temporal correlation of the seasonal cycles of the SSTA standard deviation between NESM_SMCM (NESM_CTRL, NESM_CMIP6) and the observations was 0.50 (−0.53, −0.73). The power spectrum of the time series of the Niño-3 SSTA was also calculated. The observed power spectrum exhibits a broad spectral range with four peaks, which are in the 95% confidence level, during the 2-7 year period (dashed line in Figures 5B-D). In general, the modeled Niño-3 SSTAs exhibit the observed broad spectral range during the 2-7 year period, but fewer peaks are in the 95% confidence level. For example, peaks occurred between the 3-5 year period and the 2-3 year period in the 95% confidence level in NESM_CTRL and NESM_SMCM, respectively ( Figures 5B,C). Moreover, two peaks occurred in the 95% confidence level during the 2-5 year period in NESM_CMIP6 ( Figure 5D). In order to evaluate the resemblance of the power spectra of the Niño-3 SSTA among the observations and model simulations, the correlation coefficients between the model simulations and the observations with respect to the spectral range of the 2-7 year period were computed, and the results were 0.26, 0.43, and 0.17 for NESM_CTRL, NESM_SMCM, and NESM_CMIP6, respectively. This indicates that the variation in the Niño-3 SSTA in NESM_SMCM most closely resembles the observations. The simulations of the ENSO's variability were better in the Niño-3.4 region (Supplementary Figure S1A). The temporal correlations between the modeled and observed seasonal cycle of the SSTA standard deviation were 0.60, 0.91, and 0.60 for NESM_CTRL, NESM_SMCM, and NESM_CMIP6, respectively. Moreover, the seasonality metrics in the observations, NESM_CTRL, NESM_SMCM, and NESM_CMIP6 were 1.79, 1.11, 1.35, and 1.07, respectively. It should be noted that NESM_SMCM also outperforms the others in simulating the ENSO's variability in the Niño-3.4 region. However, there were no significant changes with respect to the power spectra of the modeled time series of Niño-3.4 SSTA compared to those of the Niño-3 SSTA (Supplementary Figures S1B-D). Additionally, five members from historical experiments using the NESM3 submitted to the CMIP6 were also analyzed here. As shown in Supplementary Figure S2, differences in the ENSO's variability and period in terms of the time series of the Niño-3.4 SSTA were found, which deserves further investigation in future studies.
The ENSO's asymmetry is a salient feature, which reflects the unequivalence of the amplitudes between the two phases of the ENSO (e.g., Burgers and Stepenson, 1999;Zhang and Sun, 2014;Tang et al., 2019). Figure 6 shows the distribution histograms of the time series of the Niño-3 SSTA in the observation and model simulations. It is shown that the observed SSTA exhibits a long tail on the right, implying that the frequency (amplitude) of strong El Niño events is higher (larger) than that of strong La Niña events ( Figure 6A). This observed feature, i.e., a long tail on the right, is not captured by NESM_CTRL and NESM_CMIP6 ( Figures 6B,D), while it appears in NESM_SMCM ( Figure 6C). Moreover, the maxima of the simulated positive anomalies are smaller than that in the observations, suggesting that the simulated strong El Niño events are weaker than those in the observations. To further illustrate the asymmetry of the ENSO's Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 596442 amplitude, a metric, i.e., skewness, which measures the ENSO's nonlinearity, was calculated. The positive skewness indicates that the amplitude of the warm events is frequently larger than that of the cold events. However, the skewness in NESM_CTRL is −0.03, suggesting that the amplitude of the warm events is even weaker than that of the cold events. Moreover, the skewness in NESM_SMCM is 0.83, which is comparable to the observed skewness (0.81), while the skewness in NESM_CMIP6 is 0.29.

Atmospheric Feedback Processes
Theoretical and modeled analyses have pointed out that atmospheric processes play vital roles in shaping the activities of simulated ENSO events (e.g., Wang and Picaut, 2004;Watanabe et al., 2010;Llyod et al., 2012). As discussed above, the responses of the ENSO amplitude to the convection schemes are diverse, which is potentially dominated by the intermodel differences in the thermodynamic damping of the SSTs (e.g., Lloyd et al., 2009;Lloyd et al., 2012). Thus, two associated atmospheric feedback processes, the atmospheric Bjerknes and heat flux feedbacks, are discussed to investigate the influences of the convection schemes used in this study on the modeled ENSO. Figure 7 shows the atmospheric Bjerknes feedback µ in the observations and model simulations. Compared to the observed µ (9.76 × 10 -3 N m −2 /°C), a comparable µ was found in NESM_SMCM (9.57 × 10 -3 N m −2 /°C), while NEMS_CTRL underestimates (8.43 × 10 -3 N m −2 /°C) the coupling strength between the Niño-4 τ ' and the Niño-3 SSTA and NESM_CMIP6 overestimates it (11.88 × 10 -3 N m −2 /°C). Recall that the standard deviations of the time series of the SSTA in the Niño 3 were 0.86, 0.67, 0.65, and 0.73 in the observations, NESM_CTRL, NESM_SMCM, and NESM_CMIP6, respectively ( Figure 6). This indicates that the remote coupling strength between τ ' in the Niño 4 and the SSTA in the Niño 3 may be a major factor in shaping and modulating the ENSO's amplitude. Figure 8 shows the linear pointwise regression coefficients between the surface net heat flux anomaly and the SSTA. All of the simulations have the capability to reproduce the observed negative heat flux feedback over the equatorial Pacific; i.e., an intense α is located over the western Pacific, and a weak α is located over the central-eastern equator. The area-averaged net heat flux α over the Niño 3 was −16.45, −8.36, −10.71, and −12.79 W m −2 /°C in the observations, NESM_CTRL, NESM_SMCM, and NESM_CMIP6, respectively. It is desired that the larger thermodynamic damping leads to a weaker ENSO amplitude. This relationship does not hold true in NESM_CMIP6, potentially suggesting that, in addition to the atmospheric Bjerknes feedback, ocean processes also have great influence on the ENSO's amplitude. The surface net heat flux feedback α is composed of four components: shortwave radiation, longwave radiation, the latent heat flux, and the sensible heat flux. Figure 9 presents scatterplots of the shortwave radiation (α sw ) and latent heat flux (α lh ) as a function of the SSTA in the Niño 3. The negative thermodynamic damping of the shortwave radiation is reproduced in NESM_SMCM and NESM_CMIP6, while it is mismodeled in NESM_CTRL ( Figures 9A-D). However, the modeled α sw in NESM_SMCM and NESM_CMIP6 are underestimated. With respect to α lh , the observations and all of the simulations exhibit negative damping effects on the Niño-3 SSTAs ( Figures 9I-L). Moreover, all of the simulations underestimate the damping effect of the latent heat flux. In addition, as shown in Figures 9E-H,M-P, the feedbacks of the longwave radiation and the sensible heat flux play less significant roles in the surface net heat flux feedback, which indicates that the surface net heat flux feedback is dominated by the shortwave radiation and the latent heat flux feedbacks in the Niño-3 region (e.g., Lloyd et al., 2009).
The scatterplot of the observed shortwave radiation versus the SSTA in the Niño 3 ( Figure 9A) illustrates their nonlinear relationship; that is, the shortwave feedback is positive for cold SSTAs and negative for positive SSTAs (e.g., Zebiak and Cane, 1987;Barnet et al., 1991;Lloyd et al., 2012;Bellenger et al., 2014). Thus, the linear regression coefficient of the shortwave radiation anomaly was separately calculated for the negative SSTA (α sw− ) and for the positive SSTA (α sw+ ), which were used to evaluate the nonlinearity of the shortwave radiation feedback (α sw− − α sw+ ). As for α sw− , NESM_SMCM mismodeled it, with a negative feedback of the shortwave radiation with respect to the negative SSTAs in the Niño-3 region. In addition, NESM_CTRL failed to reproduce the negative feedback of the shortwave radiation with respect to the positive SSTAs in the Niño-3 region. It should be pointed out that NESM_CMIP6 produces the right shortwave feedback for both the negative/positive SSTAs in the Niño-3 region. Thus, the nonlinearity α sw− − α sw+ in NESM_CTRL, NESM_SMCM, and NESM_CMIP6 was 0.89, 3.14, and 8.45 W m −2 /°C, respectively, which are smaller than that of the observations (13.76 W m −2 /°C). It has been reported that the shortwave radiation feedback α sw is influenced by the modeled cloud biases, which can be highlighted by the cloud radiation forcing (CRF) (Lloyd et al., 2011(Lloyd et al., , 2012Chen et al., 2019). Because large cloud biases were found over the tropics (Figures 1, 2), the CRF is discussed here to explore the model biases for the shortwave radiation feedback. Consistently, the regression between SSTA and CRF was calculated for the positive and negative Niño-3 SSTAs separately. Figure 10 shows the shortwave CRF response to the SSTA changes with respect to the positive and negative Niño-3 SSTA. In the observations, a prominent anomalous negative shortwave CRF was found, and it was centered to the east of the date line ( Figure 10A). In NESM_CTRL, the pronounced negative shortwave CRF anomalies were weakened and shifted westward compared to the observations ( Figure 10B). A similar distribution was found in NESM_SMCM but with stronger shortwave CRF anomalies compared to NESM_CTRL ( Figure 10C). In contrast, the negative shortwave CRF anomalies in NESM_CMIP6 were comparable to the observations, but with enhanced positive shortwave CRF anomalies over the tropical western Pacific ( Figure 10D). In addition, the PCC and NRMSE scores between the observations and simulations were calculated over the tropical Pacific region (20S-20oN, 160-210oE). The results show that the PCC between NESM_CMIP6 and the observations (0.57) was larger than that between NESM_CTRL and the observational result (0.17) and between NESM_SMCM and the observational results (0.50), while NESM_CMIP6 has the largest NRMSE (1.73). With respect to the negative Niño-3 SSTA, two prominent negative shortwave CRF anomalies were located over both sides of the central equatorial Pacific in the observations ( Figure 10E). The modeled biases vary in the different simulations. Compared to the observations, the negative shortwave CRF anomalies were seriously underestimated in NESM_CTRL ( Figure 10F), while the biases were alleviated in NESM_SMCM ( Figure 10G) and NESM_CMIP6 ( Figure 10H). The PCC score also indicates that NESM_SMCM performs best in simulating the spatial pattern of the observed shortwave CRF under the conditions of a negative Niño-3 SSTA. The performance of the shortwave CRF response to the Niño-3 SSTA is consistent with the nonlinearity of the shortwave feedback α sw− − α sw+ . This implies that the modeled CRF has the potential to modulate the performance of the simulated shortwave feedback α sw (e.g., Lloyd et al., 2011).

DISCUSSIONS AND CONCLUSION
The influence of the convection schemes on the ENSO simulations in the NESM3 was explored in this study. Three convection schemes were implemented and investigated using the NESM3. First, the model results revealed that the changes in the convection schemes significantly redistribute the cloud cover and the surface heat flux, including the solar radiation. Modification with the SMCM improves the spatial pattern of the low-level Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 596442 cloud cover over the tropical Pacific, while NESM_CMIP6 improves the spatial pattern of the middle-level cloud cover compared to the default convection scheme used in the NESM3. Improvements in simulating the surface net solar radiation were obtained for NESM_SMCM and NESM_CMIP6, and improvements in simulating the longwave radiation were obtained for NESM_SMCM. Furthermore, the solar radiation was underestimated over the western Pacific in NESM_CTRL due to the intense total cloud cover there.
The annual mean and standard deviation of the SST were also altered. Specifically, the westward extension of the 28°C isotherm in NESM_CMIP6 was comparable to the observations, while the 28°C isotherm extended too far westward in NESM_CTRL and NESM_SMCM. In addition, both modified convection schemes produced better SST standard deviation spatial patterns over the tropical Pacific, while the weakest variability in the equatorial Pacific occurred in NESM_SMCM. As for the annual cycle of the SSTA along the equator, the bias was reduced in the eastern Pacific in NESM_SMCM, leading to a better performance, whereas it was slightly worse in NESM_CMIP6.
Consistent with the diverse performances in terms of the mean state of the simulated SST, the characteristics of the modeled ENSOs were also different for the different convection schemes. When the default scheme was coupled with the SMCM, the seasonality and phase locking were reproduced the best, while the seasonal cycle of the modeled ENSO variability was totally reversed; that is, the maximum SSTA standard deviation occurred in March-May and the minimum occurred in November-January, in NESM_CTRL and NESM_SMCM. However, the simulated ENSO was weakest in NESM_SMCM. Additionally, the ENSO nonlinearity measured by the skewness of the SSTA in the Niño 3 and the observed broad spectral range were simulated best in NESM_SMCM.
The analyses related to the atmospheric feedback processes revealed that the factors responsible for the ENSO's amplitude are complex. Strongest positive atmospheric Bjerknes feedback was found in NESM_CMIP6, which potentially enhances the ENSO's amplitude compared to the other two simulations. Regarding the negative surface heat flux feedback, the strongest damping in NESM_CMIP6 was desired to form the weakest modeled ENSO, but this does not hold true. This indicates that the atmospheric Bjerknes feedback may dominate the thermodynamic damping feedback and other processes, such as ocean processes, i.e., the zonal transport of the seawater temperature from the equatorial western Pacific, and it may also play central roles in determining the simulated ENSO's amplitude (e.g., Chen et al., 2019). Thus, further analyses related to ocean processes based on the upperocean heat budget and stability analysis (i.e., the BJ index) will be studied in the future. Moreover, the nonlinear relationship between the shortwave feedback and the SSTA in the Niño 3 was captured in NESM_CMIP6, while the negative shortwave feedback was misrepresented in the model containing the default convection scheme.
This study emphasizes that as a climate system with strong atmosphere-ocean interactions, ENSO's characteristics are hard to reasonably simulate by modifying only one aspect. Other physical aspects also need to be taken into consideration. For example, the resolutions of the ocean and atmospheric components in a CSM are considered to be an important factor for improving ENSO simulations (e.g., Schneider et al., 2003;Gualdi et al., 2005). In addition, ENSO simulations are also tightly linked to the ocean mixing (e.g., Richards et al., 2012;Qiao et al., 2013), the parameterization of the subsurface entrainment (e.g., Zhu et al., 2013), and so on. Thus, fully understanding the ENSO's physics is of significance for improving the model's ability to simulate ENSO and for ENSO prediction. However, when modifying the physical processes associated with ENSO, one should be careful and should consider the interactions and feedbacks with other processes, avoiding the occurrence of the "error-cancel-error" phenomena, which leads to improvements in some aspects of the modeled ENSO. It should be noted that NESM_CMIP6 had the strongest atmospheric Bjerknes feedback and thermodynamic damping. However, NESM_CMIP6 overestimates the two atmospheric feedback processes, which make opposite contributions to the ENSO's amplitude, compared to the observations. This raises questions as to whether the largest ENSO amplitude found in NESM_CMIP6 is a processcompensating product and which physical processes control and dominate the ENSO's amplitude. These issues deserve further study and are important for improving the models' performance during future model development.

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

AUTHOR CONTRIBUTIONS
LM sets and runs the numerical experiments. LM and ZJ write the manuscript and explain the results.