Indian Ocean dynamic sea level, its variability and projections in CMIP6 models

The regional sea level variability and its projection amidst the global sea level rise is one of the major concerns for coastal communities. The dynamic sea level plays a major role in the observed spatial deviations in regional sea level rise from the global mean. The present study evaluates 27 climate model simulations from the sixth phase of the Coupled Model Intercomparison Project (CMIP6) for their representation of the historical mean states, variability and future projections for the Indian Ocean. Most models reproduce the observed mean state of the dynamic sea level realistically; however, consistent positive bias is evident across the latitudinal range of the Indian Ocean. The strongest sea level bias is seen along the Antarctic Circumpolar Current (ACC) regime owing to the stronger than observed south Indian Ocean westerlies and its equatorward bias. This equatorward shift of the wind field also results in a stronger positive windstress curl across the southeasterly trade wind regime in the southern tropical basin and an easterly wind bias along the equatorial waveguide. Owing to the anomalous easterly equatorial winds, the thermocline in the eastern tropical basin is shallower in the models than observed, resulting in enhanced variability there. Such spurious variability in the eastern part of the basin causes models to become biased towards the dipole zonal mode or Indian Ocean dipole patterns in the tropics. In the north Indian Ocean, the summer monsoon winds are weak in the model leading to weaker coastal upwelling and positive sea level bias along the western Arabian Sea. Further, it is noted that the high-resolution models compare better in simulating the sea level variability, particularly in the eddy-dominated regions like the ACC regime in interannual timescale. However, these improved variabilities do not necessarily produce a better mean state likely due to the spurious enhanced mixing driven by parametrizations set in these high-resolution models. Finally, the overall pattern of the projected dynamic sea level rise is similar for the mid (SSP2-4.5) and high-end (SSP5-8.5) scenarios, except that the magnitude is higher under the high emission situation. Notably, the projected dynamic sea level change is milder when only the best-performing models are used compared to the complete ensemble.


Introduction
The ongoing global warming owing to the increase in the concentration of anthropogenic greenhouse gasses is causing a devastating irreversible impact on the global climate systems (IPCC AR6 2021). About 93% of the extra heat stored in the climate system is absorbed by the ocean, leading to a rapid sea level rise in the recent decade and therefore, posing a great threat to the densely populated low-lying coastal areas across the globe (Oppenheimer et al. 2019). As per the latest report of the IPCC (Intergovernmental Panel on Climate Change; IPCC AR6 2021), the global mean sea level has been increasing at a steady rate of 1.3 mm/year during 1901-1971, 1.9 mm/year during the period 1971-2006 to a much rapid rise of 3.7 mm/year in the recent decades (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Similar observations were echoed by others (Bindoff et al. 2007;Han et al. 2010;Church and White 2011;Church et al. 2013;Unnikrishnan et al. 2015). This rising sea level is expected to increase further at a rapid pace unless a deep reduction in CO 2 emission becomes a reality across the globe (Perrette et al. 2013;Oppenheimer et al. 2019). Moreover, this higher sea level and a warmer climate are expected to exacerbate extreme weather events such as 1 3 storm surges, cyclones, high waves and greater erosion. This will put human lives along the coastal communities at risk of oceanogenic calamities around the globe, particularly in the heavily populated south Asian countries where these effects will be devastating.
More importantly, the regional sea level change differs significantly from the global mean driven by ocean dynamics, differential heating and gravitational and solid earth processes (Lowe and Gregory 2006;Christensen and Christensen 2007;Milne et al. 2009;Yin et al. 2010;Tamisiea and Mitrovica 2011;Yin 2012;Meyssignac and Cazenave 2012;McGregor et al. 2012;Landerer et al. 2007Landerer et al. , 2014. Thus, the aim of this paper is to identify the fidelity of the CMIP6 models in simulating the mean and variability of the sea level in the Indian Ocean and relate the identified bias with possible anomalous physical processes. Recently, Harrison et al. (2021) argued that while stereodynamic processes are the leading contributors to the future sea-level change in the Indian ocean, mass changes due to the terrestrial water storage also play a key role in the spatial variations of sea level. Moreover, this non-uniform sea level change is also the manifestation of both internal variability and anthropogenic forcing (Han et al. 2010). While the contribution of natural internal variability associated with coupled climate modes such as El Niño Southern Oscillation (ENSO), Indian Ocean dipole (IOD), Pacific Decadal Oscillation (PDO), Southern Annular Mode (SAM), etc. (Phillips et al. 2021 and the references therein) is one of the dominating forcings in the interannual/decadal timescale, the contribution of the anthropogenic component is also expected to grow over time owing to future climate change (Stammer et al. 2013). Hence, at a local level, the global mean assessments are of little use to policymakers.
Phase six of the Coupled Model Inter-comparison Project (CMIP6), like its predecessors, brings standardized model outputs of the historical simulations as well as the future projections from various participating institutes in the project (Eyring et al. 2016), allowing a detailed regional specific investigation of the sea level change and associated processes. The model simulated dynamic sea level (DSL) represents the change in sea level due to the thermal and halosteric changes along with dynamic processes associated with wind forcings (Stammer et al. 2013). Since the Indian Ocean has been subjected to rapid warming in the last few decades (Alory et al. 2007;Alory and Meyers 2009;Roxy et al. 2014;Li et al. 2017;Chatterjee et al. 2020, preprint) it is important to understand the signature of the dynamic sea-level rise and its projections over the region. Earlier global mean state simulation of sea level from CMIP3 to CMIP5 has shown a marked improvement (Landerer et al. 2014). Note, however, that in the CMIP5, the future projections were driven by different climate scenarios based on the representative concentration pathways (RCPs) associated with varied radiative forcing ranging from 2.6 to 8.5 Wm −2 (Van Vuuren et al. 2011). In contrast, CMIP6 models adopted a new set of future climate scenarios based on the Shared Socioeconomic Pathway (SSP) to additionally consider land use, urbanization, economy, population, etc. (O'Neill et al. 2014;Riahi et al. 2017). Therefore, the narrative of the projected changes under the CMIP6 scenarios are similar to the CMIP5, but not the same. We expect that the projected regional change in sea level may differ under these new SSP considerations compared to the previous estimates. Lyu et al. (2020a, b) argued that CMIP6 lacks a similar improvement over CMIP5 despite its better representation of the zonal wind stress, especially the position of the southern hemisphere westerlies which improve DSL estimates in the Southern Ocean. CMIP6 models also report increased DSL projections compared to the CMIP5 models likely linked to the inclusion of models with higher climate sensitivity (Andrews et al. 2012;Lyu et al. 2020a, b;Brunner et al. 2020;Nijsse et al. 2020;Tokarska et al. 2020). In fact, Hermans et al. (2021) noted that while the increased effective climate sensitivity of CMIP6 models relative to CMIP5 models translates into substantially higher projection of twenty-first century global surface temperature (GSAT), the difference in projected global mean sea level (GMSL) is relatively modest when methods from fifth assessment report (AR5) is used.
Notably, proper mean state representation is an important requisite for models simulating future projections (Richter et al. 2017). Along with mean state simulations, an accurate representation of internal variability is also imperative for better projections. However, the ensemble spread of the climate models increases the uncertainty of the projected change. However, it is to be noted that the noise does not vary much with the choice of emission scenarios; hence, the uncertainty associated with internal variability is independent of the choice of the emission scenarios (Ferrero et al. 2021). To reduce this uncertainty, many previous studies have explored the performance-backed weightage or exclusion of models (Greene et al. 2006;Yin et al. 2010;Watterson and Whetton 2011). Most of these assessment studies for different variables and regions employed subjective ranking of models based on their ability to simulate largescale climate patterns and skill score metrics (McSweeny et al. 2015;Halder et al. 2021;Krishnan and Bhaskaran 2019). McSweeney et al. (2015) provided a rationale that the inability of the model to simulate key large-scale processes may reflect the model's deficiency in simulating climate signals induced by global warming. Nevertheless, finding a selection rationale that best suits the proposed purpose is challenging. Note, however, that while such exclusion of non-performing models for the ensemble projection reduces the uncertainty and narrows down the projected spread, it also increases the risk of mal-adaptation (McSweeney et al. 2015). However, as the cost of implementing mitigation policies to safeguard a coast from rising sea levels is expected to increase exponentially with the projected future sea level rise, lowering the uncertainty is likely to provide a sustainable solution. Further, as the oceanographic operational centres across the globe are more interested in downscaling the global sea level rise using dynamical and statistical approaches (Liu et al. 2016;Hermans et al. 2020;Jackson et al. 2022), the need to select the best-performing models for the region of interest is growing for such analysis.
In this study, we analyze historical sea level simulations from 27 CMIP6 models to identify their fidelity in simulating the observed mean state and variability of the Indian Ocean. Further, we select a subset of models based on their performance for various climatic features of the Indian Ocean. We expect that the exclusion of such leastperforming models from the ensemble may likely provide policy-relevant projections over the Indian ocean domain. Further, considering the different driving dynamical mechanisms responsible for mean sea level and its variability, we discuss the model biases and their underlying mechanisms separately for the north Indian Ocean (NIO; 0°-30° N, 20°-120° E), south tropical Indian Ocean (STIO; 0°-30° S, 20°-120° E) and Southern Indian Ocean (SIO; 30°-60° S, 20°-120° E). Note here that the regions west of 120° E meridian, but outside the Indian Ocean such as parts of South China Sea, Java Sea, etc. are excluded from the study. The rest of the paper is organized as follows: Sect. 2 describes the observation and models we use in this study, assessments of each model in simulating mean, variability and process-specific Indian Ocean features are discussed in Sect. 3, the projected changes are discussed in Sect. 4, and finally, Sect. 5 concludes our results.

Data and methodology
This study uses monthly outputs from 27 CMIP6 models (Eyring et al. 2016; Table 1). To assess the fidelity of each model in reproducing the Indian Ocean mean state and its variability, the 'historical' simulation (available from 1850 to 2014) of sea surface height above geoid (zos), the global mean thermosteric sea level change (zostoga), zonal and meridional windstress at the surface (tauu, tauv), and derived curl of the wind stress are analyzed. For future projections, we examine two Shared Socioeconomic Pathways (SSP) scenarios from Scenario Model Intercomparison Project (ScenarioMIP;O'Neill et al. 2014) for a possible midemission (SSP2-4.5) and high-emission (SSP5-8.5) radiative forcing at the end of 2100. To give equal weightage to all the models, only the first realization 'r1i1p1f1' is used for both historical and future scenarios. The dynamic sea level (DSL) is defined as the local height of the sea level above the geoid.
The DSL is supposed to have a global zero mean. Some models provide a non-zero time-dependent global mean for the variable zos, which we remove for consistency in analysis (Gregory et al. 2019;Griffies et al. 2016). Outputs from certain models need barometric correction to convert them into effective sea level to make inter-model comparison as well as model to observation comparison feasible (Griffies et al. 2016). We do this correction by converting the mass per area occupied by the sea ice into the height it would occupy as water and subsequently adding it to the modelled sea surface height for the specific models. Furthermore, these model simulations also suffer from a model drift caused by a relatively short spin-up of the subsurface water column (Gupta et al. 2013). The historical simulations are branched from a specific point in time from the pre-industrial control run. We identify this branching point in time and find out the The altimeter-derived multi-mission gridded product for 1993-2014 is used to compare the simulated model DSL. As the altimeter-derived sea level is found to have a correlation greater than 0.95 (not shown) with the mean dynamic topography obtained from Maximenko et al. (2009), the use of altimeter sea level anomaly does not affect our analysis. Since the altimetry record contains the global mean thermosteric component (Stammer et al. 2013), we remove the global mean sea level at each time step from the observation to make the observation-model comparison feasible. Noticeably though, as the altimeter sea level observation is available only for a shorter duration, the observed variability may be biased by the unresolved decadal climate variability.
The interannual variability of DSL is calculated after removing the monthly climatology signal from the total detrended DSL. The climatological DSL is used to study the annual variability. Moreover, since the models have no real state initialization, the phase of the model simulations is expected to differ from that of the altimeter. To address this issue, we use a 20-year sliding window on the entire historical simulation record and compare each of these windows with the observational data of the same length to see whether the choice of the observational period influences the results.
Since the DSL is driven by dynamical forcing, we also analyze the wind stress of each of the CMIP6 models. Wind stress curl is derived from the historical and future scenario simulations of the models as it is one of the primary drivers of the mean dynamic sea level and its variability. To compare the wind field and its anomaly European Centre for Medium Range Weather Forecasts (ECMWF) Reanalysis (ERA5) from 1993 to 2014 is also used (Hersbach et al. 2020).
In order to assess the performance of each model, basic statistical metrics such as pattern correlation, root mean square error (RMSE), Empirical Orthogonal Function (EOF), etc., are used. The pattern correlation is used to assess the fidelity of a model in simulating the spatial pattern of mean sea level and variability (standard deviation) over a specified 20-year sliding window compared to the observations. Similarly, RMSE is also calculated between spatial maps of model and observation for the field of mean or variability of DSL. In contrast, for EOF calculations, detrended monthly field anomalies (by removing the seasonal cycle) are used.Note here that standard deviation of the MMM is calculated by averaging standard deviation of each model. A similar approach is applied for the MMM EOF calculations as well.
We also use a skill score metric defined by Taylor (2001)  as given below: where ̂ f is the ratio of the spatial standard deviation of the model ( m ) and observation ( o ). R is the observed pattern correlation. R o is the maximum attainable pattern correlation and is set to 1. Spatial standard deviations are calculated using a 2-D array of the spatial map using the following formula: where xy is the spatial standard deviation for a 2-d array (for example, mean DSL over a time window), i,j is the value of the variable (say, DSL) at each grid point and is the average over the domain. n and m are the numbers of grids in the longitude and latitude axis.

Mean dynamic sea level
The mean DSL represents the geostrophic circulation driven by the overlying windstress, windstress curl, and the ocean density (Gregory et al. 2019) of the basin and therefore, plays a key role in the redistribution of heat and mass across the basin. While in the SIO, the DSL is primarily driven by the semi-permanent mid-latitude Southern Ocean westerlies, in the STIO it is driven by the southeasterly trade winds. On the other hand, DSL in the NIO is driven by the seasonally reversing monsoon winds which reverse from southwesterly in the boreal summer to northeasterly in the boreal winter (Schott and McCreary 2001;Schott et al. 2009). Hence, in the NIO, the sign of the DSL also reverses with the season. However, as the winds are much stronger during the boreal summer, the annual mean DSL and winds are biased towards the summer monsoon season of the NIO (Fig. 1). Remote influences from the Pacific Ocean and forcings associated with climate modes like ENSO, IOD and SAM also influence the sea-level in the Indian Ocean through both oceanic pathways and atmospheric teleconnections (Trenary and Han 2012;Phillips et al. 2021;Schott et al. 2009;Sprintall and Révelard 2014;Duan et al. 2021).
All the CMIP6 models could reproduce the mean DSL and wind field reasonably well across the basin. Notably, while the MMM could capture the mean sea level patterns very well, there is a consistent positive bias across the basin for most of the models (Figs. 1c and 2a). The strongest positive sea level bias is evident across the Antarctic Circumpolar Current (ACC) regime in the SIO (Fig. 1c). This region is driven by strong westerlies, causing a marked shift in the sign of the windstress curl across the wind field maxima (Fig. 2b, c). In response to the winds, the sign of DSL across the ACC regime also shifts from negative in the south to positive in the north of its axis. However, note here that the latitude of zero windstress curl does not coincide with the latitude of zero DSL (see Sect. 3.4). Most climate models tend to simulate an equatorward bias in the westerly wind field, causing an equatorward shift in the ACC regime relative to the observation (Fyfe and Saenko 2006;Russell et al. 2006;Lyu et al. 2020a, b). While this equatorward westerly wind bias across the Southern Ocean is reduced considerably in the 6th phase of the CMIP models compared to its predecessors (~ 2°; Lyu et al. 2020a), the latitude of maximum zonal windstress is still off by more than 1° in the MMM in the Indian Ocean sector of the mid-latitude Southern Ocean basin. Furthermore, the westerly winds also show a strong positive bias across the Southern Ocean with a maximum positive bias of ~ 0.1 Nm −2 along its core (Fig. 2b). This strong wind bias also causes a very strong negative (positive) windstress curl bias in the south (north) of the zero windstress curl region (Fig. 2c). A detailed analysis of the Southern Ocean sea level bias and its association with the wind field is discussed in Sect. 3.4. In the STIO, particularly in the western part, the model simulated sea level is slightly higher compared to the observation, mainly linked to the stronger southeasterly trade winds in the model. However, in the NIO, weaker modelsimulated summer monsoon winds (Figs. 1 and S1) produce weak upwelling along the coast of Somalia and Arabia, causing a positive bias in the model-simulated mean sea level in the western Arabian Sea (Fig. 1c). In most models, this also led to a contrasting east-west bias in the sea level across the NIO.
The performance of the individual models differs considerably in the various latitude bands. Figure 3 shows the spread of the sea level bias of the individual CMIP6 models for the NIO, STIO and SIO. Overall, bias is larger in the western part of the basin across all latitudes and generally increases from north to south. While in the NIO, most models show RMSE (with respect to altimeter) less than 0.15 m except for the model INM-CM5-0, MRI-ESM2-0 and GISS-E2-1-G; in the SIO, the RMSE is predominantly more than 0.15 m and reaches close to 0.3 m for the GISS-E2-1-G. The RMSE for the MMM in the NIO, STIO and SIO are 0.09 m, 0.12 m and 0.16 m, respectively, and the corresponding ratio of bias to the observation (ratio of the RMSE to the observation) are 0.19 (0.23), 0.28 (0.30), and − 0.26 (− 0.38), respectively.

Variability
In this section, we analyze the ability of the CMIP6 models to simulate the key variability features in the Indian ocean domain. The sea level over the Indian Ocean shows strong variability (standard deviations) in a large spectrum of timescales. While the variability in the seasonal timescale is most prominent due to the dominance of the seasonally reversing monsoon winds in the north Indian Ocean (Chatterjee et al. 2012), the variability in the interannual and decadal timescale is of key importance for the Indian Ocean climate variability and its impact on the global climate (Schott et al. 2009;Phillips et al. 2021).
Observed variability of the Indian Ocean shows maxima along the front of the ACC in the SIO (Fig. 4a) owing to the interannual meridional movement of the ACC current regime (Fig. 4g) driven by ENSO and SAM climate modes (Kim and Orsi 2014). The thermocline ridge region in the northern part of the STIO also exhibits stronger variability, primarily influenced by the South Equatorial Current (SEC), Bottom panels show RMSE for individual models for each latitude bands. RMSE of the MMM is marked in red overlying wind fields and the associated embedded eddies in seasonal and interannual timescale Deepa et al. 2021). The region off the west coast of Australia in the latitude band of 15°-35° S also shows strong interannual variability (Fig. 4g) driven by western Pacific winds (Trenary and Han 2012) and the eddies generated due to the baroclinic instabilities between the surface eastward flowing South Indian Counter Current (SICC) and subsurface westward flowing SEC (Jia et al. 2011;Menezes et al. 2014;Zhang et al. 2020). In the NIO, the western boundary of the Arabian Sea along the coast of Somalia and Arabia shows strong variability in seasonal and interannual timescale driven by the summer monsoon winds and its climate variability (Beal and Donohue 2013;Chatterjee et al. 2019;Lakshmi et al. 2020;Vinayachandran et al. 2021). Also, the coastal waveguide along the perimeter of the Bay of Bengal exhibits significant variability driven by coastal propagation of planetary waves forced by the equatorial winds (Vialard et al. 2009;Suresh et al. 2013;Chatterjee et al. 2017;Mukherjee et al. 2018;Phillips et al. 2021).
The historical simulations of CMIP6 models are evaluated for the Indian Ocean variability for the common period of data availability from altimeter and CMIP6 models, i.e. during 1993-2014 (Figs. 4b, e, h andS2). Notably, the phase and magnitude of the internal variability in the model may differ with observation due to freely evolving internal variability in the system. Therefore, one can expect the modelsimulated variability to differ considerably during a given  (Richter et al. 2017). In order to account for the contribution of internal variability, pattern correlation and RMSE of the model simulated sea level for a 20-year sliding window are also analyzed for each model, for all frequency bands, to see if choosing a specific period affects the models' skill in simulating the sea level variability.
Most of the CMIP6 models and, therefore, the MMM underestimate the Indian Ocean variability (Fig. 4b). Nevertheless, models perform relatively well for the seasonal scale with an almost accurate representation of the sea level variability in the western Arabian Sea associated with the Great Whirl fronts, its extension along the summer monsoon currents and along the perimeter of the Bay of Bengal (Fig. 4e). However, all the models tend to produce much stronger variability in the eastern equatorial Indian Ocean and along the coast of Java, likely due to the anomalously shallow thermocline (discussed in the next section) leading to stronger air-sea momentum exchange in the upper water column and hence, anomalously stronger excitation of lowerorder baroclinic modes.
In the interannual timescale, models severely underperform in simulating variability along the fronts of the ACC current regime. Similar weaker variability is noticed for the STIO along the band of the SICC regime as well. Notably, the high-resolution models like CMCC-CM2-HR4, CNRM-CM6-1-HR, HADGEM3-GC31-MM and MPI-ESM1-2-HR performed relatively better in reproducing this observed variability, particularly along the ACC current regime. This indicates that the eddy-dominated regions are simulated well by the high-resolution models compared to the coarser CMIP6 models. As this eddy-induced highfrequency variability causes strong low-frequency rectification through modifying the mean state and impacts the response of the internal variability, the natural variability of the Indian Ocean is likely underestimated in most CMIP6 models. For the thermocline ridge region, the variability in MMM is more spread out compared to the observation with a decrease in the magnitude of the variability eastwards. The CMCC models capture the variability in the eastern basin but overestimate in the west. On the other hand, the highresolution CNRM model does capture the variability pattern, but the simulated magnitude of the variability is lower than the observation. The HadGEM3-GC31-MM model fares better compared to others in simulating the observed location and magnitude of the variability in this region (Fig. S2). Figure 5 shows the pattern correlation calculated between the standard deviation in time for each sliding window with the observed temporal standard deviation. Interestingly, the spread of the pattern correlation and RMSE across the sliding windows are within the range of about 0.1 and 0.25 cm, respectively, indicating that the choice of a specific period (window) does not alter the conclusions. Notably, models with higher correlation tend to have smaller RMSE for total and interannual variability, which is not so evident for the annual signal.

Equatorial tilt and climate variability
Next, we compare the model simulated mean DSL for the equatorial belt. The east-west sea level gradient is an important parameter as it is dynamically linked to the depth of the thermocline and therefore, modulates the air-sea interactions over the basin. This mean state particularly becomes important for the initiation and progression of tropical climate modes such as El Niñ o, Indian Ocean Dipole, etc. (Ham and Kug 2015;Cai and Cowan 2013). Earlier, Lyu et al. (2020a) showed that the new generation CMIP6 models could not produce better results compared to the CMIP5 models in simulating the observed equatorial sea level gradient. Here, we try to assess individual models' mean state and relate that to the model performance pertaining to the simulated tropical climate variability. The mean DSL is averaged over 2° N-2° S, and the zonal mean is removed to compare with the observed sea level slope (low in the west and increasing eastward). As the equatorial westerlies primarily drive eastward equatorial upsloping of mean sea level (Wyrtki 1973), zonal windstress averaged over the equatorial band is also analyzed in Fig. 6.
The equatorward bias of the SIO westerlies and STIO trade winds leads to unrealistic anomalous equatorial easterlies in most models (Fig. S1) and thus, in the MMM wind field (Fig. 1), causing underestimation in the strength of the westerlies along the equator. The simulated bias in the equatorial sea level is consistent with the bias in the zonal wind field. The MMM sea level is primarily flat along the equator and completely failed to produce the observed eastward upsloping of the sea level. In fact, 9 models show a negative correlation in producing mean sea level relative to the observation with RMSE of more than 0.5 m. In general, the models with higher correlation show lower RMSE in equatorial sea level. But, in contrast to the variability, the high-resolution models do not exhibit any superiority over the coarser models in simulating the sea level slope. The bias in the zonal windstress corroborates well with the sea level bias (Fig. 6b). However, for winds, unlike sea level, a higher positive correlation does not correspond to a lower RMSE.
Notably, the anomalous equatorial easterlies in CMIP6 models cause an IOD like pattern with cooler SST in the east and warmer in the west (Figure not shown). In most models, these anomalously weaker westerlies cause the thermocline to remain too deep in the western equatorial IO and unrealistically shallow in the eastern side. In the western Indian Ocean, a similar warm SST bias in the CMIP5 models was noted by Fathrio et al. (2017). On the other hand, the shallow thermocline depth in the east leads to anomalously strong air-sea interactions (Cai and Cowan 2013). To understand how this underestimation of equatorial slopes affects the simulation of internal climate modes, we perform EOF on the detrended model simulations after removing the annual cycle and compare it with observations (Fig. 7). The patterns of the first and second modes of EOF are captured reasonably well by the CMIP6 models. While the MMM could capture the observed spatial patterns of the first and second modes of EOF, there is a considerable departure in explaining the spatial magnitude and variability. The first mode of EOF, which represents the dipole mode of the tropical Indian Ocean, explains about 30% of the MMM variability compared to the 16% in the observation. As noted above, the warm SST bias in the west and oversensitive air-sea interactions in the east are the likely reason for the CMIP6 models' overestimation of IOD zonal mode relative to the model simulated total variability in the tropical Indian Ocean. The variability explained by the second EOF is also high in the MMM compared to the observation. However, note here that the spatial amplitude of both modes are considerably weaker than the observation. This weak amplitude is particularly noticeable in the thermocline ridge region, explaining the weaker model-simulated variability noted in Sect. 3.2.

3
As expected, the pattern correlation for the first EOF for all the models is above 0.7 (Fig. 8). Interestingly, the correlation varies significantly for the second EOF across the models, with a maximum correlation of 0.8 for CanESM5-CanOE and close to zero correlation for EC-Earth3 and EC-Earth3-Veg. In the case of explaining variance CNRM-CM6-1-HR and HadGEM3-GC31-MM show weaker variability than observed for both the EOF modes, indicating a possible inadequate climate response in these models. In contrast, NorESM2-MM, GISS-E2-1-G and the models from the EC-Earth3 family show almost three times the observed variability for the first mode suggesting a strong bias towards the IOD mode in these models.

Remote and local forced response
The interannual and decadal variability in the Indian Ocean is majorly forced by local wind regimes and the Pacific teleconnections Han 2012, 2013;Volkov et al. 2020;Phillips et al. 2021;Li et al. 2022). The variability in the Indian Ocean is largely affected by climatic forcings associated with climate modes like ENSO, IOD and SAM through both oceanic and atmospheric pathways (Schott et al. 2009). Sea level signals forced in the western Pacific transmit to the Indian Ocean via the Indonesian throughflow (ITF) and cause the sea level variability in the STIO (Godfrey and Golding 1981; Godfrey 1996; Wijffels and Meyers 2004;Lee et al. 2002) . On the other hand, the local atmospheric forcing in the basin contributes to the variability  (2012) analyzed the relative importance of local versus remote forcing in the STIO using a suite of model experiments. They showed that in the interannual and seasonal time scale the variability over the thermocline ridge is forced by local winds, and the influence of Pacific forcing is confined mainly to the southeastern Indian Ocean (SEIO). Nagura and McPhaden (2021) extended the analysis to the midlatitudes using a simpler 1-D linear, long Rossby wave model and argued that the northern half of the STIO (10°S-18°S; an overlapping region with thermocline ridge) is mainly driven by the local winds in agreement with Trenary and Han (2012) and the southern latitudes (20°-35° S) is dominated by the long Rossby waves propagation radiated from the west coast of Australia. They further added that this variability along the Australian west coast is, in fact, driven by the SSH variability forced in the tropical Pacific.
Hence, next we analyze the fidelity of each model in simulating the local and remotely forced sea level variability in the thermocline ridge of the Indian Ocean (TRIO) and SEIO regions, respectively. Additionally, we also look at the wind stress biases over the Southern Ocean on a global scale to see the relative biases across the ocean basins.

Thermocline ridge of the Indian Ocean (TRIO)
Here we look into the sea level variability over the TRIO characterized by a relatively shallow thermocline compared to the rest of the basin. Note here that, for our analysis, the latitudinal and longitudinal extent of TRIO is defined as 5°-12° S and 50°-80° E following Trenary and Han (2012). It has been well documented that the Ekman pumping associated with prevailing negative windstress curl in the trade wind regime is the primary forcing for the shoaling of the thermocline in TRIO (McCreary et al. 1993;Tozuka et al. 2010;Trenary and Han 2012). Moreover, remotely forced Rossby waves in the eastern part of the basin also contribute to its variability (Han et al. 2014;Mukhopadhyay et al. 2022). The impact of Pacific influence through the Indonesian throughflow is, however, minimal in this TRIO latitudinal band (Han et al. 2007;Yokoi et al. 2008;Nagura and McPhaden 2021;Kersalé et al. 2022). A lead-lag correlation analysis between the observed windstress curl and sea level averaged over the TRIO box indicates about a 4 months lag in sea level from the windstress curl yields a maximum correlation of 0.5 (Fig. 9a). MMM simulates windstress curl, the associated sea level signal and the time lag reasonably well across the sliding windows with a minimum correlation of 0.34 and a maximum of 0.57. Note, however, that the individual models show a large spread across the sliding windows and from the other models (Fig. 9b). The models with higher correlation tend to show less spread compared to the models with weaker median correlation. Moreover, the lag corresponding to the maximum correlation for the respective model is around 4 months or in general longer (Fig. 9c). Among all the models, CMCC, CNRM and EC-Earth group of models perform reasonably well in simulating the correlation and the lag between the TRIO windstress curl and sea level. Interestingly, the variability in TRIO in the interannual timescale is increasingly linked to the westward propagating Rossby waves radiated by the anomalous winds over the eastern tropical Indian Ocean (Trenary and Han 2012;Zhang et al. 2019;Mukhopadhyay et al. 2022). Therefore, the large spread in correlation across the different sliding windows in a few models might be linked to the simulated interannual variability in the winds in the eastern tropical Indian Ocean and the variability in the ITF influx into the Indian Ocean.
Southeastern Indian Ocean (SEIO) The SSH variability in the western equatorial Pacific propagates through the ITF into the Indian Ocean and along the west coast of Australia (Clarke 1991;Feng et al. 2010Feng et al. , 2011Furue et al. 2017;Menezes and Vianna 2019). These signals then propagate westward to influence the interior basin of the south tropical Indian Ocean (Nagura and McPhaden 2021). Further, using sensitivity model experiments, Nagura and McPhaden (2021) demonstrated that the zonal wind filed in the western Pacific associated with ENSO drives the interannual sea level signals in the SEIO region.
In this section, we have evaluated CMIP6 models in simulating this remote forcing by the winds from the western equatorial Pacific (5° S-5° N, 160° E-170° W) in the sea level variability in the SEIO. For consistency, we define the SEIO as the region between 20°-35° S and 105-115° E, following Nagura and McPhaden (2021). Lead-lag correlation analysis of the observed western Pacific zonal winds and SEIO sea level shows a maximum negative correlation of ~ 0.9 when the Pacific zonal winds lead the sea level by 5 months (Fig. 10a). The CMIP6 models simulate the time lead between the Pacific winds and SEIO sea level reasonably well for the observed Fig. 8 a Pattern correlation and b explained variance percentage for the first and second EOF modes of the detrended interannual dynamic sea level from CMIP6 models relative to the altimeter observations. The blue and green dashed lines represent the explained variance percentage for the first and second mode relative to the observed DSL, respectively period, but the simulated correlation remains much weaker and attains a maximum negative correlation of ~ 0.7 for the MMM across the sliding windows. This weaker correlation between the Pacific winds and the SEIO sea level is consistent across most of the models. Only a few models, such as CMCC models, GISS-E2-1-G, MIROC6, MPI-ESM and NorESM2-MM, show a similar or larger correlation than the observed (Fig. 10b). Further, most of the CMIP6 models show a longer lag than the observed value (Fig. 10c), resulting in the MMM attaining the maximum negative correlation for a lag of ~ 7 months, i.e. ~ 2 months delayed than the observed (Fig. 10a). This indicates that the interbasin connection between the Pacific and the Indian Ocean though the ITF channels is not well represented in most of the CMIP6 models and hence, underestimates the ITF influence on the Indian Ocean sea level variability.
Southern Ocean The DSL over the Southern Ocean is a manifestation of the eastward geostrophic flow associated with the ACC and primarily driven by the westerly zonal winds and its meridional gradient (Nowlin and Klinck 1986;Olbers et al. 2004;Bouttes et al. 2012). As noted in Sect. 3.1, the mean DSL of the MMM show the maximum positive bias in the ACC current regime in the 35°-55° S (Fig. 11). In fact, this positive sea level bias is most strongest in the Indian Ocean sector of the Southern Ocean compared to the other basin. The midlatitude westerlies and the Fig. 9 a Lead-lag correlation between DSL and windstress curl averaged over TRIO (5°-12° S and 50°-80° E) region. Shaded curve represents the spread of the correlation for MMM for 20-year windows. b Box plot depicting maximum correlation for each model across the sliding windows and the black dashed line represents the observed maximum correlation. c Box plot showing the lag corresponding to the maximum correlation, the black line is the observed lag corresponding to the maximum correlation Fig. 10 a Lead-lag correlation between DSL in the SEIO and zonal windstress from the western equatorial Pacific (5° S-5° N, 160° E-170° W). Shaded region represents the spread of the correlation for MMM for 20-year sliding windows. b box plot depicting maximum correlation for each model across the sliding windows and the black dashed line represents the observed maximum correlation. c Box plot showing the lag corresponding to the maximum correlation and the black line is the lag of the maximum observed correlation Fig. 11 a Bias in the DSL for the MMM. Dashed black and dark green lines represent the latitude of zero crossing for the observed DSL and wind stress curl, respectively and the grey and the light green line represent the zero crossing for DSL and wind stress curl, respectively for the MMM. b Same as (a) but for zonal wind stress. c Same as (a) but for wind stress curl associated windstress curl also show the strongest bias in this part of the Southern Ocean. In fact, the equatorward shift of the maximum zonal windstress is more pronounced in the Indian ocean compared to the Atlantic and Pacific (Fig. 12). More than 20 models (~ 75%) in the ensemble show an equatorward shift of ~ 2°, resulting in a large patch of positively biased winds centered between 45° S and 50° S (Fig. 11b). Among them, CNRM-CM6-1-HR shows the largest equatorward bias of ~ 5° across the entire basin. This equatorward disposition of the maximum zonal windstress in the Indian Ocean suggests that the DSL biases present in the models are predominantly a local response.
Note that the change in the sign of the DSL (also referred to as the subtropical fronts (STF) or the southern boundary of the subtropical gyre) observed at 45° S and, therefore, does not coincide to the latitude of the observed zero windstress curl (WSC) located at ~ 52° S. In other words, STF lies ~ 7° north of the zero WSC (Figs. 2 and 11a). The MMM show a poleward bias of the STF in the Southern Ocean sectors of the Indian Ocean and the Atlantic Ocean, but, exhibits an equatorward bias for the western Pacific. This observed latitudinal difference between the latitude of STF and the zero WSC was earlier discussed by de Boer et al. (2013). They showed that the bottom pressure torque driven by the topographic gradient determines the flow in the ACC regime. Considering that the bottom pressure torque is a product of the gradient of bottom pressure and topography, the anomalous geostrophic sea level driven by westerly wind bias is the likely source in the model simulated bottom pressure torque and the corresponding bias in the ACC regime. Further, the strong positive bias in the windstress curl north of the zero WSC (Fig. 2c) also favours positive DSL in this region and thus contributes to its positive bias.

Skill score
Most of the models produce a very high skill score of more than 0.98 in simulating the mean DSL for the entire Indian Ocean basin (Fig. 13a). In fact, the minimum score is 0.95 which is associated with the INM-CM5-0 model. The skill scores for the mean windstress curl don't fare as good as the sea level and lie between 0.6 and 0.8, with a few outliers. Even though the models show a good spatial correlation for windstress curl, the larger and varying normalized standard deviations cause the inter-model spread in the skill scores. Note, however, that the skill score is dependent on the choice of the reference dataset as well. For example, the model windstress is more correlated to QuikSCAT scatterometer winds than the ERA-5 and produces a marginally higher skill score (Figure not shown). Nevertheless, while the choice of reference dataset yields marginally different skill scores, the outlying behaviour of most models remains the same. Since ERA-5 is available for a longer record covering our study period (1993( ) compared to QuikSCAT (1999, we prefer ERA-5 for this skill score analysis.
Notably, the skill scores differ significantly across the different latitude belts. For example, the skill score in simulating mean DSL is generally relatively low for the NIO ranging from 0.75 to 0.98 with only three models scoring more than 0.95. This lower skills of the models in NIO, despite relatively lower bias compared to the southern latitudes, is owing to the weaker pattern correlation and spatial standard deviation compared to the observation. In the STIO, most models score higher in simulating sea level, suggesting that most models simulate the subtropical gyre reasonably well in the Indian Ocean. Skill scores remain high in the SIO, with scores of more than 0.85 except for a few outliers. Interestingly, here the model skill in simulating windstress curl is the maximum with a usual score of more than 0.85 despite models producing predominantly anomalously strong westerlies. This discrepancy is due to the fact that the skill score depends only on correlation and its variability and not the magnitude of the variable. Hence, while skill score can provide a picture of the overall performance of each model, it should be considered along with other statistical/ physical characteristics of model simulations.

Projection
Most CMIP6 models conserve volume rather than mass and therefore, any change in the temperature of the water column does not necessarily reflect in the sea level change. In order to account for this effect, we followed the method adopted in Richter et al. (2017). The global mean thermosteric sea level (zostoga) de-drifted and referenced to the period 1993 to 2014 is added to the DSL change. This allows the regional dynamical change to reflects the changes due to winds, regional thermal expansion and the global mean steric height. Note here that the global mean thermosteric sea level change (zostoga) is defined as the change in global ocean volume due to temperature change divided by the ocean surface area (Gregory et al. 2019). MRI-ESM2-0 could not be considered for the analysis due to inconsistency in zostoga (MRI-ESM2-0 shows a discontinuity in the zostoga time-series, Figure S3) and for a few models future projections of DSL are not available at the time of this analysis (Table 1).
From the skill score analysis of the mean sea level and the mean windstress curl, we select the best 15 models by calculating the resultant magnitude of the skill score of both the mean sea level and wind stress curl. The square root of the squared sum of the skill score is calculated and the 15 best models from the ensemble are then selected for further analysis. These models are CMCC-CM2-SR5, HadGEM3-GC31-LL, IPSL-CM6A-LR, CNRM-CM6-1, CNRM-ESM2-1, NorESM2-MM, CMCC-CM2-HR4, BCC-CSM2-MR, CNRM-CM6-1-HR, EC-Earth3, NorESM2-LM, BCC-ESM1, EC-Earth3-Veg, EC-Earth3-Veg-LR and GISS-E2-1-G. After this initial selection of models, we move on to an objective ranking based on the simulated equatorial sea level slope and the total mean sea level. The RMSE of both equatorial sea level slope and basin-wide mean sea level is normalized by dividing them with the respective standard deviations of the selected models for both parameters. These are added together for the respective models to produce an RMSE index which we use to further rank the selected models (see Table 1). We use this index to discard models with RMSE values that fall in the top 25 percentile. As a result, BCC-ESM1, EC-Earth3-Veg, EC-Earth3-Veg-LR and GISS-E2-1-G were removed from the final selected models. Further, since the position of the Southern Hemisphere westerlies is an important climatic parameter, we look at the simulated position of the westerlies in the models to see whether any of the selected models is an outlier in this parameter. It is found that CNRM-CM6-1-HR shows an equatorward shift of ~ 5°, which is the highest among the ensembles, and therefore removed from the final list. After the exclusion of models based on the above-mentioned criteria, we are left with 10 best-performing models out of the 27 model ensemble, and they are CMCC-CM2-SR5, HadGEM-GC31-LL, IPSL-CM6A-LR, CNRM-CM6-1, CNRM-ESM2-1, NorESM2-MM, CMCC-CM2-HR4, BCCCSM2-MR, NorESM2-LM and EC-Earth3.
In order to calculate the projected changes over the Indian ocean for mid-emission (SSP2-4.5) and high-emission (SSP5-8.5) radiative forcing scenario, the averaged departure of projected variable for 2080-2100 relative to the historical  (Fig. 14). Interestingly, the ensemble of the best-performing models shows a weaker sea level rise in the IO compared to the full MMM solution (Fig. 15). In NIO, the Arabian Sea shows a stronger rise in sea level compared to the Bay of Bengal. The best model estimate shows a sea level rise of ~ 22-25 cm under the SSP2-4.5 scenario and ~ 35 cm rise for the SSP5-8.5 scenario. The western Bay of Bengal also shows a similar rising trend as in the Arabian Sea, but indicates a much weaker rise of DSL in the eastern part. While this increase in sea level is primarily driven by the rise in the global mean thermosteric change (~ 25 cm, ~ 85%, Fig.  S3), ~ 12-15% is contributed by the weakening of the summer monsoon wind field (Fig. 14b, e). The SIO also shows a strong sea level rise along the ACC current regime. There is a poleward shift of the Southern Ocean westerlies consistent with other studies (Kidston and Gerber 2010;Lyu et al. 2020a, b) due to the warming of the Southern Ocean driven by anthropogenic forcing. This drives the increase in sea level in this region. The warming of the Southern Ocean also has links to the weakening of the Atlantic meridional overturning circulation forced by aerosol changes under future emission scenarios (Shi et al. 2018). Under the SSP5-8.5 scenario, the maximum sea level rise in this region is expected to be ~ 35-40 cm. In contrast, the STIO latitude band show the minimum sea level rise driven by the weakening of the trade winds. Note that the projected DSL changes (without the global thermosteric effect) are negative in the southern tropical Indian Ocean (20°S-40°S) and the eastern equatorial Indian ocean near the Sumatran coast (Fig. 14a,  d). Note, however, that recently Jevrejeva et al. (2020) analysed zostoga from CMIP5 and CMIP6 models and found out that CMIP6 models only show half of the observed rate of change during the period 1940-2005 and hence caution should be exercised in the interpretation of results.
Like the mean field, the sea level variability also shows a similar pattern in the projected change under the SSP2-4.5 and SSP5-8.5 scenarios (Fig. 16a, c). However, note that, in contrast to the mean field, the best estimate of the ensemble mean variability is generally higher in the entire basin compared to the mean field (Fig. 16b, d). The western Arabian sea shows the strongest increase in the sea level variability with an increase of ~ 30% from the base period under the high emission scenario. This is a significant increase over the CMIP5 model estimates on the projected variability of this region (see Fig. 8e of Deepa et al. 2021). The ACC current regime of the SIO and the west coast of Australia also show a considerable increase in the sea level variability. The thermocline ridge region, however, shows a marked decrease in variability, which is in agreement with the CMIP5 models (Deepa et al. 2021). Note, however, that the thermocline ridge region is shown to be one of the internal variability dominated regions (Han et al. 2014). Hence, considering that the CMIP6 models underestimate the variability of this region (Fig. 7), the projected decreasing variability of this region should be treated cautiously. On the other hand, unlike the CMIP5 model, CMIP6 models project a stronger increase in sea level variability in the eastern equatorial Indian Ocean and along the coast of Java, suggesting the influence of an increase in effective climate sensitivity of the CMIP6 models as reported earlier by Hermans et al. (2021). Further, the overestimation of the dipole strength in the CMIP6 modes is also likely to impact the long-term climate signatures of this tropical basin.

Conclusion
In this study, we evaluate the skill of models participating in the sixth phase of coupled model intercomparison project (CMIP6) in simulating the DSL and the associated wind field over the Indian Ocean. Subsequently, we select a subset of models by discarding grossly biased models from the ensemble to get the best-estimated projections over the Indian Ocean. Models' ability to reproduce the mean states, variability and climate modes for the observational period is used to assess the skills of each model. We use statistical tools such as pattern correlation, RMSE, and skill score metrics to assess the performance of models. Projections of sea level and wind field are analyzed for moderate (SSP2-4.5) and high (SSP5-8.5) radiation forcing future scenarios.
Most models could reproduce the mean DSL of the IO very well with a skill score of more than 0.95. However, the skill score varies considerably across the latitude band with the lowest scores exhibited in the NIO. All the models show consistent biases in sea level and winds in a few dynamically dominated regions. Nevertheless, the zonally averaged DSL show consistent positive bias across all latitude with maximum bias in the ACC regime. Also, the equatorward Fig. 16 a Projected change in the variability of MMM DSL SSP2-4.5 scenerio. b Difference between the ensemble mean of variabilites from the best performing models and the MMM for SSP2-4.5 scenerio. c, d same as (a, b), but for the SSP5-8.5 scenerio shift of the Southern Ocean westerlies is found to be more in the Indian Ocean compared to the global averaged. In the NIO, the weaker summer monsoon winds in the model cause weak coastal upwelling along the western boundary of the Arabian Sea, resulting in a positive bias in the western IO.
The equatorward shift of the southeasterly trade winds causes anomalous easterly winds along the equator in most models. These wind biases shoal the thermocline of the eastern equatorial IO and reduce the west-east sea level (thermocline) tilt which is otherwise found in the observations. This relatively flat equatorial thermocline in the east enhances the air-sea feedback and causes a dipole pattern in SST and sea level in the tropical basin. This IOD like bias is also reflected in the EOF analysis with stronger than observed variability in the first EOF mode. Some of the models such as NorESM2-MM, GISS-E2-1-G and the models from the EC-Earth3 family exhibit almost three times of the observed variability, indicating a strong bias towards the IOD mode.
In the case of variability, model simulations generally underestimate relative to the observation across the entire basin except for the eastern equatorial IO and the Sumatran coast. The overestimated magnitude of variability along the Sumatran coast is likely due to the shallow thermocline driven by the easterly bias in the equatorial zonal winds. As expected, the high-resolution models reproduce the observed variability better compared to the coarser models, particularly in the high eddy-dominated regions like ACC current regime.
We assess the fidelity of the model simulations in reproducing the local and remotely forced sea level variability in the Indian Ocean. A large inter-model spread is observed in the correlation values between the forcing and the sea level variability over the TRIO region. This large spread might be caused by the deficiencies in simulating the inter-annual winds over the eastern equatorial Indian ocean and the ITF influx, both of which remotely influences the sea level over TRIO region. The simulated influence of the Pacific Ocean over the SEIO is weaker than observed with a lower correlation between the western equatorial pacific winds and the sea level over SEIO. Longer lag than the observation corresponding to maximum correlations indicates that the Pacific Ocean influence is not properly simulated in most models. In the Southern Ocean, DSL biases are more pronounced over the Indian ocean due to the equatorward shift of the Southern Hemisphere westerlies compared to other ocean basins. The strong bias in the windstress curl and the bottom pressure torque also contribute to the DSL bias in the Southern Ocean.
The projected sea-level change towards the end of the century is calculated as the difference between the sea level averaged over the last two decades of the projection (2080-2100) and historical simulations . Both scenarios show a similar spatial pattern of sea-level change with differences in magnitude. The spatial pattern of the sea level change shows a west-to-east gradient in the northern Indian Ocean with a significant rise in the Arabian Sea and along the east coast of India. This projected sea level rise in the Arabian Sea is primarily due to the warming of the water column and weakened summer monsoon winds. The STIO, in contrast, shows an overall dip in the projected sea level. In terms of variability, the western Arabian Sea show a 30% increase in projected sea level variability. In the SIO, the ACC current regime also shows a marked projected increase but shifted poleward as the westerly winds are also projected to move poleward owing to the anthropogenic warming of the Southern Ocean.
Finally, the best estimate of sea level projection based on an ensemble of 10 best-performing climate models indicates a slightly weaker rise in projected sea level compared to the MMM ensemble but shows higher projected variability in the entire basin.