Simulation of Indian summer monsoon rainfall, interannual variability and teleconnections: evaluation of CMIP6 models

We analyse the performance of global climate models of 6th\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{th}$$\end{document} generation of Coupled Model Intercomparison Project (CMIP6) in simulating climatological summer monsoon rainfall over India, interannual variability (IAV) of all-India summer monsoon rainfall (ISMR) and its teleconnections with rainfall variability over equatorial Pacific and Indian Oceans. The multimodel ensemble mean (MME) of 61 CMIP6 models shows the best skill in simulating mean monsoon rainfall over India compared to the MMEs of 6th\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{th}$$\end{document} generation atmosphere-only models (AMIP6) and the previous generations of Atmospheric and Coupled Model Intercomparison Projects (AMIPs and CMIPs). Systematic improvement and reduction in bias are evident from lower to higher AMIPs/CMIPs. Still, there exists dry bias over a narrow region of the monsoon zone of central India besides wet and cold bias over the surrounding oceans. The persistence of errors in atmosphere-only models hints that the source of errors could be with atmosphere models. Fifteen CMIP6 models selected through objective criteria, perform the best in simulating mean monsoon, IAV of ISMR, the strong inverse relationship between ISMR and Boreal summer El Niño-Southern Oscillation (ENSO), and the inverse relationship between all-India rainfall and north–west tropical Pacific rainfall in June. Several models reproduce the dipole structure of Equatorial Indian Ocean Oscillation (EQUINOO) with the centres over western and eastern equatorial Indian Ocean. But, ISMR-EQUINOO relationship in many of them is opposite to the observed. Our analysis implies the need for capturing ISMR-EQUINOO link to improve the simulation of IAV of ISMR which is crucial for reliable monsoon prediction and projection.


Introduction
India receives more than 75% of its annual rainfall during the Boreal summer monsoon season of June to September (JJAS). The Indian summer monsoon (ISM) rainfall is the major lifeline of the people of India which constitutes roughly one-sixth of the world's population. The ISM is a robust feature of the annual cycle but all-India summer monsoon rainfall (ISMR) manifests year-to-year or interannual variation (IAV) with a standard deviation of 9.7% of its climatological mean of ∼837 mm during 1951-2010. Sufficient ISMR is essential for bountiful agriculture and food security of India. Parthasarathy et al. (1988) reported 10-20% increase of food grain production over India, associated with positive anomalies of ISMR. Thus, IAV of monsoon rainfall is of great consequence for the socio-economic sector of India and understanding, simulating and predicting ISMR and its variability are critically important. Besides, the capability to simulate summer monsoon rainfall over India is highly essential for climate models to be useful for generating reliable climate change projection of monsoon which in turn is critical for future planning and development.
In the beginning of the modelling era, simulations of Asian summer monsoon using General Circulation Models (GCMs) addressed aspects such as the seasonal variation of tropical circulation, role of mountains, and the mean features (Manabe et al. 1974;Hahn and Manabe 1975, etc.). In a comparison of the monsoon simulations by a set of models, Gilchrist (1977) showed that the models, in general, 1 3 simulate the monsoon circulation with reasonable fidelity, but fail to simulate more detailed characteristics of rainfall distribution. It was found that the simulations are sensitive to parameterizations and most of the shortcomings appear to be near the orography. Several modelling studies were done on the impact of soil moisture and vegetation (Walker and Rowntree 1977;Shukla and Mintz 1982), and the orographic effects of the Tibetan plateau (Hahn and Manabe 1975;Kitoh and Tokioka 1987;Fennessy et al. 1994).
Another important aspect was the role of tropical sea surface temperature (SST) anomalies in the monsoon. This was studied in detail by Shukla and Fennessy (1994); Ju and Slingo (1995) etc. The lower boundary condition was shown to strongly influence the monsoon simulation (Charney and Shukla 1981;Palmer et al. 1992). As part of the TOGA Monsoon Numerical Experimentation Group (MONEG), contrasting years of ISMR viz. the drought of 1987 and the excess of 1988, were studied using GCMs forced with observed global SST (WCRP 1992;Palmer et al. 1992). Many models showed significant dynamical drift that resulted in rainfall bias but the variability was captured by most models. Many of the state-of-the-art models were capable of realistically simulating the regional and planetary scale circulation associated with ISM, by simulating the climatological features associated with orography as well as the land-sea distribution and associated differential heating, to some extent (e.g., Cess et al. 1993;Gates et al. 1993), but showed a poor simulation of spatial and temporal variability of rainfall over India. The coarse resolution of these models constrained them from correctly representing the complex topography and coastlines as well as sub-grid scale processes. Shukla and Fennessy (1994) demonstrated that the mean state of the model is of utmost importance. The ability of the models in investigating monsoon processes depended largely on their ability to simulate the mean climatology realistically (Webster et al. 1998). Soon, the ability of the models in simulating the monsoon became a critical benchmark to qualify for reliable simulation of tropical climate variability (Sperber and Palmer 1996;Gadgil and Sajani 1999).
Phase two of the MONEG project began the era of Model Intercomparison Projects (MIPs) in which the ability of the models in simulating the observed climate is assessed through a standard set of experiments. When each modelling group does a specific set of experiments, it offers a wide range of possibilities to compare the outputs with observation as well as those of other models. Through this method, a range of model behaviours can be extracted along with identifying their respective strengths and weaknesses. The 10-year integrations of atmospheric general circulation models (AGCMs) forced with observed SST of 1979-1988, were made under the phase-1 of the Atmospheric MIP (AMIP1, Gates 1992). AMIP simulations can be considered to show the potential skill of the models since they are being driven by observed boundary conditions. Sperber and Palmer (1996) followed by Gadgil and Sajani (1999) reported the merits and demerits of AMIP1 models in simulating mean ISMR and IAV of ISMR. Sperber and Palmer (1996) suggested that the models with a better simulation of climatological summer monsoon rainfall are more adept in simulating the IAV. Gadgil and Sajani (1999) emphasized the importance of simulation of mean seasonal variation of the major monsoon rainbelt over the Asia-Pacific region to achieve realistic simulation of IAV once the simulation of mean monsoon rainfall is reasonably close to observation.
Subsequently, the World Climate Research Programme (WCRP) of the Working Group on Coupled Modelling (WGCM) organized the Coupled Model Intercomparison Project (CMIP) in the mid-1990s. It was developed from the AMIP1 (Gates et al. 1999) to coordinate the coupled model simulations worldwide to advance our knowledge of climate variability and change. During the initial generations of CMIP1 and CMIP2, the design was simple with a long control simulation and analyses of different climate models helped in understanding the differences in model response.
With the evolution of CMIP, the suite of experiments also increased with more detailed and elaborate twentieth-century simulations. A major paradigm shift occurred with the introduction of CMIP3 when the model outputs from the state-of-the-art climate models were made accessible to the scientific community at large (Meehl et al. 2007). Reichler and Kim (2008) assessed the performance of three generations of CMIP models (up to 20th century simulations of CMIP3) in simulating the present-day mean climate and found that the latest generation models are much more realistic than their predecessors. The immediate yet less known phase called CMIP4 supplemented the experiments done in CMIP3.
Meanwhile, the major factors behind the IAV of ISMR were identified as the El Niño-Southern Oscillation (ENSO, Rasmusson and Carpenter 1982) over the equatorial Pacific and the Equatorial Indian Ocean Oscillation (EQUINOO, Gadgil et al. 2004). The inverse relation between ISMR and ENSO, with ISMR deficits (excesses) often associated with El Niño (La Niña), has extensively been studied since Walker and Bliss (1932) and Sikka (1980). EQUINOO has been shown to impact ISMR with its positive phase favourable for ISMR (Gadgil et al. 2004). Subsequently, a composite index derived through the bivariate analysis of ENSO and EQUINOO indices against ISMR, was found to explain 54% of ISMR variance and was better correlated with ISMR excess/deficit anomalies than ENSO or EQUINOO alone . Analysis of AGCMs in the CLIVAR project showed that models became successful in predicting the SST and associated rainfall over the El Niño region of the equatorial Pacific, but they failed over the South Asian monsoon region (Wang et al. 2004). Rajeevan and Nanjundiah (2009) analysed the CMIP3 models and found that large dry bias exists in the representation of monsoon rainfall over the continental tropical convergence zone (CTCZ) or the monsoon zone over central India, whereas overestimation of rainfall occurs over the equatorial Indian Ocean. Most of the models were found to correctly simulate the observed ENSO-ISMR relationship but fail to simulate the EQUI-NOO-ISMR relationship.
The next major phase shift happened with CMIP5 (Taylor et al. 2012) which was majorly built upon CMIP3 and included idealized processes and feedback-oriented experiments whose outputs facilitated further understanding of the climate system. Sperber et al. (2013) analysed CMIP3 and CMIP5 models and found that the multimodel ensemble mean (MME) simulation of CMIP5 models outperforms that of CMIP3 models in all diagnostics used in analysing Boreal summer Asian monsoon. However, Shashikanth et al. (2014) highlighted that even with the improvements made concerning the physics of climate, the simulation of ISMR in CMIP5 models was no better than those in CMIP3. Improved monsoon simulation necessitates finer resolution models that can capture orographic effects which are responsible for the large spatial variability of monsoon. Sperber et al. (2013) documented the lack of sufficient skill in even the most sophisticated CMIP3 and CMIP5 models in simulating important aspects of mean ISMR and its variability. Mishra et al. (2018) found that CMIP5 MME captures broad-scale precipitation features over India but with large biases, especially over north India. However, the CMIP5 models were found to simulate the surface air temperature pattern over India with better accuracy than precipitation (Jain et al. 2019). In general, the models failed to capture the contribution ratio of convective to large-scale precipitation which induces biases in the model (Yang et al. 2021).
Very recently, the WCRP has defined phase-6 of the CMIP, the CMIP6 (Meehl et al. 2014;Eyring et al. 2016). An account of the scientific gaps involved in CMIP5 models, how CMIP6 models are expected to fill them, and how CMIP6 is superior compared to CMIP5, is given by Stouffer et al. (2017). At present, more than 60 global modelling centres have contributed to CMIP6 to enhance our understanding of climate and climate change. Over the years, the ability of the climate models to simulate the mean climate including monsoon has increased. However, the models are not perfect like our theoretical understanding of climate which is still incomplete and certain simplifying assumptions are still unavoidable in the models. These imperfections introduce biases into model simulations. Over the years, we have refined our theoretical understanding, improved the physical parameterizations, increased the number and quality of observations, and increased our computational capabilities exponentially. With these advancements, the complexity of climate models has substantially increased so that the current models are far more comprehensive than the models of the 1990s. With the introduction of more powerful supercomputers, the current models are also capable of resolving finer-scale details (e.g., , though the resolution of the CMIP6 models is still not high in comparison. Along with the development in the current generation of models, another complexity is added along because of the introduction of sources of possible errors in the new processes. The uncertainty regarding these processes remains unanswered. Against this background, it is pertinent to ask how much the climate models have improved and how far we can trust the latest generation of coupled models in simulating ISM which is crucial for its prediction and future projection as well. A diagnostic assessment of the simulation of mean and variability of ISM by the latest generation of CMIP6 models is thus relevant and important. In this study, the performance of CMIP6 models is evaluated primarily by comparing them to observation, by taking care of their internal variability. Specifically, we examine if the representation of the major characteristics of mean summer monsoon rainfall over India have improved in CMIP6 models and if so, how it has evolved across the successive generations of the AMIPs and CMIPs since AMIP1. Each generation of MIP is expected to have a range in performance, but it is necessary to understand the general improvement across different phases. For this, we also compare the MMEs of AMIP1, AMIP3, CMIP3, AMIP5, CMIP5, AMIP6 and CMIP6 models. It is also imperative to assess the simulation of IAV of ISMR by CMIP6 models, which is largely modulated by ISM's teleconnection with the interannual modes of variability over the equatorial Pacific and Indian Oceans, i.e., ENSO and EQUINOO. In addition, the variation of all-India rainfall (AIR) in June is observed to have played a decisive role in the IAV of ISMR in recent years since 2007. More importantly, this variation of AIR in June is found to be inversely correlated to the variability of north-west tropical Pacific rainfall (NWTPR) in June . The simulation of this relationship by CMIP6 models is also assessed. Thus, the specific objectives addressed in our performance assessment of CMIP6 coupled general circulation models (CGCMs) in simulating important features of ISM are: (1) what are the improvements in models of CMIP6 if any, in representing climatological mean Indian summer monsoon rainfall compared to its atmosphere-only counterpart, the AMIP6, and the previous generations of AMIPs/CMIPs?, (2) how is the IAV of ISMR represented in CMIP6 models?, (3) how good are the representations of ISMR's teleconnection with ENSO and EQUINOO and (4) a brief look into how the inverse relationship between AIR and NWTPR in June, is simulated by CMIP6 models which is a new aspect investigated in this study. While our study aims to provide an overview and interpretation of the ability of the latest climate models in simulating important aspects of ISM and its variability, analysis of the processes involved is beyond the scope of this study.

MIPs datasets
We use monthly data from historical simulations of 61 CMIP6 (Eyring et al. 2016) models (Table 1) made available by the WCRP and the climate modelling groups (https:// esgf-node. llnl. gov/ search/ cmip6/). The first available realization of each model is used. Historical simulations are long-term simulations up to the present-day, carried out with observationally consistent atmospheric composition such as greenhouse gases (GHGs) and aerosols due to both anthropogenic and volcanic influences, solar forcing and land use. (e.g., Miller et al. 2014). Similarly, the protocol for the AMIP simulations has also evolved since AMIP1, with the emissions of GHGs and other forcing vary as observed in its later generations, instead of fixed values. To verify the skill of CMIP6 models, we also analyse the 6 th generation AMIP simulations of 40 atmosphere-only models forced with prescribed SST since 1979 (AMIP6, Table 2). For the assessment of mean monsoon simulation, we compare available outputs of 43 CMIP5 models (Taylor et al. 2012) and 25 CMIP3 models (Meehl et al. 2007) as listed in Tables 3 and 4 respectively. The atmosphere-only simulations of 30 AMIP5 models (Table 3), 13 available AMIP3 models (Table 4), and 31 AMIP1 models (Table 5) are also compared.
To assess the improvement in the simulation of climatological ISM by CMIP6 models, the MME simulations of AMIP1, AMIP3, AMIP5, AMIP6 and CMIP3, CMIP5, and CMIP6 models are intercompared against the observation. All model datasets are bilinearly regridded to 2 • × 2 • grid for intercomparison and the ISM season is defined as June to September (JJAS). To derive statistics, 1951To derive statistics, -2010To derive statistics, (2000 is taken as the base period for the CMIP6 (CMIP5 and CMIP3) simulations. The base period for the AMIP simulations starts from 1979 up to the respective period available for each generation of AMIP. All MIPs' datasets except AMIP1 are taken from https:// esgf-node. llnl. gov/. AMIP1 datasets are from the study of Gadgil and Sajani (1999).

Validation datasets
Primary observation dataset used for precipitation over India is 0.25 • × 0.25 • gridded rainfall of India Meteorological Department (IMD) provided by Pai et al. (2014), for the period 1951-2019. Daily rainfall estimates from ∼ 6950 rain gauge stations across the country with varying periods of observation were regridded to fixed spatial grids after applying standard quality checks by Pai et al. (2014). The large scale climatological features of rainfall over India derived from this dataset were found to be consistent with those derived from IMD's previous rainfall data sets at 0.5 • × 0.5 • and 1 • × 1 • resolutions (Rajeevan et al. 2008).
The EQUINOO index is defined as the negative of the anomaly of the zonal wind close to the surface (10 m) over the equatorial region (60 • -90 • E; 2.5 • S-2.5 • N) of the Indian Ocean (IO). Surface wind from the fifth-generation European Centre for Medium Range Weather Forecast Reanalysis (ERA5) data (Hersbach et al. 2020) which is regridded to 2.0 • grid for 1979-2019 (https:// cds. clima te. coper nicus. eu/ cdsapp# !/ datas et/ reana lysis-era5-press ure-levelsmonth ly-means? tab= form) is used to estimate the EQUI-NOO index. This is same as the EQWIN index defined by Gadgil et al. (2004). A slightly different region (70 • -90 • E; 5 • S-5 • N) had been used by Saji et al. (1999) for estimating the surface wind index to correlate with Indian Ocean Dipole (IOD) events. However, the results are not sensitive to this small change in the region over which the wind is averaged and the variation of the wind over the two regions is highly correlated with a correlation coefficient of 0.96. The two poles of EQUINOO are taken as western equatorial Indian Ocean (WEIO, 50 • -70 • E; 10 • S-10 • N) and eastern equatorial Indian Ocean (EEIO, 90 • -110 • E; 10 • S-0 • ). The difference between SST anomalies over WEIO and EEIO is taken as the Dipole Mode Index (DMI) which is the index for IOD (Saji et al. 1999). ERA5 data is used for the validation of 850 hPa circulation as well.

Methodology
To assess the skill in simulating the climatological summer monsoon rainfall pattern over the Indian region, Taylor diagram (Taylor 2001) showing the skill scores for the mean spatial pattern simulation ('Taylor statistics') are plotted. These statistics are derived with respect to the observed IMD rainfall over India.Taylor statistics used for assessing the matching of the spatial pattern of JJAS mean climatological rainfall simulated by the models with the observed, are the spatial pattern correlation coefficient (PCC), their centred root-mean-square difference as the percentage of observed standard deviation (RMSE), the ratio of their standard deviations (SD) and the climatological bias (the difference in the climatological mean JJAS rainfall of model minus the observed as the percentage of the observed mean) referred to as BIAS to distinguish from the mean bias at each grid (the difference between climatological mean JJAS rainfall of model and the corresponding observation at each grid-point).
To analyse the IAV of individual models in detail, we have selected a group of models which have better skill in simulating primarily the spatial distribution of JJAS mean rainfall climatology over India. We have also considered their skill in simulating climatological JJAS rainfall over the extended Indian region including the surrounding oceans (68 • E-98 • E, 10 • S-38 • N) and the amount of all-India total climatological JJAS rainfall and temporal coefficient of variation of ISMR (CV in percentage of mean) estimated as: RAIN JJAS is the all-Indian total rainfall during JJAS, is the all-India total climatological JJAS mean and is the standard deviation of all-India total JJAS rainfall during the analysis period (N years).
For assessment of IAV simulation, basic statistics such as temporal mean, SD, CV and simultaneous and time-lagged correlations are used. NWTPR is estimated as the standardized rainfall anomaly over 120 • -150 • E; 20 • -30 • N region. All anomalies of indices are standardized with their standard deviation.

Simulation of mean Indian summer monsoon
First, the performance of the MMEs of CMIP6 and AMIP6 models in simulating the spatial distribution of climatological summer mean (JJAS) rainfall, low-level circulation at 850 hPa and SST over the Asia-Pacific region ( Fig. 1), is analysed. Over the Indian region, intense solar heating in summer results in differential heating between land and surrounding oceans. This leads to the development of crossequatorial southwesterly flow transporting increased moisture flux onto the continent which constitutes the lower limb of the deep overturning meridional circulation associated with the regional Hadley cell. The Himalayan and west coast orographic structures act as barriers, adjacent to which rainfall tends to be maximum. Similarly, climatological seasonal rainfall also tends to be concentrated along the orographic regions of East Asia. The tilted band of concentrated rain (rainbelt) over the seasonal Indian monsoon trough characterised with strong cyclonic circulation, stretches across a wide longitudinal range up to the West Pacific. The large latitudinal extent of high rainfall over Indian longitudes is associated with a secondary rainbelt over the equatorial Indian Ocean. Crossequatorial monsoon flow is maintained by the latent heat release associated with ISM rainfall. The dynamics of ISM can thus be regarded as an extension of the dynamics of the Intertropical Convergence Zone (ITCZ). Thus, the maximum latitudinal extent of the ITCZ occurs over the Indian region. This makes ISM the strongest of global monsoons.
The simulations capture the prominent features of mean monsoon such as the location and orientation of major rainbelts over the Indo-Gangetic plain of central India, the orographic Himalayas, the Western Ghats, the East Asian regions, the maritime continent, South China Sea, the West Pacific and the narrow rain band over the central Pacific slightly north of the equator. The zonally oriented secondary rainbelt over the equatorial Indian Ocean (EIO) has a relative maximum over the east (Fig. 1b). Between the two poles that characterise the IOD (denoted by solid boxes in Fig. 1), climatologically, the eastern EIO (EEIO) is more active with organised convection and cyclonic circulation than the western EIO (WEIO). Both simulations capture the strong cross-equatorial flow over the Arabian Sea, monsoon westerlies, strong trade wind easterlies over the Pacific and the anticyclonic circulation over the northern tropical Pacific. Mostly the rainfall maxima over the oceans are collocated with warm SSTs. The distribution of warm SSTs is well simulated by the MME of CMIP6 but with a shrinking of their area over the warm pool regions of the West Pacific and the Bay of Bengal. There exists warm bias over a narrow region of the tropical Pacific east of 160 • W.

Climatological bias
A deeper look into the differences in simulated climatological summer monsoon climate parameters from the respective observations (referred to as biases) reveals the existence of significant biases (Fig. 2). The wet bias over the southern peninsular India and surrounding oceans, the South China Sea, the warm pool region of the West Pacific and both sides of the equatorial central Pacific, are common in the MMEs of CMIP6 and AMIP6 which are stronger in the MME of AMIP6 especially over the Indian region. These regions show significant cold bias in the MME simulation of CMIP6 which could be associated with the wet bias produced by the atmosphere model. Assessing simulations of 18 pairs of atmosphere-only and coupled models from the CMIP5,  Yang et al. (2019) had reported that slightly improved monsoon precipitation in coupled models is largely due to compensation from SST biases that originate from atmosphere model biases. Over India, the monsoon rainbelt extending northwestward from the head Bay of Bengal along the Indo-Gangetic plain is much weaker in CMIP6 MME. AMIP6 MME also shows a narrow dry bias over this region. The wet bias over the rest of the Indian landmass is much stronger in AMIP6. Correspondingly, there are opposite to the observed (slightly enforcing) lower tropospheric wind bias in the cross-equatorial westerlies over the Arabian Sea in CMIP6 (AMIP6). There is strong dry bias over the oceanic rainbelt of EIO, particularly over EEIO which is more conspicuous in CMIP6 with associated easterly wind bias. Thus, both MMEs show dry bias in simulating Indian monsoon rainbelt and the oceanic rainbelt over EIO. At the same time, both show wet bias over the tropical convergence zones (TCZs) of maritime continent and West Pacific which are collocated with warm SSTs. Another striking feature is that the Meiyu/Changma/ Baiu rainbelt which extends from eastern China to Japan is too weak in the MME of AMIP6. Over oceans, wet (dry) biases are mostly associated with westerly (easterly) wind biases. There is wide-spread cold bias over the oceans except over a part of WEIO and tropical Pacific east of 160 • W.

Improvement in CMIP6 models
Since Knutti (2010) had called out for equal weighting of models in constructing MME, several studies explored the potential of weighting the models by their skill (e.g., Knutti et al. 2017) or grouping them based on their merit in simulating certain important aspect of regional climate (e.g., Jayasankar et al. 2015). So for effective intercomparison, we have selected a group of CMIP6 models with greater skill in simulating ISMR, based on two aspects; i) primarily the spatial pattern of climatological JJAS mean rainfall over the Indian region and ii) the interannual variation of ISMR. The criteria used with regard to the first aspect are: i. PCC for mean summer monsoon rainfall over the Indian continent against the corresponding IMD observation is greater than or equal to 0.75 (i.e., PCC IMD ≥ 0.75), ii. PCC for mean summer monsoon rainfall over the extended Indian region including surrounding oceans (68 • E-98 • E; 10 • S-36 • N) against the corresponding GPCP observation is greater than or equal to 0.65 (i.e., PCC GPCP ≥ 0.65). It is found that the classification is not sensitive to the minor changes in the selected region, iii. SD is between 0.8 and 1.2, iv. The centred pattern RMSE for mean summer monsoon rainfall over the Indian continent is below 85%, v. BIAS in simulating rainfall climatology over the Indian continent is within ± 22.5% of the observed and vi. Climatological JJAS mean rainfall over India is within ± 25% of the observed value.
And the additional criterion used for the second aspect is: vii. The coefficient of variation (CV) is between 8% and 16%.
Based on the aforementioned objective criteria, we have selected 15 CMIP6 models which have better skill in simulating primarily the spatial distribution of JJAS mean rainfall climatology over India and secondarily the IAV of ISMR, for performance assessment in simulating IAV of ISMR. The group of 15 CMIP6 models is hereafter referred to as 'CMIP6_Grp' and the MME of 15 CMIP6_Grp models is referred to as the 'MME_Grp'. Figure 3 shows the comparison of observed and simulated patterns of climatological JJAS mean rainfall over the Indian region from different observations, the MMEs of 61 CMIP6 models and 15 CMIP6_Grp models (MME_Grp) and 15 individual models. Spatial distribution of mean summer monsoon rainfall in GPCCv7 and GPCP datasets is well comparable with IMD rainfall with PCCs of 0.83 and 0.74 respectively. The major rainbelt extends west-northwestward from the head of the Bay of Bengal, across the Indo-Gangetic plain that constitutes the artery of the ISM (Fig. 3a). Variation in the location and intensity of this monsoon rainbelt largely gives rise to the IAV of ISMR. Thus it is important for models to capture this monsoon rainbelt in order to achieve realistic IAV simulation. There is a zonally oriented secondary rainbelt over the equatorial IO between 10 • S-0 • , stretching from 70 • E to 100 • E. The interplay of different processes and the presence of multiple rainbelts or TCZs pose a grand challenge for climate models to simulate the spatial distribution of summer monsoon rainfall over the Indian region. As can be gleaned from Fig. 3, the ISMR simulation is still not perfect in the selected CMIP6 models. ECMWF-IFS-LR, EC-Earth3-Veg, ECMWF-IFS-MR, GFDL-CM4C192 and MIROC6 which have PCCs comparable to that of the MMEs of CMIP6 and CMIP6_Grp (MME_Grp), underestimate the monsoon rainbelt over the continental TCZ with strong dry BIAS. The corresponding simulations of mean monsoon rainfall from the rest of the CMIP6 models are shown in Fig. S1. Several simulations show deficiency in simulating the oceanic TCZ as well. Many models overestimate oceanic rainfall over the equatorial Indian Ocean including EEIO where the observed rainbelt is located. The models show better capability in simulating the orographic rainbelts.
Next, we examine if there is any improvement in simulating mean ISMR pattern and reducing the model bias with the progression across 4 generations of MIPs. Taylor diagram (Taylor 2001) effectively summarises matching of a set of simulated climatological ISMR patterns with the IMD observation by plotting PCC, SD and RMSE. The centred pattern RMSE is equal to the standard deviation of model error. BIAS in model climatology is denoted by an upward (downward) symbol for positive (negative) BIAS classified into 5 categories of magnitudes (represented by symbols of 5 different sizes). Figure 4 shows the Taylor diagram for 61 CMIP6 models, 5 observation datasets and the MMEs of the models of CMIP6, CMIP6_Grp (MME_Grp), AMIP6, CMIP5, AMIP5, CMIP3, AMIP3 and AMIP1. With the improvement in resolution and incorporation of physical processes, the monsoon features get better captured by the models as they progressed from one generation to the next.
Another striking feature is the disparity among the observation datasets or their lack of consensus with IMD observation. This could be due to the diverse algorithms, resolutions and measurement methods of these datasets. It is important to assess the skill of the models against the backdrop of the disparity among 4 land (GPCCv7, UDel, PREC-L and CRU) and the oceanic (GPCP) observations. In this manner, this observational skill estimates help to compare and understand the significance of the corresponding values from the simulations. The model skill is predominantly assessed using mean statistics of the models against the primary observation, the IMD rainfall. The MMEs of most of the MIPs are in very good agreement with IMD statistics. Climate models predominantly show a tendency for negative or dry BIAS in simulating ISMR climatology except for the MME of AMIP6 which shows a wet bias. There is a systematic improvement in the simulation of the spatial distribution of mean monsoon rainfall in terms of statistics from lower to higher generation of MIPs (from AMIP1 to AMIP3/CMIP3 to AMIP5/CMIP5 to AMIP6/CMIP6) which is more evident in the coupled models. Among the MIPs, the MME of CMIP6 models is having the higher skill in terms of PCC and RMSE, although there are models like EC-Earth3, EC-Earth3-Veg and ECMWF-IFS-MR which are slightly better than the MME of all CMIP6 models. The MME of CMIP6_ Grp (MME_Grp) shows the highest skill, as it has the largest PCC, the least RMSE and a realistic standard deviation. The negative BIAS or 'dry' bias is either comparable or lower in Fig. 3 a-d Spatial distribution of climatological JJAS mean rainfall from observations of a IMD, b GPCC v7 and c GPCP, and d the MME of 15 CMIP6 models (MME_Grp) selected as the best group for ISMR simulation (CMIP6_Grp) (left to right in top row). e-t The corresponding distributions from the MME of all 61 CMIP6 models and 15 members of CMIP6_Grp. The pattern correlations PCC IMD ( PCC GPCP ) with IMD (GPCP) rainfall over India (the extended Indian region) for simulations are given as the first (second) value in the lower-left side of each panel from d-t. PCC IMD for GPCC v7 and GPCP are written in the lower-left of b and c ◂ comparison to other MMEs. Also, MME_Grp is better than the best model of CMIP6.
In the 6 th generation, the MME of CMIP6 models becomes slightly better than that of AMIP6 models. Noticeably, the climatological biases of the two MMEs are of opposite sign. To understand if there is a systematic difference between coupled and atmosphere-only models in simulating mean monsoon rainfall over India, we analysed the simulations of 10 pairs of CMIP6_Grp models for which both coupled and atmosphere-only versions are available (Fig. 5). It can be seen that both versions mostly have wet bias over the southern peninsular India which is much stronger in atmosphere-only models (AGCMs). This overestimation of rainfall over the southern peninsular India seems to have contributed to the positive bias in the simulation of monsoon rainfall over India, in the MME of AMIP6 models (also can be seen from Figs. 2 and 4). But, the dry bias over the monsoon zone of central India (monsoon rainbelt), is not reversed in most of the AGCMs. Two models, INM-CM5-0 and MIROC6 show significant wet bias over the monsoon zone which are however associated with strong dry bias over the oceanic rainbelt of EIO, especially over EEIO. Corresponding circulation bias features stronger cross-equatorial winds over the northern Arabian Sea and significant easterly bias over EEIO. Since Sikka and Gadgil (1980), several studies have been pointing out that there is a competition between the oceanic and continental rainbelts over India. However, the persistence of dry bias over the monsoon zone in atmosphere-only models forced with observed SSTs hints that the source of errors could be with the atmosphere model. Figure 6 shows the box and whisker plots of PCC, all-India total summer monsoon rainfall ( ) and its standard deviation ( ) for all MIPs. For each MIP, the ends of the box correspond to the upper and lower quartile values among the PCCs of the member models and the box spans the interquartile range with the median marked by the line inside the box (Fig. 6a). The whiskers outside the box extend from the highest to the lowest PCC among the member models. The corresponding PCC of the MME of each MIP is written along side the box. The minimum, maximum, interquartile Taylor diagram comparing the statistics from 61 CMIP6 models, 5 observation datasets of rainfall (GPccv7, UDel, PREC-L, CRU and GPCPv2 denoted by black filled symbols), the MME of 15 CMIP6 models (MME_Grp) with better skill in simulating ISMR climatology (CMIP6_ Grp) and the MMEs of models of CMIP6, AMIP6, CMIP5, AMIP5, CMIP3, AMIP3 and AMIP1. Statistics plotted are i) pattern correlation coefficient (PCC) given by the azimuthal angle, normalised standard deviation as the radial distance and the normalised standard deviation of errors as semicircles centred around the IMD label. The lower symbol attached to each symbol shows the BIAS in simulating the mean rainfall with respect to IMD observation range (first to the third quartile) and median (second quartile) of PCCs and the PCCs for the MME show clear increase from CMIP3 to CMIP5 to CMIP6. The atmosphere-only counterparts also show improvement beginning from AMIP3 to AMIP6. The range of PCCs remain somewhat similar for AMIP1 and AMIP3. The CMIPs outperform the corresponding AMIPs very clearly in capturing the spatial pattern of climatological summer monsoon rainfall over India. The box and whisker plot for manifest the general tendency for the climate models to underestimate the all-India total monsoon rainfall (Fig. 6b). Dry bias systematically reduces with higher generations of CMIP though it is still stronger than the dry bias in their atmosphere-only counterparts. As explained earlier, the MME of AMIP6 models tends to overestimate , possibly due to the overestimation of rainfall over peninsular India.

Interannual variability of ISMR and monsoon teleconnections in CMIP6 models
It is seen that the dry bias over the monsoon zone of central India is a common problem of climate models which still exists in CMIP6 models. This can have a decisive impact on the simulation of IAV of ISMR. In an apparently corresponding manner, the ISMR standard deviation is also underestimated by the previous generations of models (Fig. 6c), except CMIP6 and AMIP6 models. The variability systematically increased from previous generations of CMIPs till CMIP6. In comparison, AMIPs tend to overestimate the IAV of ISMR. Stronger than observed variability remains as a deficiency in 6th generation MIPs particularly AMIP6. Figure 7 shows against CV from observation, CMIP6 and AMIP6 models. The values of the 15 models of CMIP6_Grp and 10 of their atmosphere-only counterparts are shown in coloured symbols. The CV for the MME gets reduced and the interannual variability gets smoothed out if we average across models. Hence we analyse the simulations of individual models. Majority of CMIP6 models underestimate possibly due to their dry bias over monsoon zone and overestimate the or the IAV as also can be inferred from Fig. 6b and c. Thus, the CV is much larger than the observed. Models of CMIP6_Grp have comparatively lesser bias (Table 6) in simulating not only but the CV as well (Fig. 7). The models such as EC-Earth3-Veg, ECMWF-IFS-LR and MPI-ESM1-2-HR have CV comparable to the observation. INM-CM5-0, MIROC6, NESM3 and GFDL-CM4C192 show good skill in simulating the mean and CV which are close to observation.
The atmosphere-only models (AMIP6) mostly overestimate the mean rainfall and hence the CV is not as high as that of the coupled models. We find that the models still have deficiency in simulating the IAV of ISMR even in the latest generation of MIP. There are more number of coupled models than atmosphere-only models in the two circles representing 15% and 30% deviations from the observations (Fig. 7). When 14 out of 15 of CMIP6_Grp models are clustered around the observations, there are only 4 out of 10 atmosphere-only models inside the circles. This shows the higher skill of CGCMs in simulating the proportion of seasonal mean and IAV compared to the AGCMs. This aspect is seen in the simulations of the 5 th generation MIPs as well (Fig. S2). The coupled (atmosphere-only) models of CMIP5 (AMIP5) tend to underestimate (slightly overestimate) the mean and slightly underestimate (overestimate) IAV compared to observation (Fig. 6). Hence, the CV is not as large as that in CMIP6/AMIP6 models (Fig. S2). As seen in the 6 th generation MIPs, CMIP5 coupled models show better skill with a comparatively larger fraction of them closer to the observation in simulating mean versus IAV, than in the atmosphere-only AMIP5 models. But, the difference in performance between AGCMs and CGCMs is conspicuous in 6 th generation MIPs, possibly contributed by the improvement from CMIP5 to CMIP6.

ISMR-ENSO teleconnection
Several observational (Walker 1924;Troup 1965;Sikka 1980;Shukla and Paolino 1983, etc.) and modelling studies (Shukla and Wallace 1983;Soman and Slingo 1997, etc.) have shown that an out of phase relationship exists between monsoon deficits over India and abnormal rainfall over central-eastern Pacific and vice versa. This suggests that large-scale teleconnection operates through east-west atmospheric circulation. The link between ENSO and ISMR was established in the 1980s with the warm phase of ENSO (El Niño) associated with droughts and its cold phase (La Niña) associated with excess rain in India (Sikka 1980;Ropelewski and Halpert 1987, etc.). Later, KrishnaKumar et al. (2006) showed that El Niño events with warmest SST anomalies in the central equatorial Pacific are more effective in causing subsidence over India than events with warmest SSTs in the eastern equatorial Pacific. The correlation of ISMR with rainfall over Indian and Pacific Oceans shows its teleconnections (Fig. 8a). The largest negative correlations occur over equatorial central-east Pacific regions. Weak negative correlations are seen over a small region of eastern equatorial IO slightly west-northwest of EEIO. The highest positive correlations are seen over the maritime continent of West Pacific, western IO including WEIO and southeast EEIO. We find that CMIP6 models show good skill in simulating the interannual variability of NINO3.4 rainfall. Figure S3 shows the box and whisker plots of the standard deviation of Boreal summer season (JJAS) rainfall over the NINO3.4 region simulated by CMIP6, AMIP6, CMIP5 and AMIP5 models. AMIP6 models simulate the standard deviation closer to the observed, with a smaller spread, possibly because of the prescribed SST boundary forcing over NINO3.4 where the air-sea coupling is strong (e.g., Rajendran et al. 2019). The skills of CMIP5 and AMIP5 models are comparable to those of CMIP6 and AMIP6 models respectively.
The distribution of correlation of ENSO index with tropical rainfall (Fig. 8b) and SST over the oceans (Fig. 8c) resembles Fig. 8a, but with an opposite sign, except over WEIO. El Niño (La Niña) induces suppression (enhancement) of ISMR. Both ISMR and ENSO are positively correlated with WEIO rainfall (Fig. 8a and b) and ENSO is positively correlated with the SST over WEIO (Fig. 8c) as well. At the same time, ENSO is strongly anticorrelated with EEIO rainfall and SST. This suggests that El Niño (La Niña) is associated with warming (cooling) and enhancement (suppression) of rainfall over WEIO and NINO3.4 regions and cooling (warming) and suppression (enhancement) of rainfall over EEIO to west Pacific (Fig. 8b and c). The impact of ENSO on WEIO is reinforcing and on ISMR is suppressing. During ENSO years, its impact can elicit a negative or lack of positive relationship between ISMR and WEIO rainfall (Fig. 8b). Strong or stronger than the observed relationship between ISMR and ENSO in models can also end them up in having the opposite relationship between ISMR and WEIO rainfall if they fail to simulate the basic positive relationship between them as observed.
The basic favourable relationship between WEIO rainfall and ISMR could be due to the fact that enhanced cross-equatorial flow is associated with enhanced off-equatorial diabatic heating along the continental TCZ over the monsoon zone resulting in strong correlation between convection over WEIO and ISMR (negative correlation with OLR as the proxy for convection over WEIO in Fig. 9). The strong correlation between ISMR and WEIO convection could be also aided by the fact that a warmer WEIO allows more moisture to be supplied to the monsoon. But if ENSO induced large-scale subsidence prevails over India, this positive role of WEIO might not be effective. As shown in Fig. 9, some of the CMIP6_Grp models successfully simulate the strong relationship between ISMR-WEIO convection with intense inverse correlation between ISMR and OLR over WEIO, but the correlations are much weaker in their atmosphere-only counterparts (Fig. 9 inset). The correlations reduce and inter-model spread increases if all the CMIP6/ AMIP6 simulations are considered (Fig. S4). AMIP6 models show lower skill in capturing this relationship compared to CMIP6 models. This remains the same with CMIP5 and AMIP5 models as well (Fig. S4). On the contrary, over the eastern part of EEIO, the impact of ENSO is significant and much stronger than its impact over WEIO which is captured well by the majority of the models (Table 7). This strong inverse relationship between ENSO and EEIO rainfall (Fig. 8b) is not out of phase or at odds with the negative or rather weak relationship between ISMR-EEIO rainfall (Fig. 8a). Figure 10a shows the correlation coefficient between normalized AIR anomaly with ENSO index from 1 year before to 1 year after the monsoon season for observation and CMIP6_Grp models. The AIR-ENSO index inverse correlations are the largest during ISM season characterised by the lag-0 correlation between anomalies of AIR and ENSO index. The evolution of the relationship suggests the strong link of ENSO with ISMR during the peak monsoon months of JJAS and the possible role of ISMR variation on the subsequent tropical circulation through associated fluctuations in location and intensity of diabatic heating fields. Many models capture the evolution of the inverse link between ISMR and ENSO, though with a tendency to build-up the relationship before the monsoon season. It can be seen that several models tend to overestimate this negative correlation and the coupled models are more successful in capturing this strong inverse correlation (Fig. 10a inset). Many CMIP6 models have problems in simulating the timing and magnitude of this teleconnection which are at odds with the observation (Fig. S5a). This is consistent with the findings of Annamalai et al. (2007) who have examined the simulations of ENSO-ISMR relationship by IPCC AR4 models. The tendency for coupled models to build-up the relationship much before the monsoon season could be associated with the length of the ENSO events in the model. For example, in a comparison study of ENSO-monsoon relationship in CMIP3 and CMIP5 simulations, Jourdain et al. (2013) suggested that the Indian summer monsoon-ENSO relationship is affected by overly persistent ENSO events in many CMIP models. McKenna et al. (2020) reported the existence of large spread of ENSO periodicity and broad peaks in CMIP6 simulations though these aspects are slightly improved from CMIP5 simulations.
Correlation for JJAS average indices in observation is −0.56 and it ranges from −0.34 to −0.75 with the median around −0.5 in 15 CMIP6_Grp models (Fig. 10a inset). Atmosphere-only models are not as successful in getting this relationship with the correlations of 10 AGCMs ranging from −0.04 to −0.46 with the median around −0.2 . The improved skill of coupled models in capturing the inverse relationship between ISMR and ENSO compared to the atmosphere-only models and the improvement in CMIP6 models compared to CMIP5 models are reflected in the box and whisker plots showing the correlation between ISMR and ENSO index in Boreal summer season of JJAS for CMIP6, AMIP6, CMIP5 and AMIP5 models (Fig. S6) and Table 6 Diagnostic statistics for 15 CMIP6 models selected as a group of models with greater skill in simulating the spatial pattern of climatological JJAS mean rainfall over India; mean rainfall averaged over India (column 3), temporal standard deviation (column 4), Coefficient of Variation (CV) (column 5), PCC IMD / PCC GPCP (column 6), SD (column 7) and BIAS (column 8) in the spread and strength of the correlations (Text S1). The observed correlations of AIR during JJAS with EMI are not that strong suggesting lack of relationship between monsoon and El Niño Modoki (Fig. 10b). However, CMIP6_Grp models show much stronger relationship between ISMR and El Niño Modoki as well. Figure S5b shows that the relationship of ISMR with ENSO/El Niño Modoki is much stronger in CMIP6 models than the AMIP6 models.

ISMR-Equatorial Indian Ocean relationship
Climatologically, EEIO is more active during monsoon season with organized convection/rainfall which is much stronger than that over WEIO (Fig. 1a). The SST is higher and consistently the mean sea level pressure (not shown) is lower over EEIO than WEIO, which results in westerly winds over central equatorial IO (CEIO). But, some years have these reversed with negative SST anomalies over EEIO and positive over WEIO. As a result, convection gets suppressed over EEIO and enhanced over WEIO. This mode of variability is termed as the Indian Ocean Dipole/Zonal Mode (IOD/IOZM) in the tropical Indian Ocean (Saji et al. 1999;Webster et al. 1999). Saji et al. (1999) confirmed the existence of this dipole like structure coupled with surface wind anomalies of CEIO. The condition with low SSTs over EEIO and high over WEIO, is the positive phase of IOD. On the other hand, enhancement of the normal state of the tropical Indian Ocean with below normal SSTs over WEIO and above normal SSTs over EEIO is the negative phase of IOD. Due to the background seasonal changes in the equatorial IO, the IOD is phase-locked with developing during June/ July/August, peaking during September to November and thereafter decaying in December to February (e.g., Fan et al. 2017). Climate models are known to simulate overly intense IOD variability (e.g., Wang et al. 2017), possibly due to the biases in SST, thermocline and wind feedbacks (Cai and Cowan 2013). The strength and seasonality of IOD simulated by CMIP5 and CMIP6 models based on the monthly standard deviation of DMI are depicted in Fig. 11a and b. It can be clearly seen that many models of CMIP5 and CMIP6 show deficiency in representing the strength, onset, peak and decay of IOD. Majority of the models simulate overly strong IOD variability. The changes in the basic characteristics of the IOD from CMIP5 to CMIP6 are very modest. There is an increased tendency for a shift in the seasonality of IOD towards an earlier peak in September in CMIP6 models compared to CMIP5 models.
ENSO is observed to be connected to some of the IOD events and have triggered them or modulated their evolution (e.g., Fan et al. 2017). IOD is thus due to a combination of internally forced and ENSO driven phenomenon and its onset, peak and decay could vary to some extent due to ENSO. Correspondingly, ENSO and DMI are observed to be well correlated (0.45) and CMIP6 models show a slight improvement in capturing this relationship compared to CMIP5 models (Fig. 11c). However, the relationship between ISMR and DMI is rather weak as reported in earlier studies , with a correlation coefficient not significantly different from zero in which only about 3% of the variance of ISMR explained by the variance of DMI (Fig. 11c). The models tend to simulate higher correlation between ISMR and IOD. Gadgil et al. (2004) linked the extremes of ISMR to the variation of convection over EIO termed as the Equatorial Indian Ocean Oscillation (EQUINOO). Its positive (negative) phase is defined as the state of the equatorial IO with enhancement (suppression) of convection over WEIO and suppression (enhancement) over EEIO along with associated changes in the zonal wind anomaly over CEIO. Thus, during the positive phase, the southwesterly cross-equatorial monsoon flow as well as the low-level jet stream carries the excess clouds off to the southern as well as central India complementing the ISM. In contrast, during the negative phase, convection is enhanced over EEIO and westerly would dominate over EIO carrying the excess cloud further to the east and gets drifted off to the Pacific region pushing India to drought condition. Figure 12 shows the schematic diagrams of the positive and negative phases of EQUINOO with the corresponding anomalies of SST, convection/cloudiness and surface wind directions.
EQUINOO is often considered as the atmospheric counterpart of IOD. In ENSO, there is a tight coupling between the Southern Oscillation (SO) in the atmosphere and the fluctuations between El Niño and La Niña in the ocean with the correlation coefficient between their indices being 0.86. However, the correlation between EQUINOO index and the index of IOD is relatively low ( ∼ 0.53) with both having opposite phases for June-September in almost onethird of the years. Also, it is important to note that for the summer monsoon season of JJAS, the EQUINOO index is very poorly correlated with the ENSO index whereas the IOD index is significantly correlated with ENSO index (e.g., Table 2 of Sajani et al. 2015).  . 6 a Box and whisker plots showing the PCCs of mean summer monsoon rainfall over India against the corresponding IMD observation for CMIP6, AMIP6, CMIP5, AMIP5, CMIP3, AMIP3 and AMIP1 models. Whiskers show the maximum and minimum values among the PCCs of the member models and the box shows the interquartile range or the first to the third quartile (second quartile or the median is shown as the solid line inside the box) of the PCCs (top panel). The PCC for the MME of models from each MIP is written in the right side of its box. b Same as a, but for all-India total summer monsoon rainfall ( ) in mm (middle panel). The range of observed estimated for the periods corresponding to different MIPs are shown in the background as grey shaded box. c Same as b, but for the standard deviation of all-India summer monsoon rainfall in mm (bottom panel) The positive phase of EQUINOO with enhanced rainfall over WEIO is favourable for ISMR. The opposite happens in the negative phase. Correspondingly, surface zonal wind anomaly over the central equatorial IO is easterly (westerly) in the positive (negative) phase (Ashok et al. 2004). So CEIO wind-based index is used as the EQUINOO index (as explained in section 2). The correlation between the ISMR and EQUINOO index is found to be significant with the EQUINOO index explaining 19% of ISMR variance by Sajani et al. (2015). They have shown that a composite index derived through the bivariate analysis of indices of ENSO and EQUINOO, explains 54% ISMR variance. The composite index is found to be negative for all the deficit ISMR years and positive for most of the excess rainfall years. In strong El Niño years like1997, its effect on ISMR was offset by the positive EQUINOO during the season and the ISMR was near normal. But in years like 1987 with unfavourable EQUINOO (negative) and ENSO (El Niño), ISMR was deficit by 27%. It is therefore of utmost importance for models to capture not only ISMR's relation with ENSO but with EQUINOO as well to improve the simulation of IAV of ISMR. Thus, modelling the combined effect of ENSO and EQUINOO on ISMR is a topic with high relevance. As seen from Fig. 10, the models capture the opposite relationship between ISMR and ENSO though with an early build-up and overestimation. Still, it is important because the major source of ISMR variability is ENSO forcing. Next, we assess the models' skill in simulating two characteristics of EQUI-NOO; i) the dipole-like structure of EQUINOO in the EIO and ii) the positive link between ISMR and WEIO rainfall/ EQUINOO.
Correlations of WEIO rainfall and EQUINOO Index with tropical rainfall at every grid point show their teleconnections over Indian and Pacific Oceans (Fig. 13). There is a remarkable similarity between the two correlation maps for the two indices of EQUINOO, i.e., the WEIO rainfall and the EQUINOO index. Marked characteristics of EQUINOO are (i) the opposite relationship between rainfall over WEIO and EEIO or the dipole structure of EQUINOO and (ii) the positive correlation between WEIO rainfall and ISMR or EQUINOO Index and ISMR. Another feature is the positive correlation between WEIO rainfall/EQUINOO Index and  4). Correlation coefficients significant at 95% are stippled rainfall over NINO3.4 region. For the dipole structure, the observed correlation is −0.54 and some of the CMIP6_Grp models capture this inverse relationship between WEIO and EEIO rainfall realistically (Fig. 14a). But the correlations are much weaker in their atmosphere-only counterparts (Fig. 14a inset). The correlations reduce and inter-model spread increases if all the CMIP6/AMIP6 simulations are considered (Fig. S7). AMIP6 models show lower skill in capturing this relationship compared to CMIP6 models. This remains the same with CMIP5 and AMIP5 models as well (Fig. S7).
The observed correlation of 0.43 between ISMR and the EQUINOO index for the period 1958-2010 was found to be significant and explaining 19% of ISMR variance by Sajani et al. (2015) based on NCEP reanalysis winds. Using ERA-5 data the correlation is found to be 0.4 which explains 16% of ISMR variance (Fig. 14b). However, most of the models show serious deficiency in capturing the second feature, i.e., EQUINOO's relationship with ISMR. None of the models is capable of simulating the significant positive correlation between EQUINOO index and ISMR and many tend to simulate an inverse relationship. This suggests that several models reproduce the dipole-like structure in the equatorial Indian Ocean, but with opposite phase relationship with ISMR. This might be one of the major reasons why models fail to correctly capture the observed interannual variation of ISMR.
Regarding the relation between WEIO rainfall and ISMR, six of the coupled models of CMIP6_Grp simulate the positive correlation. We compare the coupled and Table 7 Correlation coefficient between time series of ENSO index and EEIO (east pole of EQUINOO) rainfall for June to September and JJAS, in 15 models of CMIP6 group and observation (italics). Models with the best simulation of the distribution of the correlation coefficient between ENSO index and grid rainfall over the tropical Asia-Pacific region are highlighted in bold Fig. 9 Scatter plot of standardised anomalies of ISMR versus WEIO OLR from observation (IMD rainfall and NOAA CDR denoted by the last two digits of the years) and 5 representative models of CMIP6_ Grp. The corresponding correlations are given in top-right corner. Inset: Box and whisker plots showing the correlation coefficient between ISMR and WEIO OLR for 15 CMIP6_Grp models and 10 of their corresponding atmosphere-only models from AMIP6. Whiskers show the maximum and minimum values among the correlations of the member models and the box shows the interquartile range or the first to the third quartile (second quartile or the median is shown as the solid line inside the box). The horizontal black line denotes the observed correlation atmosphere-only versions of these six models which show some skill in capturing correctly the sign of the correlation between WEIO rainfall and rainfall over India (Fig. 15). We can see that many of the coupled versions are more adept in capturing not only ISMR-EQUINOO relationship, but EQUINOO's relationship with EEIO rainfall and rainfall over NINO3.4 region, as well. The AMIP runs are not as skilful in simulating these features. In this context, it is important to note that CMIP6 models are faring better in simulating the mean monsoon rainfall over India, its interannual variability and monsoon teleconnections with ENSO and EQUINOO, than the AMIP6 models (Figs. 4 and 6). But, the lack of skill in capturing the link between ISMR and convection over the EIO consistently in the majority of models, could be a crucial detrimental factor limiting the ability of the coupled models in simulating IAV of ISM.

June AIR link with north-west tropical Pacific rainfall
In recent decades, the extremes in ISMR are often associated with large interannual variation of AIR in June which has an influence on the seasonal rainfall. However, it is not yet understood if June rainfall plays a decisive role behind the rainfall of the rest of the season. Studies are suggesting that rainfall in June has direct links with convection and rainfall over the northwest tropical Pacific (Gadgil and Francis 2016;Sajani et al. 2019). In June, the strongest inverse correlation over the Asia-Pacific is seen between AIR and NWTPR with most of the AIR excesses (deficits) associated with suppression (enhancement) of NWTPR ). This teleconnection is reasonably well simulated by some of the CMIP6_Grp (Fig. 16) models with a strong inverse relationship between the simulated June AIR and NWTPR. Sajani et al. (2019) proposed that NWTPR variation is associated with a shift in the latitudinal location of subtropical westerly jet over the region from West of Tibetan Plateau to NWTP and the phase of Rossby wave steered in it. This feature is also captured well by these models. The skill of atmosphere-only models are comparatively lower in simulating the inverse correlation between AIR and NWTPR in June (Fig. 16 inset). The correlations weakens and inter-model spread increases if all the CMIP6/AMIP6 simulations are considered (Fig. S8). AMIP6 models show lower skill in capturing this relationship compared to CMIP6 models. This aspect remains the same with CMIP5 and AMIP5 models as well. There is a slight improvement in the simulated correlations from CMIP5 to CMIP6.

Concluding remarks
Our study is a performance assessment of global climate models of 6 th generation of Coupled Model Intercomparison Project (CMIP6) in simulating the climatological summer monsoon rainfall over India, the interannual variability (IAV) of all-India summer monsoon rainfall (ISMR) and monsoon teleconnections with the interannual modes of (a) (b) Fig. 10 a Correlation coefficients between Indian summer monsoon rainfall (ISMR) anomaly and NINO3.4 SST anomaly (ENSO index) at various monthly lags for observation (black thick line) and 15 models (blue thin lines) of CMIP6_Grp (top panel). JJAS season is highlighted in green. Top panel inset: Box and whisker plots showing the correlation coefficient between ISMR and ENSO index during JJAS for 15 CMIP6_Grp models and 10 of their corresponding atmosphere-only models from AMIP6. Whiskers show the maximum and minimum values among the correlations of the member models and the box shows the interquartile range or the first to the third quartile (second quartile or the median is shown as the solid line inside the box). The horizontal black line denotes the observed correlation of −0.56 . b Same as a, but for correlation coefficients between ISMR anomaly and El Niño Modoki Index (EMI) at various monthly lags (bottom panel) with the inset panel showing the box and whisker plots of the correlation coefficient between ISMR and ENSO index during JJAS for 15 CMIP6_Grp models and 10 of their corresponding atmosphere-only models from AMIP6  (4) a brief look into how the inverse relationship between AIR and NWTPR in June, is simulated by CMIP6 models. The MME of 61 CMIP6 models is found to have the best skill in simulating climatological summer monsoon rainfall pattern over India compared to the MMEs of the atmosphere-only models of AMIP6 and the models of the previous generations of Atmospheric and Coupled Model Intercomparison Projects (AMIPs and CMIPs). Systematic improvement and reduction in bias are seen from lower to higher AMIPs/CMIPs, more evidently in coupled models. Still, there exists slight dry bias over a narrow region of the . The anomalies of SST (red: warmer and blue: colder), direction of prevalent surface winds (yellow arrows) and cloudiness (silver shade) as proxy for organized convection/rainfall associated with the two phases, are also shown Fig. 13 Distribution of correlation coefficient of a WEIO rainfall versus rainfall at every grid over the tropical Asia-Pacific sector (top panel) and b EQUINOO index versus rainfall at every grid over the tropical Asia-Pacific sector (bottom panel). Black dashed boxes show regions of western and eastern equatorial Indian Ocean (WEIO and EEIO respectively) and equatorial central Pacific (NINO3.4). Correlations significant at 95% are stippled monsoon zone of central India along with wide-spread wet and cold bias over the Indian and Pacific Oceans which is the lowest in CMIP6 models.
A subgroup of 15 CMIP6 models (CMIP_Grp) selected through objective criteria based on their ability to simulate climatological summer monsoon rainfall over India and the IAV of ISMR, perform better in simulating the IAV of ISMR in comparison to the atmosphere-only models. Majority of these coupled models capture the strong inverse correlation between ISMR and El Niño-Southern Oscillation (ENSO) index and the evolution of the inverse link between ISMR and ENSO, though with a tendency to build-up the relationship before the monsoon season with some models slightly overestimating the simultaneous negative correlation in the Boreal summer (JJAS) season. Atmosphere-only models are least successful in simulating this relationship which could be a factor affecting their skill in simulating the IAV of ISMR. The CMIP6 models are also found to simulate well the strong inverse relation between all-India rainfall and north-west tropical Pacific rainfall (NWTPR) in June.
Although the coupled models are successful in simulating the dominant monsoon teleconnections associated with the inverse relationship between ISMR and ENSO in the Boreal summer season, they lack skill in simulating another important teleconnection existing between ISMR and EQUINOO. Many models reproduce the dipole structure associated with the EQUINOO with the out of phase relationship between WEIO and EEIO rainfall. However, the majority of the models fail to capture the link between ISMR and EQUINOO. Or in other words, several models reproduce EQUINOO mode as such, i.e., the dipole-like structure in the equatorial Indian Ocean, but with opposite phase relationship with ISMR. Regarding the relationship of ISMR with rainfall over WEIO which constitutes one pole of EQUINOO, six coupled models of CMIP6_Grp simulate the positive correlation between WEIO rainfall and ISMR. We compared the coupled and atmosphere-only versions of these six models in simulating ISMR-EQUINOO relationship. The coupled versions appear to be more adept in capturing not only WEIO rainfall-ISMR relationship but its relationship with EEIO rainfall and rainfall over NINO3.4 region as well. The AMIP runs are not as skilful in simulating these features.
The lack of skill in capturing the link between ISMR and convection over the equatorial Indian Ocean consistently in the majority of CMIP6 models could be a crucial detrimental factor limiting the ability of the coupled models in simulating IAV of ISMR. Although the CMIP6 models fare better in simulating the climatological monsoon rainfall and its pattern over India, interannual variability of ISMR, the proportion of coefficient of variation to all-India summer monsoon rainfall, the simulation of EQUINOO and monsoon teleconnections with ENSO and NWTPR, failure to simulate the ISMR-EQUINOO link could be one of the major reasons 14 Scatter plot of standardized anomalies of a EEIO rainfall versus WEIO rainfall from observation (GPCP rainfall denoted by the last two digits of the years) and 4 representative models of CMIP6_ Grp (top panel). Inset: Box and whisker plots showing the correlation coefficient between EEIO rainfall and WEIO rainfall during JJAS for 15 CMIP6_Grp models and 10 of their corresponding atmosphereonly models from AMIP6. Whiskers show the maximum and minimum values among the correlations of the member models and the box shows the interquartile range or the first to the third quartile (second quartile or the median is shown as the solid line inside the box). The horizontal black line denotes the observed correlation of −0.54 . b Scatter plot of standardised anomalies of ISMR versus EQUINOO index from observation (IMD rainfall and ERA-5 data). The corresponding correlations are given in top-right corner 1 3 why coupled models fail to capture the observed interannual variation of ISMR, correctly.
The impact of ENSO elicits a response over the equatorial Indian Ocean with strong positive (negative) correlation with rainfall and SST over the WEIO (EEIO) and suppression of monsoon rainfall over the Indian subcontinent. This results in an opposite relationship between WEIO rainfall and ISMR in ENSO years. Strong or stronger than the observed Correlations significant at 95% are stippled relationship between ISMR and ENSO in models can also end them up having the opposite relationship between ISMR and WEIO rainfall if they fail to simulate the basic positive relationship between them as observed. Thus, our analysis implies the need for capturing ISMR-EQUINOO relationship, by models to improve the simulation of the IAV of ISMR which is crucial for reliable monsoon prediction.

Fig. 16
Scatter plot of standardized anomalies of June all-India rainfall (AIR) versus June NWTPR from observation (denoted by the last two digits of the years) and 4 representative models of CMIP6_Grp. The corresponding correlations are given in top-right corner. Inset: Box and whisker plots showing the correlation coefficient between ISMR and NWTPR during JJAS for 15 CMIP6_Grp models and 10 of their corresponding atmosphere-only models from AMIP6. Whiskers show the maximum and minimum values among the correlations of the member models and the box shows the interquartile range or the first to the third quartile (second quartile or the median is shown as the solid line inside the box). The horizontal black line denotes the observed correlation