1 INTRODUCTION

El Niño is referred to anomalous warming of sea surface temperature (SST) in the eastern-central equatorial Pacific Ocean, and the Southern Oscillation (SO) is referred to the differences in sea level pressure (SLP) across the tropical Pacific Ocean. Their combination coins the El Niño-Southern Oscillation (ENSO), which occurs irregularly every 2–7 years (Chao, 1993). As the largest interannual signal originating from the coupled ocean-atmosphere interactions in the tropical Pacific, ENSO significantly affects weather and climate all over the world.

In the past several decades, ENSO has been extensively and intensively investigated; significant progress has been made in ENSO studies (e.g., McCreary and Anderson, 1991; Neelin et al., 1992; Philander, 1999; McPhaden et al., 2006). Before the 1960s, the El Niño and SO phenomena were separately analyzed as two independent signals in the ocean and atmosphere, respectively; some basic characteristics of interannual anomalies for SST and atmospheric surface fields (e.g., SLP) were identified based on limited observations. In the late 1960s, it has been recognized that the El Niño and SO phenomena are two closely related aspects of the same large-scale ocean-atmosphere interaction processes in the tropical Pacific (Bjerknes, 1969), setting up a new era when these two phenomena are studied as a whole from simple statistical analyses to dynamical diagnostic analyses. The interactions among SST, sea surface winds and the thermocline were illustrated in the tropical Pacific (the so-called thermocline feedback or Bjerknes feedback; Bjerknes, 1969; Wyrtki, 1975). In the 1970s and 1980s, highly idealized conceptual models and simplified ocean-atmosphere coupled models were formulated and used to understand ENSO processes, including oceanic equatorial waves (the Kelvin waves and Rossby waves), instability theories of coupled ocean-atmosphere interactions in the tropics, and so on (Pacanowski and Philander, 1981; Battisti and Hirst, 1989). As extensive observational data became available at that time, they were analyzed to characterize a canonical pattern of ENSO evolutions and cycling from the onset, to the developing, mature, and decaying stages (Rasmusson and Carpenter, 1982). Several theories have been proposed to explain the ENSO cycling nature involving positive and negative feedbacks within the tropical Pacific climate system, including the delayed action oscillator (Suarez and Schopf, 1988; Battisti and Hirst, 1989), the recharge/discharge oscillator (Wyrtki, 1975; Jin, 1997), the western Pacific paradigm (Weisberg and Wang, 1997), and the convergence and advection oscillators (Picaut et al., 1997). The advances in the theoretical understanding of ENSO processes and mechanisms lead to model developments; essential processes underlying the ENSO phenomenon have been identified and incorporated into coupled ocean-atmosphere models.

A series of international collaboration programs and initiatives have greatly advanced ENSO studies, including the Tropical Oceans and Global Atmosphere (TOGA) Program and the Climate Variability and Predictability (CLIVAR) Programs etc. (e.g., Trenberth et al., 1998; Wallace et al., 1998). For example, a hierarchy of ENSO models has been developed to synthesize observations and process understanding, forming a solid mathematical basis for ENSO simulations. In terms of complexity, coupled models can be classified into highly idealized conceptual models, simplified intermediate coupled models (ICMs) and hybrid coupled models (HCMs), and fully coupled general circulation models (CGCMs), respectively. In particular, various coupled models have been successfully used for ENSO simulation and prediction (Zebiak and Cane, 1987; Barnett et al., 1993; Chen et al., 1995; Zhang et al., 2003, 2005b; Barnston et al., 2012; Zhang and Gao, 2016). The combined progresses in observations, process understanding, and model developments lead to realizations of ENSO real-time predictions using coupled models; now more than 20 models have been routinely used to make real-time ENSO predictions 6 months to 1 year in advance (see https://iri.columbia.edu/our-expertise/climate/forecasts/enso/current/). Overall, the current coupled models with data assimilation techniques can provide effective predictions of ENSO-related warm and cold SST anomalies in the tropical Pacific 6-12 months ahead (Chen et al., 1995; Fedorov et al., 2003; Jin et al., 2008; Barnston et al., 2012; Zhu et al., 2012; Zhang and Gao, 2016).

As our understanding of ENSO has been getting deepened and more practice of real-time ENSO predictions has been carried out extensively, new challenges emerge in ENSO modeling. For example, ENSO is observed to exhibit great diversity and complexity (Zhang et al., 1998; Yu and Kao, 2007; Yeh et al., 2009, 2014; Capotondi et al., 2015; Chen et al., 2015; Feng et al., 2015; Timmermann et al., 2018; Xie et al., 2018), and it can be modulated by various forcing and feedback effects (e.g., Kang et al., 2017a). It has been recognized that ENSO can have different flavors with ENSO events evolving differently from canonical events as described in Rasmusson and Carpenter (1982). In 1982–1983, for instance, an unexpected strong El Niño event took place in the tropical Pacific; this event evolved strikingly differently from previous canonical El Niño events, including its onset time and the way SST anomalies propagated across the equatorial Pacific. Also, the recent strong 2015 El Niño evolved strikingly differently from the previous extreme events of 1997–1998 and 1982–1983 in terms of formation and warming pattern. For example, the 2015 El Niño event was preconditioned by a weak warming in the central equatorial Pacific in early 2014, a pause in the warming in mid and late 2014, and then amplified dramatically in spring 2015 as second-year warming, which peaked in late 2015 (e.g., Zhang and Gao, 2016; Zhu et al., 2016). Currently, the occurrences of these different types of El Niño cannot be satisfactorily explained by classical ENSO theories; yet current coupled models still have difficulties in depicting ENSO diversities (Zhang et al., 1998; Cai et al., 2018). Many factors have been identified that can affect the way ENSO evolves, including westerly wind bursts, the global warming effects and the interactions between natural climate variability and human-induced climate changes (Lian et al., 2014; Chen et al., 2015; Tan et al., 2020). However, specific processes responsible for ENSO changes have not been understood well, presenting a great challenge to process understanding and accurate modeling of ENSO.

Indeed, large biases and uncertainties still exist in ENSO simulations and predictions. For several decades, some systematic biases have been persisted in climate models, as indicated in the recently released Coupled Model Intercomparison Project Phase 6 (CMIP6) simulations participating in the Intergovernmental Panel for Climate Change (IPCC) Assessment Report (Eyring et al., 2016). For example, the state-of-the-art coupled models still experience noticeable large model biases in simulations of mean state, seasonal variations and interannual variability associated with ENSO. In terms of ENSO prediction, coupled models indicated substantial discrepancies in real-time prediction. For instance, there was a false alarm in predicting the 2014–2016 warm conditions in the tropical Pacific, which embarrassed the scientific community. As described above, in reality, the weak warm SST anomalies observed in early 2014 did not evolve into a strong El Niño event in late 2014; but a strong El Niño event was predicted to occur in late 2014 by many coupled models. For the actual occurrence of one of the strongest events in history, the predictions of the 2015 El Niño event using coupled models exhibited large biases in its onset throughout spring 2015 and the rapid warming in mid-2015 (Zhang and Gao, 2016; Zhu et al., 2016). Additionally, the predicted intensity of the 2015 event in summer and fall 2015 spread widely across coupled models, which indicate that there exist large intermodal differences and uncertainties in ENSO simulations. Further investigations reveal that there are large differences in El Niño prediction skills from one decade to another, with low predictability in the 2000s compared with the 1980s, which may reflect decadal changes in ENSO predictability (Chen et al., 1995; Kirtman and Schopf, 1998; Tang et al., 2008; Hu et al., 2013; Zheng et al., 2016; Luo et al., 2017). However, specific responsible processes have not been understood well.

In summary, various coupled models, being as the most powerful tool for simulating and predicting ENSO, have been developed for uses in ENSO studies, which differ in dynamical formulations, processes parameterizations, resolution, and so on. However, current models still have difficulty in depicting ENSO diversity, whose processes have not been identified clearly. Systematic biases and large uncertainties still exist in ENSO simulations and predictions in coupled models. The inability for current coupled models to make real-time ENSO predictions presents a great challenge to ENSO modeling, indicating an urgent need for model improvements. All these challenges have stimulated a continued interest in ENSO research in the world. Indeed, dedicated efforts have been devoted to model developments for ENSO in the past several decades. Specifically, various ICMs, HCMs, and CGCMs have been carefully tuned so that they have been able to realistically depict ENSO evolution and the related climate variability. In particular, CGCMs without bias and/or flux corrections have been now available and used for different purposes with a variety of applications to climate modeling as evident in the CMIP6 simulations. Furthermore, the integrated efforts of observing, process understanding and model developments of ENSO have led to capacity building of real-time ENSO and related climate predictions (Mu and Ren, 2017).

The motivation of this article is to review ENSO model developments and its simulations. The large numbers of ENSO-related modeling studies make it impossible for us to do so in this brief review. Rather, our primary focus here is placed on the progress recently made at some selected institutions in China, which have been continually and actively pursuing ENSO model developments with a variety of successful applications. As evident recently, the model performances for ENSO and climate simulations have been steadily improved, including physical parameterizations, fining spatial and temporal resolutions, and others. In the following, several coupled models having their reasonably good ENSO simulations at a few institutions in China are reviewed. Note that there are several other ENSOrelated review papers focusing on different aspects of ENSO in the literatures, including McCreary and Anderson (1991), Webster and Yang (1992), Delecluse et al. (1998), Neelin et al. (1998), Stockdale et al. (1998), Philander (1999), Latif et al. (2001), AchutaRao and Sperber (2002), Guilyardi et al. (2004), Jin et al. (2008), Barnston et al. (2012), Capotondi et al. (2015), Tang et al. (2018), Timmermann et al. (2018), Wang (2018), Fang and Xie (2020), Ren et al. (2020), and others.

The rest of the paper is organized as follows. Section 2 presents an overview of the coupled modeling of ENSO. In terms of complexity, ICMs, HCMs, and CGCMs are described in Sections 3, 4 and 5, respectively. Summary and discussion are given in Section 6.

2 AN OVERVIEW OF COUPLED OCEANATMOSPHERE MODELING FOR ENSO

Numerical models provide the mathematical basis for ENSO studies. As ENSO simulation skills are strongly model dependent on dynamical formulations, process parameterizations, resolution, and so on, various coupled models in the past several decades have been developed to represent ocean-atmosphere interactions responsible for the existence of ENSO. The most complicated models used for ENSO studies are primitive equations-based atmosphere general circulation models (AGCMs) and ocean general circulation models (OGCMs), which include comprehensive process parameterizations. Meanwhile, considering the ENSO-related variability characteristics, atmospheric and oceanic models can be effectively simplified for use in ENSO studies not only for computational efficiency, but also for facilitations of physical understanding and process illustrations.

On the atmospheric side, statistical and dynamical models have been used. Statistical atmospheric models can be adequately adopted in ENSO modeling because the ENSO-related atmospheric anomalies can be treated as a feedback response to interannual SST variability in the tropical Pacific. For example, statistical models, constructed from historical data using statistical methods, have been successfully applied to the coupled modeling of ENSO. Another justification of using statistical atmospheric model for ENSO studies lies in the difficulties in representing atmospheric convection processes in an AGCM that could bring uncertainties and errors in wind simulations and correspondingly in ENSO simulations and predictions (Zhu et al., 2017); so the related problem could be implicitly avoided by using a statistical atmospheric model. In terms of atmospheric dynamical models, AGCMs are primitive equationbased comprehensive models, which are very timeconsuming to run. Historically, simplified atmospheric dynamical models have been widely used for representing atmospheric characteristic responses to ENSO-induced SST variability in the tropical Pacific, including the shallow water equation-based model. In particular, considering the vertical structure of atmospheric anomalies during ENSO evolution, a two-layer approximation can be taken in the vertical direction to depict the characteristic atmospheric responses to SST anomalies, which can be solved by using the shallow water equation-based model; this has been clearly illustrated by Gill (1980), Zebiak and Cane (1987), and others. Practical simulations indicate that these simplified atmospheric models work extremely well for ENSO modeling, being able to capture essential atmospheric responses associated with ENSO evolution.

On the oceanic side, dynamical models need to be considered for representing ENSO processes because of the importance of oceanic dynamics to ENSO; various dynamical ocean models have been formulated to study ENSO, including comprehensive OGCMs that are complicated in its formulation and expensive to run. Basically, ocean models can be divided into level ocean models (e.g., Modular Ocean Model (MOM) Version 3 (MOM3); Pacanowski and Griffies, 1998) and layer ocean models (e.g., Gent and Cane, 1989), in terms of vertical coordinates adopted in their dynamical formulations (e.g., Zhang and Zebiak, 2002; Zhu and Zhang, 2018). Different simplifications for ENSO modeling have been effectively made for the ocean, including approximations of reduced gravity, shallow water and linearization etc. (e.g., Moore and Philander, 1977; Busalacchi and O'Brien, 1980). For example, the existence of a well-defined vertical structure of upper-ocean thermal responses makes a two-layer approximation work well for simplifying the complicated ocean responses to surface winds in the vertical. In particular, a modal decomposition in the vertical can be performed to simplify the problem solving for upper-ocean responses to atmospheric wind forcing, with only a few vertical modes being retained to adequately depict main response features (e.g., Zebiak and Cane, 1987; Keenlyside and Kleeman, 2002). Again, such simplified considerations in the ocean modeling for ENSO are not only for computational efficiency, but also for being amenable for clear illustrations of processes involved.

These individual oceanic and atmospheric models are used to perform ocean-only and atmosphere-only experiments for ENSO-related simulations. When forced by prescribed atmospheric forcing fields (surface wind stress, heat flux and freshwater flux), ocean models can capture SST anomalies associated with El Niño and La Niña events and their transitions. When forced by prescribed oceanic forcing fields (i.e., SST), atmosphere models can depict surface wind stress fields, etc. in association with the SO. As ENSO originates from the ocean-atmosphere coupling in the tropical Pacific, coupled models are required for use in representing ENSO. When the ocean and atmosphere models come to be coupled together; however, good performances of each oceanic and atmospheric component model cannot guarantee that ENSO can be adequately captured well in their coupled cases. In fact, the way ENSO can be depicted in coupled models depends not only on the performances of individual component model of the ocean and atmosphere, but also collectively on coupled interactions arising from the tropical Pacific. Various factors have been identified that can affect ENSO performance in coupled models. Practically, atmospheric and oceanic models need to be very carefully tuned when they are coupled together so that responsible processes for ENSO cycling can be adequately represented in a balanced way, including the well-defined phase relationships between interannual anomalies, various positive and negative feedbacks, and underlying processes etc. In particular, as a critically important component, great cares have been devoted to ocean models that are tuned specially for ENSO process representations, being able to sustain ENSO cycles in the climate system of the tropical Pacific.

Combinations of these differently formulated oceanic and atmospheric models lead to different types of coupled models in the past several decades, which are widely used to investigate oceanatmosphere interactions and ENSO. In terms of complexity, coupled models can be categorized into ICMs, HCMs, and CGCMs (e.g., Zebiak and Cane, 1987; Kleeman, 1993; Latif and Yamagata, 1993; Chen et al., 1995; Latif et al., 1998; Stockdale et al., 1998; Zhang and Gao, 2016). As the most complicated coupled models, CGCMs are referred to as a coupled model in which primitive equation-based AGCMs and OGCMs are both used for the atmosphere and ocean. Relative to CGCMs which are computationally expensive, ICMs and HCMs are simplified models, in which a kind of simplified anomaly model is adopted for the atmospheric or oceanic component, respectively. ICMs are formed by coupling both simplified oceanic and atmospheric models. For example, Zebiak and Cane (1987) formulated an ICM in which the Gill-type atmospheric steady dynamical model is coupled to a linear two-layer oceanic dynamical model; this is the first dynamical model used for ENSO prediction experiments (Cane and Zebiak, 1985; Cane et al., 1986). In a more complicated configuration, HCMs are also simplified models, in which one component is taken as a simple anomaly model, whereas the other component is general circulation model (GCM)-based. For instance, an HCM can be formed by coupling an atmospheric statistical model with an OGCM (e.g., Zhang, 2015a, b). On the other way around, another type of HCMs can be formulated by coupling a simplified dynamical oceanic model with an AGCM (Hu et al., 2019).

There are obvious advantages and disadvantages when using these various coupled models with different degrees of complexity, ranging from simplified to CGCM-based models. Computing efficiency is one aspect to be considered, as well as model performances and facilitating clear and easy ways for physical understanding and process illustrations. For example, consisting of OGCMs and AGCMs, CGCMs include comprehensive processes and are susceptible to climate drift because the full coupling is realized between the ocean and atmosphere. In ICMs and HCMs, one component is constructed to be an anomaly form and anomaly coupling is adopted, being very helpful for avoiding climate drift problems and having good model performances of ENSO simulations and predictions. So, simplified coupled models of ENSO have been widely used not only for time efficiency, but also for being amenable to process understanding and interpretations as reflected in Gill (1980), Zebiak and Cane (1987), and others.

Intensively dedicated efforts have been devoted to model developments for ENSO in China in the past several decades. Currently, various coupled models have been available with a variety of successful applications, including physical understanding and realistic simulations. For instance, coupled models have been used successfully for ENSO simulations and participations in the CMIP6. Now, the state-ofthe-art models are used to examine a variety of forcing and feedback mechanisms responsible for ENSO modulations, depict the spatial structure and time evolution of ENSO, the well-defined relationships among oceanic and atmospheric fields, and even to make real-time ENSO predictions, and climate predictions and projections. In the following, several examples of ENSO models that currently have been used actively for ENSO simulations and predictions in China are selected and described in detail.

3 ICM SIMULATIONS OF ENSO

Intermediate coupled models (ICMs) are referred to as a particular type of simplified ENSO models, which are distinguished by being intermediate in complexity between highly idealized conceptual ENSO models and comprehensive CGCMs (Zebiak and Cane, 1987; Kang and Kug, 2000; Zhang et al., 2003; Song et al., 2018). One basic feature of such an ICM is its anomaly formulation in which only interannual perturbation fields of the atmosphere and ocean are calculated directly, whereas seasonally varying climatological fields are specified from the corresponding observations. An apparent advantage of such an anomaly-focused approach lies in the fact that many processes with primary importance to the mean climatology (e.g., deep ocean processes, small and meso-scale processes etc.) may be not critical to ENSO and thus there is no need to explicitly represent their relative contributions to ENSO modeling. That is to say, for ENSO-related interannual modeling, there is no explicit need for caring much of processes responsible for the mean climatology; rather, only interannual relationships between atmospheric and oceanic anomalies and the underlying processes are represented in these simplified models.

An ICM can be formed by coupling a simple atmospheric model with a simple oceanic model. For the atmosphere, dynamical and statistical models are adopted. For instance, a simple two-layer approximation can be taken to represent the vertical structure of atmospheric responses to SST anomalies as in the Gill (1980) model; also, statistical atmospheric models can be used for ENSO modeling by constructing statistical relationships between interannual anomalies of SST and other atmospheric fields from historical data using common statistical methods (i.e., singular value decomposition (SVD)). For the ocean, ocean models need to be dynamical because oceanic dynamics are critically important to ENSO. Considering the response characteristics in the upper ocean of the tropics, simplified dynamical ocean models have been adopted. For example, the upper ocean can be approximated as a two-layer structure, whose dynamics have been understood well (Gill, 1980) and numerical solutions can be easily obtained by using the shallow water equation-based model. In addition, within the ICM framework, full ocean-depth thermodynamics are not necessary to be considered; instead, only SST is taken into account to represent the surface layer thermodynamics. Further, an SST anomaly model can be formulated to depict interannual SST variability, which is embedded into a simplified dynamical ocean model.

These simplified component models are coupled through an anomaly coupling for the ocean and atmosphere: SST anomalies, directly produced by its anomaly model, force the atmosphere to produce surface wind variations, which in turn affect the ocean (including the ocean circulation and thermocline), leading to changes in SST. As crucially important to SST, subsurface and thermocline effects on SST need to be adequately represented. Within the ICM context, the subsurface entrainment temperature (T_e) into the surface mixed layer is an explicitly expressed variable, which can be parameterized in terms of the thermocline variability (as represented by sea level (SL)) due to their close relationships (Wyrtki, 1975).

One well-known ICM is that developed by Zebiak and Cane (Zebiak and Cane, 1987; ZC87 thereafter), the first dynamical model that is used for ENSO prediction (Cane et al., 1986) and has been pioneering in ENSO studies. More specifically, the ZC87 adopted a two-layer linear shallow water approximation to represent the ocean-atmosphere system in the tropical Pacific. Atmospheric component is a steady dynamical model based on the Gill model; the ocean component is based on the linear shallow water dynamics. A nonlinear SST submodel is incorporated into the dynamical ocean model to represent the effects of various processes on SST; the subsurface entrainment temperature (T_e) anomaly field is parameterized by its analytical relationship with the thermocline variability. As an anomaly model, the interactions between the ocean and atmosphere are realized by anomaly coupling, in which the mean seasonal climatologies of SST and vertical gradient of temperature etc. are prescribed from observations. Therefore, a modeling focus is placed on primary processes important to SST interannual variability. Such a simplified coupled model can bypass the socalled climate drift problems, which are commonly seen in CGCMs. As have been demonstrated, the ZC87 can very successfully depict SST variability associated with El Niño (e.g., Zebiak and Cane, 1987), and thus has been widely applied to ENSO studies (e.g., Mu et al., 2007). Moreover, the ZC87 has been further simplified into a type of conceptual models, in which the governing equations of ENSO can be further reduced to the maximum degree of simplicity so that semi-analytical solutions can even be obtained, in which only essential features of oceanatmosphere anomalies and their interactions are retained to represent ENSO. These derived conceptual models have been extremely helpful for theoretical analyses and process illustrations (e.g., Jin, 1997; Jin and An, 1999). More recently, the tropical Pacific ZC87 model has been extended to include the Indian Ocean to examine interannual variability in the Indian Ocean and its effects on that in the tropical Pacific (Song et al., 2018).

In China, such a type of ICMs was formulated at the Institute of Oceanology/Chinese Academy of Sciences (IOCAS) for ENSO-related modeling studies (denoted as the IOCAS ICM; Zhang and Gao, 2016). The IOCAS ICM is an anomaly model consisting of an intermediate ocean model (IOM) and a statistical model for atmospheric wind stress anomalies (Fig. 1). Its dynamic ocean component is an intermediate-complexity model developed by Keenlyside and Kleeman (2002) based on a baroclinic modal decomposition method in the vertical direction (McCreary, 1981). In contrast to the commonly used ZC87, this relatively newly developed ICM accounts for the spatially varying stratification and the realistic vertical structure of the upper ocean (e.g., ten vertical modes are included). In addition, nonlinear effects of the momentum equations for ocean currents are partially considered so that the zonal and meridional currents in the equatorial oceans can be more realistically depicted. An SST anomaly model is incorporated into the ocean dynamical model to represent the thermodynamic processes of the surface mixed layer. One striking feature of the IOCAS ICM is the way the temperature of subsurface water entrained into the mixed layer (T_e) is parameterized. An optimized procedure is developed to depict T_e using an inverse approach based on the SST anomaly equation of the ICM and its empirical relationship with SL variability (Zhang et al., 2005a). The ICM has been carefully tuned so that it can reproduce sustained ENSO, involving the interactions between the ocean and atmosphere in the tropical Pacific. The ICM has its big computational advantage relative to HCMs and CGCMs, and yet the performance of ENSO simulations is satisfactorily good (Fig. 2).

The IOCAS ICM has already been successfully used for ENSO simulations and predictions. For example, this model is one of the coupled models that made a good prediction of the cold SST conditions during 2010-2012 in the tropical Pacific (Zhang et al., 2013). After optimizing the model performance in terms of ENSO simulations and retrospective ENSO predictions (Gao et al., 2016; Zhang and Gao, 2016), the IOCAS ICM has been routinely used to predict the SST evolution in the tropical Pacific since August 2015 (for more detail, see https://iri.columbia.edu/our-expertise/climate/forecasts/enso/current/). A realtime ENSO prediction example is given in Fig. 3. Also, this ICM has been used for ENSO-related data assimilation and predictability studies (Zheng et al., 2006; Gao et al., 2016, 2018; Tao et al., 2018; Tao and Duan, 2019; Mu et al., 2019).

Historically, simplified ICMs have played a very important role in the developments of ENSO understanding and predictions. Such a level of simplified model not only allows for realistic ENSO simulations and straightforward comparisons against observations, but also is amenable to model tuning and process illustrations. Certainly, such simplified ICMs with their anomaly formulation have limited applicability, particularly in such a situation where both mean climatology and interannual anomaly fields for the ocean and/or atmosphere need to be determined by the models themselves. Therefore, more comprehensive coupled models need to be used, in which total fields of oceanic and atmospheric states (both climatology and interannual anomalies) can be adequately determined by the model itself. Also, there are many important multiple processes responsible for ENSO modulations that need to be taken into account, such as atmosphere stochastic forcing and ocean biology-induced heating (e.g., Zhang et al., 2008; Lian et al., 2014; Zhang, 2015a; Kang et al., 2017a); so it is necessary to use AGCMs and OGCMs for comprehensive ENSO modeling.

4 HCM SIMULATIONS OF ENSO

As categorized to be the so-called HCM, its level of complexity lies between ICMs and CGCMs (Neelin et al., 1992; Neelin and Jin, 1993; Barnett et al., 1993; Syu et al., 1995; Chang et al., 2001; Tang, 2002; Zhang and Busalacchi, 2009; Zhu et al., 2009, 2013; Zhang, 2015b). Basically, two types of HCMs can be formulated, depending on whether an AGCM or OGCM is used for the atmospheric component or oceanic component, respectively. Here, HCM^AGCM is referred to as an HCM in which an AGCM is coupled with a simplified dynamical ocean model; HCM_OGCM is referred to as an HCM in which an OGCM is coupled with a simplified atmospheric model (statistical or dynamical). Because one component of such an HCM is taken as an anomaly form in which only interannual perturbation fields are considered for the atmosphere or ocean; an anomaly coupling is adopted for representing interactions between the atmosphere and ocean, in which seasonally varying observed climatological fields are prescribed for use. For the HCM^AGCM, the AGCM produces total surface wind fields, whose interannual anomalies are calculated relative to their corresponding climatological fields of the AGCM; the resultant wind anomalies are used to force the simple ocean model to produce interannual SST anomalies directly, which are combined with observed SST climatologies to force the AGCM. For the HCM_OGCM, the OGCM produces total SST fields, whose interannual anomalies are calculated relative to the corresponding climatological SST field of the OGCM; the resultant interannual SST anomalies are used to force a simple atmospheric model to produce interannual anomalies of wind stress, which are combined with observed climatologies of wind stress to force the OGCM. Evidently, the HCM shares some advantages with the ICMs, as described above. For example, because one component for the ocean or atmosphere is taken as an anomaly model, an anomaly coupling is realized between the ocean and atmosphere, in which observed climatologies are prescribed; this approach helps to bypass climate drift problem that is commonly faced in CGCM-based simulations. In addition, such a simplified HCM is computationally efficient to run in time integration.

This type of coupled models makes it suitable and highly efficient for addressing some of the fundamental questions regarding coupled behaviors of air-sea interactions in the tropics. For example, as interannual anomaly fields are explicitly separated from total fields, the HCM allows for interannually varying coupled forcing and feedback effects on ENSO to be isolated and examined in a clean way. Neelin et al. (1992) formulated such an HCM_OGCM, which consists of a level OGCM (The MOM type of OGCMs developed by Geophysical Fluid Dynamics Laboratory (GFDL)/NOAA) and an atmospheric steady dynamical model; this HCM was used to investigate fundamental behaviors of coupled interannual variability associated with ENSO, including the relationships between the changes in coupling parameters (representing forcing and feedback intensities) and the ways interannual oscillations are sustained. In a very early stage, this type of model is also used for ENSO predictions (e.g., Barnett et al., 1993). Historically, such type of coupled models has made great contributions to ENSO understanding and prediction studies. Currently, this type of simplified HCMs remains valuable and powerful tools for understanding the tropical air-sea interactions associated with interannual and decadal variability of the climate system. Two types of such simplified HCMs have been formulated in China, which have been used for ENSO modeling; some results are briefly described in this subsection as follows.

4.1 An HCM^AGCM consisting of an AGCM and an intermediate ocean model (IOM)

One HCM has been recently formulated (the right panel in Fig. 1), denoted as HCM^AGCM, in which a simplified intermediate ocean model (IOM) is coupled with an AGCM (Hu et al., 2019). The ocean component used is the oceanic component model of the IOCAS ICM as described above (Fig. 1); the atmospheric component is the ECHAM5, the fifth version of the European Centre for Medium-Range Weather Forecasts (ECMWF) model developed by the Max Planck Institute for Meteorology (Roeckner et al., 2003). These two component models are coupled by employing an anomaly coupling strategy because the ocean model is formulated as an anomaly form (Fig. 1): the IOM directly produces SST anomalies, which are added onto observed SST climatology (SST_clim) to form total SST fields, which are used to force the AGCM; the AGCM produces total wind stress (τ) fields, whose interannual anomalies are extracted relative to the AGCM's climatology (τ_inter=τ–τ_clim), and then are used to force the ocean. In addition to its computational advantage relative to CGCMs, the HCM^AGCM is also climate drift-free because of using anomaly coupling, with which observed seasonal climatology of SST is prescribed from observations and only oceanic perturbation fields are calculated directly from the ocean model; this configuration presents a constraint on the HCM, which has less freedom to evolve. The use of the AGCM allows for considerations of atmospheric stochastic forcing effects and atmospheric teleconnection due to remote dynamical responses. Note that the anomaly ocean model is adopted in the HCM^AGCM configuration, which allows us to focus on primary processes that are only important to SST interannual variability in the tropical Pacific, without the explicit need for caring much of those responsible for the mean climatology in the ocean.

Currently, a regional HCM^AGCM is formulated, in which active ocean-atmosphere coupling is allowed only in the tropical Pacific although a global AGCM is used, which can depict global atmospheric responses to SST anomalies in the tropical Pacific. A long-term simulation (more than 100 years) has been performed using the HCM^AGCM; its simulation skills for ENSO are evaluated using SST over the tropical Pacific (Fig. 2). The HCM^AGCM produces irregular ENSO events with a broad spectrum of oscillation periods between 2 and 5 years, which is attributed to the use of an AGCM. The amplitude and evolution of ENSO events and the phase locking of interannual SST anomalies to the annual cycle are reproduced realistically. Also, the ENSO-related SST effects on interannual variability of the global atmosphere can be well depicted (Hu et al., 2019). Despite the slightly stronger amplitude of SST variability over the central Pacific in the HCM^AGCM compared with the observed, the simulated ENSO-related response patterns in the atmosphere are in good resemblance with the corresponding observations. In particular, the use of the global AGCM in the HCM^AGCM allows for examining remote effects of tropical SST anomalies on weather and climate in the subtropics and midlatitude regions, which are realized through the Hadley cell and wave trains in the atmosphere. Thus, teleconnection patterns of climate anomalies around the globe can be investigated, which are affected by SST anomalies in the tropical Pacific, including the Pacific and North America (PNA) patterns. These simulated atmospheric response patterns in the HCM^AGCM are similar to what is observed in nature. Therefore, this HCM^AGCM can not only depict the ENSO-related interactions over the tropical Pacific, but also reproduce ENSO-induced global atmospheric variability, thereby providing a useful modeling tool for studying ENSO-related coupling in the tropical Pacific and remote influences on the atmosphere globally.

4.2 An HCM_OGCM consisting of an OGCM and an atmospheric statistical model

Another type of HCM is formulated for the tropical Pacific, consisting of an OGCM and a simple statistical model for the atmosphere (Fig. 4), denoted as HCM_OGCM (Zhang and Busalacchi, 2009; Zhang, 2015a). The atmospheric statistical model for wind stress, which is constructed using SVD analyses, is the same as that used in the IOCAS ICM (Fig. 1). The OGCM is based on Gent and Cane (1989), which is a reduced gravity, primitive equation-based layer ocean model: the uppermost layer is considered to be a mixed layer (ML) and the layers below are divided according to sigma coordinates (the ratio of each layer to the total layer thickness is fixed to be constant). Chen et al. (1994) embedded a bulk ML model into the ocean model, in which the depth of the ML (H_m) is considered as a prognostic model variable that is directly predicted using the embedded bulk ML model. These atmospheric and oceanic component models are coupled by employing an anomaly coupling strategy: SST anomalies are obtained from the OGCM simulations relative to the OGCM's climatology (SST_inter=SST–SST_clim); an atmospheric statistical model is used to calculate the corresponding wind stress anomalies (τ_inter), which are combined with observed climatological wind stress (τ=τ_clim+τ_inter) to force the OGCM. In addition to its computational advantage relative to CGCMs, the model is climate drift-free because anomaly coupling is taken with observed climatology of wind stress being prescribed from observations. That is, the use of seasonally varying climatological fields from observations (e.g., τ_clim) in the anomaly coupling acts to impose a constraint on the coupled system, which has less freedom to evolve in HCMs.

Note that this layer ocean model-based HCM has some important features that are different from other level ocean model-based HCMs (e.g., Neelin et al., 1992; Tang, 2002). For example, our layer OGCM explicitly takes into account the mixed layer dynamics: the first layer is treated as a mixed layer, whose depth is predicted explicitly as a model prognostic variable by using a bulk mixed layer model (Kraus and Turner, 1967). One consequence of such a configuration in the layer ocean model is that the way atmospheric forcing is applied to the upper ocean is different from level ocean models (Zhang and Zebiak, 2002). For instance, in the layer ocean model-based HCM, the ML depth is directly affected by wind stress and buoyancy flux (Q_B; the combination of heat flux and freshwater flux). In contrast, the effect of Q_B on the depth of the ML is not explicitly affected in the level ocean model based-HCM (e.g., the MOM3). The effects of freshwater representations (including Q_B) on ocean and coupled simulations have been examined (Kang et al., 2017b; Gao et al., 2020)

The use of the comprehensive OGCM for the HCM_OGCM allows for ocean processes-focused modeling studies. For example, we formulated one regional HCM to represent interactions among the atmosphere, ocean physics and ocean biogeochemistry (AOPB) in the tropical Pacific (Zhang et al., 2018). As shown in Fig. 4, the atmospheric component is a statistical feedback model representing responses to SST anomalies, and the oceanic component model also includes ocean biological processes. Based on this HCM AOPB, various experiments have been performed to illustrate ENSO modulations induced by oceanic processes, including ocean chlorophyllinduced modulating effects on the penetration of shortwave radiation in the upper ocean. The HCM_OGCM has been used to systematically investigate ocean biology-induced heating feedback onto the climate, including biological effects on ENSO at different feedback intensities. Figure 5 shows an example for the simulated interannual anomalies of some fields from the HCM AOPB (Tian et al., 2018; Zhang et al., 2018). As evident from the HCM_OGCM, ENSO simulated is quite regular because the atmospheric model adopted is a statistical model, in which stochastic wind forcing is not taken into account. In addition, wind stress is considered to be a feedback response to SST variability, without considering atmospheric dynamical processes explicitly and remote teleconnetion effects in the atmosphere. In addition, the active coupling is only allowed within the tropical Pacific; so, such a model configuration does not permit the investigation of global atmospheric teleconnection that can be affected by SST variability in the tropical Pacific. In this regards, AGCMs need to be used for the coupling, in which the influences of atmospheric stochastic forcing and the comprehensive atmospheric dynamics and teleconnection can be adequately represented. These processes can be taken into account by using CGCMs.

5 CGCM SIMULATION

Moving to the most complicated category of coupled models, fully coupled GCMs (CGCMs) are formed by coupling an AGCM with an OGCM in a straightforward way. Total fields of SST and surface winds etc. are completely determined by the CGCM itself, which are directly used for the coupling between the ocean and atmosphere: wind stress, heat, and freshwater fluxes at the interface are used to drive the ocean, whereas total SST fields are used to force the atmosphere. The CGCM-based simulations can provide comprehensive descriptions of the oceanic and atmospheric states. Compared with other types of simplified models (e.g., ICMs and HCMs as described above), one obvious disadvantage in using CGCMs which include comprehensive processes is that CGCM-based simulations are very expensive to perform in time integration. In particular, the enormous computer resources demanded by CGCMs make it difficult to conduct numerous simulations needed to examine interannual behaviors associated with ENSO. Obviously, the complexities in the AGCM, OGCM and their fully coupled configuration also make it difficult to clearly delineate processes involved in ENSO evolution and the coupled interactions among different components of the climate system. Practically, it is quite difficult for CGCMs to have a good performance for ENSO simulations because climate drift problems can occur to mix up the numerical solutions. Indeed, intensive cares are highly required to tune CGCMs for having reasonable ENSO simulations.

In the early developmental stage, CGCMs suffer from serious problems in coupled climate modeling, including climate drift and systematic biases in simulations of mean state, seasonal, and interannual variations (e.g., Zhang et al., 1995). Even at present, some major challenges are still faced in many CGCM simulations, including the double Intertropical Convergence Zone (ITCZ) problem, cold SST biases in the eastern equatorial Pacific, and others. The existence of these problems can be attributed to the fact that CGCMs are fully and directly coupled between the ocean and atmosphere; no constraint is imposed on the coupled system, in which the climate states in the CGCMs can have the freedom to evolve. To avoid climate drift, correction or anomaly coupling can be adopted. In terms of ENSO modeling, the performance of CGCMs is still not satisfactory at present. These simulation discrepancies indicate that some natural processes associated with forcings and feedbacks are not represented well in the coupled system.

In the past, intensive and extensive developmental efforts have been devoted to improvements in CGCMs at some institutions in China, making a continual and steady progress. In the past several decades, CGCMs have been carefully tuned so that they are able to capture self-sustained interannual oscillations associated with ENSO. Recently, great progress has been made in CGCMs, as represented by the recent release of CMIP6 in which full coupling between the ocean and atmosphere is executed without flux and bias corrections. The state-of-the-art CGCMs have a variety of successful applications to climate modeling, including interannual variability studies associated with El Niño, climate simulations and projection studies associated with future climate change, and so on. Such modeling efforts have been underway at several institutions in China and CGCM-based ENSO studies have been conducted, which will be briefly described below.

5.1 CGCM-based activities at LASG/IAP

In the past several decades, dedicated model developments have been pursued at the state key Laboratory of numerical modeling for Atmospheric Sciences and Geophysical fluid dynamics (LASG), Institute of Atmospheric Physics (IAP), Chinese Academy of Sciences (e.g., Zeng, 1979). OGCM and AGCMs have been developed at LASG/IAP, which have some unique characteristics different from others in terms of the dynamical formulation and numerics, including the calculation of the departures of thermodynamic variables by subtracting the standard stratification in the vertical, total available energy consideration, and so on.

For example, the IAP OGCM is the first free surface OGCM that was developed in the world, with successful applications to simulations of the basinscale ocean general circulation and climate (e.g., Zeng, 1979). Note that before the late 1980s, a rigid lid approximation was taken in basin-scale OGCMs used for the ocean and climate modeling (e.g., Haney, 1974; Bryan et al., 1975; Cox, 1975; Manabe et al., 1979; Han, 1984; Philander and Seigel, 1985; Latif, 1987; Philander et al., 1987; Rosati and Miyakoda, 1988; Gordon and Corry, 1991; Semtner and Chervin, 1992). That is, the ocean surface is considered to be fixed, which is thus called as rigid-lid. Correspondingly, the sea surface boundary condition for vertical velocity in the rigid-lid OGCMs is set to zero and thus gravy waves are filtered out. In reality, however, the sea surface moves up and down, changing with time and space. Physically, as pointed out by Zeng et al. (1991), the rigid-lid approximation is artificially introduced into the ocean models, which thus can induce errors in ocean simulations. Setting vertical velocity to be zero is equivalent to removing the divergence of vertical mean flow from all scales of motion, including the large-scale barotropic Rossby waves. As a result, this artificial constraint excludes the sea surface height-related available surface energy part and its conversion to kinetic energy, possibly leading to a distortion of the energy conversion and cycle in the ocean. Consequently, it might introduce errors in the computation of surface elevation and currents, propagation of very long waves, and so on. Mathematically, the artificially introduced rigid-lid approximation presents a constraint on the ocean model with non-divergence of the vertically averaged mean flow; a special computational procedure is needed to obtain barotropic component of the solution of OGCMs with the rigid-lid approximation.

To address these issues existing in the rigid-lid OGCMs used at that time, a free surface OGCM was formulated and successfully used for basin-scale ocean simulations in the late 1980s at LASG/IAP, CAS (e.g., Zeng et al., 1991). The free surface IAP OGCM allows for removing the rigid-lid surface condition and thus permits explicit prediction of sea level topography as a prognostic variable. To achieve stability and efficiency in solving the primitive equation-based OGCM without the rigid-lid approximation, several time integration schemes were also developed for the free surface IAP OGCM, including a barotropic-baroclinic mode-splitting algorithm, the semi-implicit scheme for time integration. Also, the spatial finite-difference schemes were formulated, which keep total potential available energy conservation and its conversion to other energy components (Zeng, 1979), which presents effective constraints on long-term time integration for such a free surface ocean model. The performance of free surface OGCMs had been assessed by conducting ocean-only and coupled ocean-atmosphere modeling experiments (e.g., Zhang and Liang, 1989; Zeng et al., 1991; Zhang et al., 1992; Zhang and Endoh, 1992, 1994). At first, such a free surface OGCM is thought to be enormously time consuming; but it turns out that using numerical schemes mentioned above, the computation of the free surface OGCM could be carried out at least as efficiently as that of the rigid-lid one. The IAP OGCM has had a variety of applications, including simulations of annual mean circulation of the world ocean (Zhang and Liang, 1989), of annual mean and seasonal cycle and interannual variability in the tropical Pacific Ocean (Zhang and Endoh, 1992, 1994), and so on. The simulated results show that the IAP OGCM is able to reproduce not only the observed current system and thermal structure and their variations, but also the actual sea level in the tropical Pacific, which shares many satellite-derived sea surface features.

Following the success in the free surface OGCM activities at LASG/IAP, other ocean modeling groups began to develop free surface OGCMs in the early 1990s. Some techniques developed for the freesurface IAP OGCM are adopted by other institutions to develop OGCMs without the rigid-lid approximation (e.g., Killworth et al., 1991; Dukowicz and Smith, 1994). For example, in the release of the MOM3 at GFDL/NOAA (Pacanowski and Griffies, 1998), one new improved feature in MOM3 relative to its previous version (MOM2) includes the implementation of an explicit free surface treatment. Now, OGCMs with the free surface have been widely used in the world.

Based on the free surface OGCM framework, great progress has been made in coupled climate modeling at LASG/IAP. A series of CGCMs have been developed, which have contributed to the multiple phases of the Coupled Model Inter-comparison Project (CMIP; Taylor et al., 2012; Eyring et al., 2016). A family of Flexible Global OceanAtmosphere-Land System Model (FGOALS) CGCMs has been developed for climate modeling at LASG/IAP. Specifically, three generations of the FGOALS CGCMs have been developed at LASG/ IAP since 2002. The first generation is FGOALS-g1.0 (Yu et al., 2004, 2008, 2011); its oceanic component is the IAP OGCM used in the LASG/IAP Climate System Model version 1.0 (LICOM1.0; Liu et al., 2004) and its atmospheric component is the Gridpoint Atmospheric Model of IAP/LASG version 1.0 (GAMIL1.0; Wang et al., 2004). Then, the second generation of FGOALS CGCMs includes two versions, respectively denoted as FGOALS-g2 (Li et al., 2013b; Yu et al., 2013) and the other denoted as FGOALS-s2 (Bao et al., 2013), differing in the atmosphere; here, -g stands for the grid-point model and -s stands for the spectral model, respectively. Here, the oceanic component is the same as LICOM2.0 (Liu et al., 2014), but two different atmospheric components are used, one called as GAMIL2.0 (Li et al., 2013a) for FGOALS-g2 and the other called as the Spectral Atmospheric Model of IAP LASG version 2.0 (SAMIL2.0) (Bao et al., 2013; Liu et al., 2014) for FGOALS-s2, respectively. Now, the third generation has been recently developed for two parallel versions FGOALS-g3 (Li et al., 2020) and FGOALS-f3-L (Guo et al., 2020), in which the same oceanic component model (denoted as LICOM3.0; Yu et al., 2018) and two atmospheric component models (denoted as GAMIL3.0 and the Finite-volume Atmospheric Model of the IAP/LASG version 3.0 (FAMIL3.0); Li et al., 2019, 2020) were respectively employed.

These different versions of the FGOALS CGCMs developed at LASG/IAP since 2002 have had a variety of applications to climate simulations and projections. In particular, they have contributed to the multiple phases of the CMIP. For example, the coupled model FGOALS-g1.0 contributed to CMIP3 (Yu et al., 2004, 2008, 2011), the two versions FGOALS-g2 and FGOALS-s2 contributed to CMIP5, and the recent two versions FGOALS-g3 and FGOALS-f3-L contributed to CMIP6 (Lin et al., 2019), respectively.

As an interest for this review article, we here focus on ENSO simulations and some examples for ENSO modeling studies are given below. To synthesize the overall performance of ENSO simulations in the three generations of the FGOALS CGCMs, we present comparisons for some basic ENSO features between the observation and simulations. To facilitate the comparison, g1, g2, s2, g3 and f3-L were denoted to as FGOALS-g1.0, FGOALS-g2, FGOALS-s2, FGOALS-g3, FGOALS-f3-L, respectively; here, the number indicates the version of the FGOALS CGCMs, and -g stands for grid-point AGCM and -s stands for spectral AGCM, respectively.

Figure 6 shows the spatial pattern of the standard deviation of SST variability over the tropical Pacific. The degree to which ENSO simulations are evolutionarily improved can be clearly seen from early version FGOALS-g1.0 to the next version. For example, although the early versions can depict interannual variability associated with ENSO, large biases exist as shown in Fig. 6. FGOALS-g1.0 shows extremely strong and regular ENSO, characterized by a single peak at approximately 3 years. The SST variability simulated by the first generation of FGOALS (i.e., as indicated by g1 in Fig. 6) shows severely overestimated interannual variability over the central and eastern equatorial Pacific (Fig. 6b), indicating the overly large ENSO amplitude.

After its improvements, FGOALS-g2 shows improved performance in simulating the overall ENSO behaviors. As clearly seen in Bellenger et al. (2014) for comparisons among the CMIP5 models, FGOALS-g2 exhibits excellent performance in simulating ENSO features, including its amplitude, the spatial structure and spectrum of SST variability, and phase-locking with seasonal cycles. Chen et al. (2016) analyzed the reasons for the improvement of ENSO simulation in FGOALS-g2 compared with its predecessor FGOALS-g1.0 in detail. The improved simulation of ENSO amplitude can be mainly attributed to the reasonable representation of the thermocline and thermodynamic feedbacks, which could be traced back to the more realistic basic state simulations in FGOALS-g2 than those in FGOALS-g1. In particular, the simulated deeper mean thermocline in FGOALS-g2 results in a weakened thermocline response to the zonal wind stress anomaly, and the looser vertical stratification of mean temperature leads to a reduced response of anomalous subsurface temperature to anomalous thermocline depth, both of which cause the more realistic thermocline feedback in FGOALS-g2. On the other hand, the alleviated cold bias of mean SST leads to more reasonable thermodynamic feedback in FGOALS-g2. In addition, the change from the regular ENSO oscillation in FGOALS-g1.0 to the somewhat irregular ENSO oscillation in FGOALS-g2 primarily arises from the fact that the atmospheric highfrequency "noise" (e.g., westerly wind bursts) is simulated relatively well in FGOALS-g2.

In general, all the five of FGOALS CGCMs capture the spatial distribution features of ENSO as observed. The specific standard deviations of the Niño3 index and Niño4 index derived from the observation and the five simulations are listed in Table 1. The improvements from one version of the FGOALS CGCMs to the next version are striking. For example, the overestimation of ENSO amplitude is alleviated to some extent in s2 and f3-L (Fig. 6c & e). It is evident that the ENSO amplitude reproduced by both g2 and g3 (Fig. 6d & f) is now comparable to the observation, compared with the early version.

We further present the standard deviation of Niño3.4 for each calendar month (Fig. 7), which can be used to indicate the ENSO phase locking feature. As shown in Fig. 7a, ENSO-related SST variability usually peaks during boreal winter. It is exciting to note that all the five FGOALS CGCMs capture the feature of ENSO phase locking very well. Note that more than half of the CGCMs in CMIP5 archive have difficulty in reasonably representing the ENSO phaselocking feature (Ham and Kug, 2014). Figure 7b shows the power spectra of the Niño3.4 index for the observation and five simulations from the corresponding versions of the FGOALS CGCMs. The power spectra of the observed Niño3.4 index are characterized by multiple peaks at around 3–5 years, while the counterparts in g1 and f3-L are characterized by a single peak at approximately 3 years. Such deficiency in ENSO periodicity simulation has been alleviated in the other three versions (i.e., s2, g2, and g3) to some extent.

Some biases still exist and further improvements are underway. For example, the spatial structure of SST variability indicates some discrepancies compared with observations. Also, all the five simulations cannot reproduce the observed positive skewness well, as shown in Table 1. This indicates that the ENSO asymmetry is underestimated in all the simulations. It is worth mentioning that the underestimation in ENSO asymmetry still remains prevalent in almost all the CGCMs in CMIP5 (Yu et al., 2013). At the same time, the FGOALS CGCMs have had a variety of successful applications, which can be found in the related publications (Yu et al., 2013); the CMIP-related applications and products are available from the CMIP website for the FGOALS CGCM simulations.

5.2 CGCM-based activities at FIO

One CGCM was developed at the First Institute of Oceanography (FIO), the Ministry of Nature Resources of China; its Earth System Model (ESM) version 1.0 (FIO-ESM v1.0) has been successfully used for climate modeling, including participation in the CMIP5 (Qiao et al., 2013). One striking feature with the FIO-ESM v1.0 is that the ocean surface wave effects occurring at the air-sea interface are explicitly taken into account, which is realized by embedding an ocean surface wave model into an OGCM to represent the non-breaking surface wave-induced vertical mixing (Qiao et al., 2004). Specifically, the FIO-ESM v1.0 consists of the Community Atmosphere Model Version 3 (CAM3), the Community Land Model Version 3.5 (CLM3.5), the Los Alamos National Laboratory sea ice model Version 4 (CICE4), the Parallel Ocean Program Version 2.0 (POP2.0), and the Marine Science and NUmerical Modeling (MASNUM) surface wave model; see Qiao et al. (2016) for the related references in more detail. The horizontal resolutions of the model are as follows: T42 (about 2.875°) for both CAM3 (with 26 vertical layers) and CLM3.5, a nominal 1° (about 1.125° in longitude and 0.27°−0.54° in latitude) with the North Pole displaced to Greenland for POP2.0 (with 40 vertical layers) and CICE4, and 2° for the MASNUM surface wave model (with a resolution of 30° for wave direction). The detailed descriptions of the FIO-ESM v1.0 can be found in Qiao et al. (2013).

Previously, the role of ocean surface waves has not been explicitly considered in climate and earth system models. It is the first time that in 2013 the nonbreaking surface wave vertical mixing effects are adequately incorporated in the FIO-ESM v1.0 for climate modeling (Qiao et al., 2013). Comparisons with and without the wave effects indicate that climate simulations in the FIO-ESM v1.0 have been reasonably improved, such as the sub-tropical mixed layer depth, heat content in the upper ocean, and so on (Chen et al., 2018, 2019). Moreover, surface wave effects act to significantly reduce tropical biases that have been long standing in all climate models, including too cold simulations of annual mean SST field and its seasonal cycle in the eastern equatorial Pacific (Song et al., 2012, 2014).

In terms of ENSO simulations, the FIO-ESM v1.0 exhibits good performance, including ENSO irregularity, spatial distributions of SST variability, and ENSO prediction skill (Song et al., 2015). For example, the wavelet analysis of Niño 3.4 index shows that the FIO-ESM v1.0 can depict a broad spectral peak in the range of 2-7 years with double peaks around 3 and 5 years (Fig. 1b), which is basically consistent with what is observed in the HadISST (Fig. 1a). However, there are obvious discrepancies compared with observations, which is demonstrated to be related to the simulated biases of SST mean state in the tropical region (Chen et al., 2019). For example, the simulated amplitude of SST variability (about 2℃) is much stronger than the observed (about 1℃). In addition, a spurious amplitude peak emerges in boreal summer, although the largest SST anomalies occur in boreal winter during ENSO evolution (Fig. 8a & b). The FIO-ESM v1.0 can capture the basic characteristics of ENSO, but still suffers from common simulation biases that are similar to other CMIP5 simulations. Therefore, the model needs to be improved to reduce these simulation biases.

As previous studies indicated, simulation skills of the OGCMs can be substantially improved by including non-breaking wave-induced vertical mixing. Therefore, three more distinctive physical processes-related air-sea interactions were considered; a new FIO-ESM v2.0 version was configured by including the effects of the Stokes drift on momentum and heat fluxes, the sea spray on air-sea heat flux and the SST diurnal cycle, respectively. Besides incorporating these new physical processes, the model components and resolutions were also upgraded relative to the FIO-ESM v1.0. For example, the AGCM of the FIO-ESM v2.0 is the CAM5, but with the Finite-Volume (FV) dynamical core; its horizontal resolution is 0.9°×1.25°, with 30 vertical layers, which has a much higher resolution than that of the FIO-ESM v1.0. The land surface model is upgraded from CLM3.5 to CLM4.0 with the same horizontal resolution as in the CAM5. The OGCM is the POP2 with the horizontal resolution of 1.1°×(0.27–0.54°), which is the same as that in the FIO-ESM v1.0, but the vertical layers are increased from 40 to 61 layers with the first layer being defined at the sea surface (0 m) for the SST diagnoses that are based on the parameterization of SST diurnal cycle (Yang et al., 2017). The ocean surface wave model is the MASNUM surface wave model (Qiao et al., 2016), with its horizontal resolution increasing from 2° to 1.1°×(0.27–0.54°), which is the same as that of the POP2 and CICE4. The coupling frequency is also increased for each component model. The atmosphere, land surface, and sea ice components exchange information with the coupler every 30 min, while the ocean and wave models exchange information with the coupler at 3 h intervals. The details of the FIOESM v2.0 can be found in Song et al. (2019). Currently, FIO-ESM v2.0 has participated in the CMIP6.

Compared with the FIO-ESM v1.0, the first preliminary inspection indicates that the simulation biases in the FIO-ESM v2.0 are effectively reduced, especially in the tropical region, where the SST biases can be reduced by half relative to the FIO-ESM v1.0 (Song et al., 2019). Moreover, the simulation skills of ENSO events are also greatly improved. As shown in Fig. 8c, the FIO-ESM v2.0 not only can capture the broad spectral peak in the range of 2–7 years, but also the ENSO irregularity and amplitude that are much more reasonable than those in the FIO-ESM v1.0. For instance, the simulated amplitude of SST variability decreases from 2℃ in the FIO-ESM v1.0 to 1.2℃ in the FIO-ESM v2.0, a definite improvement compared with the observed SST variability (which is 1.0℃). In particular, the spurious peak of SST variability in boreal summer is eliminated in the FIO-ESM v2.0. These preliminary results are encouraging for the improvements brought by incorporating the surface wave processes and SST diurnal cycle in the air-sea interaction in the CGCM. The related dynamical reasons for the improvements in terms of ENSO simulations in the FIO-ESM v2.0 need to be investigated in the future.

The FIO-ESM v1.0 has been further used for ENSO prediction experiments. Using an ensemble adjusted Kalman filter assimilation scheme, an operational short-term climate prediction system has been set up to predict ENSO events with ten ensemble members (Song et al., 2015). The hindcast results for the period of 1992–2014 showed that this system has high predictability for ENSO events (Fig. 9). For instance, the correlation coefficients/root mean square errors (RMSEs) of Niñ o-3.4 index between hindcast results and observations respectively are 0.96/0.23℃, 0.87/0.46℃, and 0.77/0.60℃ at the 1-, 3-, and 6-month lead times (Song et al., 2015), which is comparable to other dynamical CGCMs. Moreover, this climate prediction system also shows good skills in real-time ENSO prediction (e.g., 2015/16 El Niño event); the prediction results are routinely provided to the National Marine Environmental Forecasting Center (NMEFC) since November 2016. Currently, this climate prediction system is being upgraded to the FIO-ESM v2.0; it is expected that the FIO-ESM v2.0 should have a better performance in terms of ENSO prediction because its simulation skills have a significant improvement compared with the FIOESM v1.0.

The extent to which the FIO-ESM v2.0 has been improved in climate modeling indicates that smallscale ocean surface wave can play important roles in shaping large-scale climate system. Besides the effects of non-breaking wave-induced vertical mixing, the Stokes drift and sea spray, ocean surface waves also affect the climate system in combinations with other processes, such as the changes in the atmospheric bottom friction and albedo through the sea surface roughness, the atmospheric boundary layer aerosol through wave breaking, and the sea ice distribution through the surface wave-sea ice interactions (Chen et al., 2019). Therefore, the roles of ocean surface waves in the climate modeling should be paid more attention for better simulation and prediction of the climate and particularly ENSO.

5.3 CGCM-based activity at CMA

One CGCM was formulated at the Beijing Climate Center (BCC) in China Meteorological Administration (CMA), and the model has been continually improved for climate simulations, including ENSO. Started in the 1990s, the first version (BCC-CGCM1.0) was built up, including an AGCM and OGCM, with fully coupling without any flux correction, and was used for climate simulation and prediction (Ding et al., 2004). However, the performance of ENSO simulation with the BCC-CGCM1.0 was not good enough and further fundamental improvements have been pursued. Then, the second generation of BCC climate system models, named as BCC-CSM1.1 with a coarse resolution in the atmosphere (approximately 2.812 5°×2.812 5°) and BCC-CSM1.1m with a medium resolution (approximately 1.125°×1.125°) was developed (Wu et al., 2010, 2014). The BCC-CSM1.1 and BCCCSM1.1m are fully coupled global climate-carbon models, including the atmospheric component BCCAGCM2.2 (the Beijing Climate Center Atmosphere General Circulation Model version 2.2; Wu et al., 2010), land component BCC-AVIM1.0 (Wu et al., 2014), oceanic component MOM4-L40, and sea ice component. These components interact with each other through fluxes of momentum, energy, water, and carbon at their interfaces. Both the oceanic and sea ice components were developed by GFDL (Griffies et al., 2008), using a tripolar grid with which the latitudinal resolution is 1° longitude and the meridional resolution ranges from 1/3° latitude between 10°S and 10°N to 1° latitude at 30°S/30°N poleward.

It turns out that ENSO can be well reproduced in BCC-CSM1.1m, including the main features such as ENSO intensity, spatial pattern, feedbacks, and teleconnections (Ferrett et al., 2020), as well as some particular properties such as the ENSO combination mode (Ren et al., 2016), and ENSO persistence barrier (Tian et al., 2019). However, there are model biases which are common to other climate models. For example, BCC-CSM1.1m has a simulated ENSO scale that is much shorter than the observed (Lu et al., 2018); an improving effort is made later to significantly reduce this bias by enhancing the convection entrainment process in the atmosphere model, leading to ENSO periods that are elongated (Lu and Ren, 2016, 2019). In addition, BCC-CSM1.1m has difficulty in distinguishing the two types of El Niño (Ren et al., 2019), although this inability can be empirically corrected with statistical calibrations (Wang et al., 2020). Since 2015, this model has been applied to operational ENSO prediction and routine products are issued in BCC and IRI/CPC every month (Fig. 3; Ren et al., 2016). Overall, BCC-CSM1.1m has shown its good ability in predicting the 2015–2016 super El Niño event, as represented in the China multi-model ensemble system (Ren et al., 2017, 2019). Further improvements are clearly needed for BCC-CSM1.1m for reducing simulation biases of SST and tropical rainfall associated with ENSO.

In recent years, BCC-CSM1.1m has been updated to the newly developed version of BCC-CSM2-MR (Wu et al., 2019), which has been used to conduct most experiments of the CMIP6 in BCC (Xin et al., 2019). The atmospheric component in BCC-CSM2- MR is updated to BCC-AGCM3-MR (Wu et al., 2019) at a horizontal resolution of T106 and 46 vertical layers with its top at 1.459 hPa, and the land component is the newly improved version BCCAVIM2.0 (Li et al., 2019). The ocean and sea ice components in BCC-CSM2-MR are the same as those in BCC-CSM1.1m. According to the CMIP historical simulations from both models of BCC-CSM1.1m and BCC-CSM2-MR, the simulated SST anomaly standard deviation in the ENSO key regions from BCC-CSM2-MR is weaker than that in BCCCSM1.1m. Therefore, the amplitude of the ENSO simulated by BCC-CSM2-MR is closer to the observation than that simulated by BCC-CSM1.1m (Wu et al., 2019). BCC-CSM2-MR improves the positive feedback of the thermocline and zonal advection, which result from stronger air–sea interaction and the amplification of the effect of atmospheric Bjerknes feedback part. Through coupling, the Niño3 SST anomaly rapidly results in a zonal wind stress anomaly near the Niño4 region and further leads to the change in the depth of the thermocline in the equatorial eastern Pacific. Meanwhile, the Niño3 SST anomaly also results in a zonal current anomaly in the equatorial eastern Pacific, which leads to zonal advection of the mean change in SST.

6 DISCUSSION AND CONCLUSION

ENSO is the strongest interannual variability in the tropical Pacific, and yet has been identified as the most predictable interannual signal in the climate system. Extensive and intensive studies for ENSO have been conducted in the past, reaching a stage where coupled models can now be routinely used for real-time ENSO predictions. Indeed, great achievements in ENSO studies (from its discovery, to process understanding and to real-time prediction) have been recognized as one of the most successful examples for scientific research and applications in Earth sciences (McPhaden et al., 2006). Such progresses can be attributed to several integrated efforts. International observing networks for the coupled ocean-atmosphere states in the tropical Pacific have been developed and implemented in a real-time manner, with in situ and satellite-based measurements such as TOGA, Global Ocean Observation System (GOOS) and Argo, etc. Physically, ENSO processes have been understood very well.

As one of the most powerful tools for ENSO studies, various coupled models have been developed to represent ocean-atmosphere interactions in the tropical Pacific, differing in dynamical formulations and process parameterizations. In terms of complexity, CGCMs consist of complicated AGCMs and OGCMs, which are primitive equation-based and include comprehensive processes. Meanwhile, simplified coupled models have been formulated for uses in ENSO modeling because they can well depict characteristic atmospheric and oceanic responses to ENSO, including linear two-layer representations in the vertical and shallow water equation model, etc. For instance, in simplified ICMs and HCMs, one component is taken as a simple anomaly form (i.e., anomaly fields are directly produced), whereas the other component is taken to be AGCM or OGCM; then, an oceanatmosphere anomaly coupling is adopted, in which observed climatological fields are prescribed for use in the exchange of fluxes at the ocean-atmosphere interface; this procedure acts as a flux correction, which helps alleviate climate drift problem. When these individual oceanic and atmospheric models are coupled, great cares are needed for the coupled modeling so that ENSO cycles can be self-sustained by the model itself resulting from coupled interactions between the atmosphere and ocean in the tropical Pacific. Indeed, intensive and extensive efforts have been devoted to ENSO model developments and improvements in the past several decades. Now, various models with different degrees of complexity have been available and widely used in ENSO modeling, which is strongly model dependent. There are advantages and disadvantages of differently formulated models. For instance, CGCMs are expensive to run and may experience the so-called climate drift problems due to the fully coupled model configuration. Simplified coupled models (e.g., ICMs and HCMs) are efficient to run and are free from climate drift problems because some observed climatological fields are used in the anomaly coupling, acting as a flux correction. Evidently, simplified models have been widely used not only for computational efficiency in time integration, but also for being amenable for process illustrations and understanding.

Coupled models provide a mathematical basis for ENSO simulations and in the past dedicated efforts have been devoted to ENSO model developments. The model performances for ENSO and climate simulations have been continually improved in the past three decades, including physical parameterizations, fining spatial resolutions, and others. For example, comprehensive CGCM-based simulations without flux corrections can now well depict ENSO evolution and climate variability and change. Furthermore, coupled models have been continually upgraded to the degree that current stateof-the-art coupled models can be used for ENSO predictions. Indeed, the integrated efforts of observing, process understanding and modeling lead to successful implements of the real-time prediction system for ENSO. Importantly, comprehensive CGCMs provide a basis for short-term climate prediction and future climate projection under the influence of global warming as indicated in CMIP6.

Nevertheless, great challenges remain in ENSO modeling, including long-standing model biases that exist in simulations and predictions of ENSO. For example, the available CMIP products still exhibit systematic model biases and uncertainties in ENSO simulations, as indicated in ENSO intercomparison project (Latif et al., 2001). The performance of CGCMs is still not satisfactory in terms of ENSO modeling, with large biases and uncertainties in ENSO simulations and predictions. In terms of representing ENSO behaviors by using current CGCMs, the main biases can be summarized here. First, the simulated ENSO events exhibit a large spread in the amplitude across differently formulated models, ranging from half to more than two times of the observed one; these large uncertainties can be attributed to model dynamical formulations, process representations, and so on. Secondly, more than half of the current CGCMs have difficulty in reasonably representing the ENSO phase locking feature (Ham and Kug, 2014). Thirdly, the simulated ENSO events tend to have a too regular oscillation with a period of near-biennial, rather than the observed period of 2–7 years (AchutaRao and Sperber, 2002; Guilyardi et al., 2004; Capotondi et al., 2015). Fourthly, observations indicate that the amplitude of El Niño events is usually larger than that of La Niña events; yet, the current CGCMs show poor performance in reproducing the asymmetry of ENSO amplitude between the two phases of ENSO, as reported by some recent studies. Fifthly, the current CGCMs have difficulties in capturing the ENSO diversity. As recognized and observed, El Niño can have two flavors in terms of distinct spatial patterns of SST variability: one is the canonical El Niño that exhibits a SST warming center in the eastern equatorial Pacific, and the other is the so-called the central Pacific (CP) El Niño whose warming center tends to occur more dominantly in the central Pacific. Recent studies documented that the current CGCMs still have difficulty in adequately reproducing the two types of El Niño events (e.g., Ham and Kug, 2012; Yeh et al., 2014; Capotondi et al., 2015).

These discrepancies in ENSO simulations using CGCMs limit their applicability to ENSO prediction. Indeed, in terms of ENSO predictions, CGCMs indicate no obvious advantage in prediction skills compared with simplified ICMs and HCMs. Bias corrections and/or flux corrections are often made to enhance prediction skills when using CGCMs for ENSO real-time predictions. In fact, bias corrections are necessarily needed to extract predictable signals from drifted mean climatological fields. For example, model output statistics (MOS) corrections are often used (e.g., Ji et al., 1996). These problems in ENSO simulation and prediction remain to be addressed for ENSO model developments; clearly, coupled models need to be continually improved.

There are clear rooms for further improvements in ENSO modeling, which has been actively underway. For example, coupled model biases can be of oceanic origin, atmospheric origin and of coupled origin. Previously, various forcings and feedbacks have been identified that can affect ENSO, involving in processes at multiple scales, deterministic and stochastic processes, etc. These processes need to be adequately included to represent their effects on ENSO modulations. In particular, as ocean processes play a critically important role in successfully depicting ENSO, ocean processes need to be adequately represented in ENSO modeling, including ocean vertical mixing, freshwater flux forcing, tropical instability waves, ocean biology-induced heating feedback and so on (Qiao et al., 2004; Zhang et al., 2012; Zhang, 2014; Zhi et al., 2015, 2019; Kang et al., 2017a; Zhu and Zhang, 2018; Gao et al., 2020). Putting all these efforts together, ENSO simulations and predictions can be expected to improve greatly and thus satisfactorily meet various scientific and societal demands for climate studies.

This work provides a brief review of coupled model development and prediction for ENSO studies in China. The coupled models with different complexities reviewed in the paper are limited primarily to those from four institutions only, which have been most active in ENSO modeling in China, with progress being made continually and steadily in history. Indeed, there are other coupled models relevant to ENSO studies in China and thus the current review is incomplete. That is to say, our selections of models for the review are rather objective, but as described above, these reviewed models are most representative of ENSO coupled modeling studies in China because they typically exhibit some unique characteristics and obvious advantages in terms of model dynamical formulations or process representations. Note that this short review has several weaknesses that need to be taken care. For example, it is necessary to compare the strengths and weaknesses of ENSO modeling studies in China to those of international studies so that useful information can be provided about what is the current level of ENSO studies in China. Here, we just show the performances of the objectively selected models in different aspects. It is highly desirable to provide a synthetic view of biases and uncertainties of the reviewed models in specific aspects of ENSO properties, including the amplitude and periodicity, the seasonal phase locking, asymmetry and irregularity, spatial-temporal evolution and prediction skill. In doing so, common measures with the same criterion and metrics should be applied to all the models reviewed in order to quantify the model performances in a coherent and uniform way. In other words, a parallel analysis should be performed to provide information about some important intermodal differences in ENSO simulations and predictions. Fortunately, recently released CMIP6 products have been available, which can be systematically compared with each other for such purposes. Some related detailed analyses and modeling experiments are underway and the results will be presented elsewhere (e.g., Zhu et al., 2020).

7 DATA AVAILABILITY STATEMENT

Except for Fig. 3, all figures and tables in the paper are created by the authors (the figures plotted by using the Grid Analysis and Display System (GrADS) which is available at http://www.iges.org/grads/grads.html). The data and computer codes used in the paper are available from the authors (rzhang@qdio.ac.cn).

8 ACKNOWLEDGMENT

The authors wish to thank the three anonymous reviewers for their comments that helped to improve the original manuscript.