A review of progress in coupled ocean-atmosphere model developments for ENSO studies in China*

ZHANG Rong-Hua , YU Yongqiang , SONG Zhenya , REN Hong-Li 1 TANG Youmin 1 1 QIAO Fangli , WU Tongwen 1 GAO Chuan HU Junya TIAN Feng ZHU Yuchao CHEN Lin , 1 LIU Hailong , LIN Pengfei , WU Fanghua WANG Lin

1 Key Laboratory of Ocean Circulation and Waves, Institute of Oceanology, and Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China

2 Laboratory for Ocean and Climate Dynamics, Pilot National Laboratory for Marine Science and Technology, Qingdao 266237, China

3 Center for Excellence in Quaternary Science and Global Change, Chinese Academy of Sciences, Xi’an 710061, China

4 University of Chinese Academy of Sciences, Beijing 100049, China

5 LASG, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100049, China

6 First Institute of Oceanography, Ministry of Natural Resources, Qingdao 266061, China

7 Laboratory for Regional Oceanography and Numerical Modeling, Pilot National Laboratory for Marine Science and Technology, Qingdao 266237, China

8 Key Laboratory of Marine Science and Numerical Modeling (MASNUM), Ministry of Natural Resources, Qingdao 266061, China

9 National Engineering Laboratory for Integrated Aero-Space-Ground-Ocean Big Data Application Technology, Qingdao 266061, China

10 State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, Beijing 100081, China

11 Laboratory for Climate Studies and CMA-NJU Joint Laboratory for Climate Prediction Studies, National Climate Center, China Meteorological Administration, Beijing 100081, China

12 Environmental Science and Engineering, University of Northern British Columbia, British Columbia V2N 4Z9, Canada

13 State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, Ministry of Natural Resources, Hangzhou 310012, China

14 Key Laboratory of Meteorological Disaster, Ministry of Education (KLME), and Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters (CIC-FEMD), Nanjing University of Information Science and Technology, Nanjing 210044, China

15 Beijing Climate Center, China Meteorological Administration, Beijing 100081, China

Abstract El Niño-Southern Oscillation (ENSO) is the strongest interannual signal that is produced by basinscale processes in the tropical Pacifi c, with signifi cant effects on weather and climate worldwide. In the past, extensive and intensive international efforts have been devoted to coupled model developments for ENSO studies. A hierarchy of coupled ocean-atmosphere models has been formulated; in terms of their complexity, they can be categorized into intermediate coupled models (ICMs), hybrid coupled models (HCMs), and fully coupled general circulation models (CGCMs). ENSO modeling has made signifi cant progress over the past decades, reaching a stage where coupled models can now be used to successfully predict ENSO events 6 months to one year in advance. Meanwhile, ENSO exhibits great diversity and complexity as observed in nature, which still cannot be adequately captured by current state-of-the-art coupled models, presenting a challenge to ENSO modeling. We primarily reviewed the long-term efforts in ENSO modeling continually and steadily made at different institutions in China; some selected representative examples are presented here to review the current status of ENSO model developments and applications, which have been actively pursued with noticeable progress being made recently. As ENSO simulations are very sensitive to model formulations and process representations etc., dedicated efforts have been devoted to ENSO model developments and improvements. Now, different ocean-atmosphere coupled models have been available in China, which exhibit good model performances and have already had a variety of applications to climate modeling, including the Coupled Model Intercomparison Project Phase 6 (CMIP6). Nevertheless, large biases and uncertainties still exist in ENSO simulations and predictions, and there are clear rooms for their improvements, which are still an active area of researches and applications. Here, model performances of ENSO simulations are assessed in terms of advantages and disadvantages with these differently formulated coupled models, pinpointing to the areas where they need to be further improved for ENSO studies. These analyses provide valuable guidance for future improvements in ENSO simulations and predictions.

Keyword: El Niño-Southern Oscillation (ENSO); coupled ocean-atmosphere models; simulations and predictions; model biases and uncertainties

1 INTRODUCTION

El Niño is referred to anomalous warming of sea surface temperature (SST) in the eastern-central equatorial Pacifi c Ocean, and the Southern Oscillation (SO) is referred to the differences in sea level pressure (SLP) across the tropical Pacifi c 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 Pacifi c, ENSO signifi cantly affects weather and climate all over the world.

In the past several decades, ENSO has been extensively and intensively investigated; signifi cant 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 fi elds (e.g., SLP) were identifi ed 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 Pacifi c (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 Pacifi c (the so-called thermocline feedback or Bjerknes feedback; Bjerknes, 1969; Wyrtki, 1975). In the 1970s and 1980s, highly idealized conceptual models and simplifi ed 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 Pacifi c 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 Pacifi c 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 identifi ed 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 classifi ed into highly idealized conceptual models, simplifi ed 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 Pacifi c 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 Pacifi c; 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 Pacifi c. 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 Pacifi c in early 2014, a pause in the warming in mid and late 2014, and then amplifi ed 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 identifi ed 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, specifi c 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 Pacifi c, which embarrassed the scientifi c 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, specifi c 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 identifi ed 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. Specifi cally, 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, fi ning 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 simplifi ed 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 Pacifi c. For example, statistical models, constructed from historical data using statistical methods, have been successfully applied to the coupled modeling of ENSO. Another justifi cation 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, simplifi ed atmospheric dynamical models have been widely used for representing atmospheric characteristic responses to ENSO-induced SST variability in the tropical Pacifi c, 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 simplifi ed 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 simplifi cations 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-defi ned 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 simplifi ed 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 fi elds (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 fi elds (i.e., SST), atmosphere models can depict surface wind stress fi elds, etc. in association with the SO. As ENSO originates from the ocean-atmosphere coupling in the tropical Pacifi c, 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 Pacifi c. Various factors have been identifi ed 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-defi ned 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 Pacifi c.

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 simplifi ed models, in which a kind of simplifi ed anomaly model is adopted for the atmospheric or oceanic component, respectively. ICMs are formed by coupling both simplifi ed 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 fi rst dynamical model used for ENSO prediction experiments (Cane and Zebiak, 1985; Cane et al., 1986). In a more complicated confi guration, HCMs are also simplifi ed 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 simplifi ed 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 simplifi ed 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, simplifi ed 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-defi ned relationships among oceanic and atmospheric fi elds, 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 simplifi ed 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 fi elds of the atmosphere and ocean are calculated directly, whereas seasonally varying climatological fi elds are specifi ed 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 simplifi ed 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 fi elds 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, simplifi ed 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 simplifi ed dynamical ocean model.

These simplifi ed 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 ( Te) 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 fi rst dynamical model that is used for ENSO prediction (Cane et al., 1986) and has been pioneering in ENSO studies. More specifi cally, the ZC87 adopted a two-layer linear shallow water approximation to represent the ocean-atmosphere system in the tropical Pacifi c. 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 ( Te) anomaly fi eld 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 simplifi ed 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 simplifi ed 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 Pacifi c 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 Pacifi c (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 stratifi cation 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 ( Te) is parameterized. An optimized procedure is developed to depict Teusing 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 Pacifi c. The ICM has its big computational advantage relative to HCMs and CGCMs, and yet the performance of ENSO simulations is satisfactorily good (Fig.2).

Fig.1 Schematic diagrams showing the structures of the IOCAS ICM and HCM AGCM

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 Pacifi c (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 Pacifi c 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, simplifi ed ICMs have played a very important role in the developments of ENSO understanding and predictions. Such a level of simplifi ed 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 simplifi ed ICMs with their anomaly formulation have limited applicability, particularly in such a situation where both mean climatology and interannual anomaly fi elds 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 fi elds 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

Fig.2 Time series of SST anomalies (°C) averaged over the Niño3.4 region (the left panels of (a)–(c)) and the standard deviation of SST interannual variability (°C) (the right panels of (a)–(c)) for observation (a), ICM (b), and HCM (c), respectively; power spectra of the Niño3.4 SST index in observation and model simulations (d); seasonal variations of the standard deviation of the Niño3.4 SST index (°C) as a function of calendar month (e)

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, HCMAGCMis referred to as an HCM in which an AGCM is coupled with a simplifi ed dynamical ocean model; HCMOGCMis referred to as an HCM in which an OGCM is coupled with a simplifi ed atmospheric model (statistical or dynamical). Because one component of such an HCM is taken as an anomaly form in which only interannual perturbation fi elds 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 fi elds are prescribed for use. For the HCMAGCM, the AGCM produces total surface wind fi elds, whose interannual anomalies are calculated relative to their corresponding climatological fi elds 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 HCMOGCM, the OGCM produces total SST fi elds, whose interannual anomalies are calculated relative to the corresponding climatological SST fi eld 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 simplifi ed HCM is computationally efficient to run in time integration.

Fig.3 The predictions (colored lines) of the Niño 3.4 SST anomalies for the period 2020-21, which are made from mid-May 2020 as initial condition using different coupled models, including the IOCAS ICM and BCC-CSM1.1m

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 fi elds are explicitly separated from total fi elds, 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 HCMOGCM, 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 simplifi ed 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 simplifi ed 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 HCMAGCM, in which a simplifi ed 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 fi fth 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 (SSTclim) to form total SST fi elds, which are used to force the AGCM; the AGCM produces total wind stress ( τ) fi elds, 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 HCMAGCMis also climate drift-free because of using anomaly coupling, with which observed seasonal climatology of SST is prescribed from observations and only oceanic perturbation fi elds are calculated directly from the ocean model; this confi guration 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 HCMAGCMconfi guration, which allows us to focus on primary processes that are only important to SST interannual variability in the tropical Pacifi c, without the explicit need for caring much of those responsible for the mean climatology in the ocean.

Currently, a regional HCMAGCMis formulated, in which active ocean-atmosphere coupling is allowed only in the tropical Pacifi c although a global AGCM is used, which can depict global atmospheric responses to SST anomalies in the tropical Pacifi c. A long-term simulation (more than 100 years) has been performed using the HCMAGCM; its simulation skills for ENSO are evaluated using SST over the tropical Pacifi c (Fig.2). The HCMAGCMproduces 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 Pacifi c in the HCMAGCMcompared 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 HCMAGCMallows 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 Pacifi c, including the Pacifi c and North America (PNA) patterns. These simulated atmospheric response patterns in the HCMAGCMare similar to what is observed in nature. Therefore, this HCMAGCMcan not only depict the ENSO-related interactions over the tropical Pacifi c, but also reproduce ENSO-induced global atmospheric variability, thereby providing a useful modeling tool for studying ENSO-related coupling in the tropical Pacifi c 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 Pacifi c, consisting of an OGCM and a simple statistical model for the atmosphere (Fig.4), denoted as HCMOGCM(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 fi xed to be constant). Chen et al. (1994) embedded a bulk ML model into the ocean model, in which the depth of the ML ( Hm) 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 (SSTinter=SST–SSTclim); 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 fi elds 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.

Fig.4 A schematic diagram illustrating a hybrid coupled model (HCM OGCM), including the atmosphere and ocean physics and biology (AOPB) in the tropical Pacifi c

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 fi rst 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 confi guration 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 QBon 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 QB) on ocean and coupled simulations have been examined (Kang et al., 2017b; Gao et al., 2020)

Fig.5 Longitude-time sections along the equator for interannual anomalies of SST (a), zonal wind stress (b), and Chl (c) simulated from the HCM OGCM

The use of the comprehensive OGCM for the HCMOGCMallows 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 Pacifi c (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 HCMOGCMhas 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 fi elds from the HCM AOPB (Tian et al., 2018; Zhang et al., 2018). As evident from the HCMOGCM, 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 Pacifi c; so, such a model confi guration does not permit the investigation of global atmospheric teleconnection that can be affected by SST variability in the tropical Pacifi c. 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 fi elds 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 fi elds 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 simplifi ed 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 confi guration 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 Pacifi c, 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 stratifi cation in the vertical, total available energy consideration, and so on.

For example, the IAP OGCM is the fi rst 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 fi xed, 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 fi ltered 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 artifi cially 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 artifi cial 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 artifi cially 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 fi nite-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 fi rst, 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 Pacifi c 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 Pacifi c, 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 Ocean-Atmosphere-Land System Model (FGOALS) CGCMs has been developed for climate modeling at LASG/IAP. Specifi cally, three generations of the FGOALS CGCMs have been developed at LASG/IAP since 2002. The fi rst 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 Pacifi c. 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 fi rst generation of FGOALS (i.e., as indicated by g1 in Fig.6) shows severely overestimated interannual variability over the central and eastern equatorial Pacifi c (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 stratifi cation 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 fi ve of FGOALS CGCMs capture the spatial distribution features of ENSO as observed. The specifi c standard deviations of the Niño3 index and Niño4 index derived from the observation and the fi ve 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 fi ve 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 fi ve 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 defi ciency in ENSO periodicity simulation has been alleviated in the other three versions (i.e., s2, g2, and g3) to some extent.

Fig.6 Spatial distributions of the standard deviation for SST interannual anomalies (°C), derived from observation (a) and the fi ve simulations using different versions of the FGOALS CGCMs for g1 (b), s2 (c), g2 (d), f3-L (e), and g3 (f), respectively

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 fi ve 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 theunderestimation 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.

Table 1 The standard deviation of Niño SST indices and the skewness, estimated from observations and simulations using different versions of the Flexible Global Ocean-Atmosphere-Land System Model (FGOALS) CGCMs developed at LASG/IAP

Fig.7 Seasonal variations in the standard deviation of Niño3.4 index as a function of calendar months, derived from the observation and the fi ve simulations using different versions of the FGOALS CGCMs for g1, s2, g2, f3-L, and g3, respectively (a); the power spectra calculated for observation and the fi ve simulations (b)

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). Specifi cally, 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 fi rst 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 signifi cantly reduce tropical biases that have been long standing in all climate models, including too cold simulations of annual mean SST fi eld and its seasonal cycle in the eastern equatorial Pacifi c (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ñ o3.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°C) is much stronger than the observed (about 1°C). 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.

Fig.8 Wavelet analyses of historical Niñ o3.4 SST (5°N−5°S, 170°W−120°W) anomalies from HadISST observation (during the period 1870-2014) (a) and simulations for FIO-ESM v1.0 (during the period 1870-2005) (b), and FIO-ESM v2.0 (during the period 1870-2014) (c), respectively

Fig.8 Continued

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 confi gured 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 fi rst layer being defi ned 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 fi rst 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 thebroad 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°C in the FIO-ESM v1.0 to 1.2°C in the FIO-ESM v2.0, a defi nite improvement compared with the observed SST variability (which is 1.0°C). 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.

Fig.8 Continued

The FIO-ESM v1.0 has been further used for ENSO prediction experiments. Using an ensemble adjusted Kalman fi lter 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°C, 0.87/0.46°C, and 0.77/0.60°C 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 signifi cant 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 fi rst 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 signifi cantly 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 amplifi cation 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 Pacifi c. Meanwhile, the Niño3 SST anomaly also results in a zonal current anomaly in the equatorial eastern Pacifi c, which leads to zonal advection of the mean change in SST.

6 DISCUSSION AND CONCLUSION

ENSO is the strongest interannual variability in the tropical Pacifi c, and yet has been identifi ed 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 scientifi c 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 Pacifi c 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 Pacifi c, 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, simplifi ed 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 simplifi ed ICMs and HCMs, one component is taken as a simple anomaly form (i.e., anomaly fi elds 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 fi elds 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 Pacifi c. 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 confi guration. Simplifi ed coupled models (e.g., ICMs and HCMs) are efficient to run and are free from climate drift problems because some observed climatological fi elds are used in the anomaly coupling, acting as a flux correction. Evidently, simplifi ed 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, fi ning 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 Pacifi c, and the other is the so-called the central Pacifi c (CP) El Niño whose warming center tends to occur more dominantly in the central Pacifi c. 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 simplifi ed 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 fi elds. 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 identifi ed 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 scientifi c 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 specifi c 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 fi gures and tables in the paper are created by the authors (the fi gures 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.