This study evaluated the historical experiments from 1850–2014 of 17 CMIP6 ESMs (Table 1). We select these ESMs because they provided full periods of the following three standard CMIP6 experiments: piControl, historical, and SSP5-8.5 experiments (Bourgeois et al., 2022). Although this study focuses on evaluating historical experiments, the assessments of these ESMs can provide benchmark results for future sensitivity studies of climate change. Table 1 shows an overview of the ESMs used in this study.

Since this study aims to evaluate the ESMs performance in the STZ, the spatial domain is confined within 30°S-55°S. Previous studies have documented the important roles of the atmosphere and oceans in the STZ in the climate system. Westerlies in the Southern Ocean exert wind stress on the sea surface and drive the ACC and downwelling of surface water in the STZ (Rintoul et al., 2001; Rintoul and Garabato, 2013). As surface waters are transported northward, the subduction of mode water and intermediate water can greatly contribute to the carbon sink in the STZ (Gruber et al., 2019; DeVries et al., 2017). Meanwhile, the carbon uptake in the STZ is also modulated by the SST and overlying atmospheric temperatures via air-sea heat fluxes (Frölicher et al., 2015). Overall, air-sea interactions in the STZ, including heat and freshwater fluxes, can affect the ocean stratification and water masses formation, with implications for the global climate system (IPCC, 2021). Therefore, we evaluate and compare the performance of 17 CMIP6 ESMs in terms of the zonal and meridional 10 m winds (u_10m, v_10m), 2 m temperature (t_2m), 2 m specific humidity (q_2m), precipitation (P), surface downwelling shortwave radiation (rsds), surface downwelling longwave radiation (rlds), and the oceanic potential temperature (θ) and salinity (S).

The analysis in this study utilizes a single ensemble member (r1i1p1f1 or equivalent) of the CMIP6 historical experiments. For the ESMs that have not provided the r1i1p1f1 ensemble member, equivalent ensemble members are selected.

To evaluate the simulation performance of the CMIP6 ESMs, we adopt the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF; ERA5) atmospheric reanalysis data set to evaluate the atmospheric fields (Hersbach et al., 2020, 2023), and we adopt the objective analysis data set of the World Ocean Atlas 2023 (WOA23) to evaluate θ and S (Reagan et al., 2023).

The performance of the ERA5 reanalysis has been comprehensively evaluated by Hersbach et al. (2020). Produced by using the 4D-Var data assimilation and the ECMWF Integrated Forecasting System (IFS), the atmospheric model is coupled to a land-surface model and an ocean wave model. With a spatial resolution of 31 km and hourly outputs, ERA5 provides comprehensive atmospheric data on 37 vertical pressure levels. The ERA5 provides the reanalysis data from 1940 to the present. To assess the CMIP6 performance in the STZ, we use the ‘Surface or single level’ data category, which has been provided as 2D parameters, including u_10m, v_10m, t_2m, q_2m, P, rsds, and rlds.

The WOA23 provides climatological annual means and monthly climatology of in situ temperature and S. Building upon the fundamental framework of the Climatological Atlas of the World Ocean (Levitus, 1983), WOA23 is the latest advancement of these oceanographic climatological analyses. Compared to the last version WOA18 that is based on oceanographic casts during 1955-2017, WOA23 incorporates more hydrographic observations during 1955–2022 from the World Ocean Database 2023 (Mishonov et al., 2024). World Ocean Database 2023 includes in situ measurements from ships, autonomous floats and gliders (e.g., Argo program), and moored buoys. Based on an objective analysis of observations, WOA23 provides a climatological analysis of in situ temperature and salinity on the 0.25° × 0.25° horizontal grids, with 102 vertical layers ranging from 5 m at the sea surface to 100 m at 5500 m depth. Note that the climatological annual mean of WOA23 provides data with the full 5500 m depth range, while monthly climatology data are only provided in the upper 1500 m layers. In this study, we employ the climatological mean data with the label of ‘all’, an average of all available data, on the 0.25° × 0.25° grid resolution.

To evaluate the ESMs representation in the STZ, the climatological annual mean of every ESM is calculated, and monthly climatology data are also calculated for the evaluations across seasons. Generally, the seasonal cycle in the Southern Hemisphere is defined as the austral spring (September, October, and November), the austral summer (December, January, and February), the austral autumn (March, April, and May), and the austral winter (June, July, and August), respectively.

To show the inter-model spread, the standard deviation (SD) of climatological annual means across the 17 CMIP6 ESMs was calculated as follows (Huang et al., 2020):

where Pi is the climatological annual mean of an ensemble from individual ESM, PMME is the multi-model ensemble (MME) mean, and N=17 is the number of ESMs employed. The SDESMEquation 1 indicates the dispersion of CMIP6 ESMs relative to their MME.

To evaluate the model abilities of ESMs, we adopt the MVIE method. The development of the MVIE method undergoes three phases. First, Xu et al. (2016) devised the vector field evaluation (VFE) diagram, which is a generalized Taylor diagram (Taylor, 2001). The VFE diagram quantifies model skill in simulating vector fields through three statistical metrics: the root-mean-square length (RMSL), the vector similarity coefficient (VSC), and the root-mean-square vector difference (RMSVD). The RMSL can measure the magnitude and variance of vector lengths, the VSC can assess the pattern similarity of normalized vector pairs, and the RMSVD can represent overall deviations. The RMSL, VSC, and RMSVD are calculated as follows:

where A denotes a vector field that can be written as a pair of vector sequences; In Equations 2–4, Ai=(xi,yi), i = 1, 2,…, N. N means the number of vectors in the sequence.

where B denotes a reference vector field, and the • symbol denotes the inner product. The value of the VSC ranges from -1 and 1, with the larger value corresponding to the higher similarity between vector fields.

where the RMSVD becomes smaller when two vector fields approach more alike in both direction and vector length. The VFE diagram overcomes the limitations of scalar-based Taylor diagrams, enabling the assessments of both the directions and amplitudes of vector fields.

Second, Xu et al. (2017) have proposed the MVIE method. By normalizing and grouping scalar fields into a multidimensional vector, Xu et al. (2017) introduced the multivariable integrated evaluation index (MIEI) to summarize the model performance across simulated fields. For multiple scalar fields, the MIEI is calculated as follows:

where M is the number of scalar fields, LAm* in Equation 5 denotes the ratio of the root mean square (RMS) of the m-th standardized scalar field to the reference value (Equation 6). When the RMS of the simulated scalar fields approaches the RMS of the reference data, LAm* approaches 1. F is a weighting factor that adjusts the relative importance of data amplitude and data similarity in the MIEI. When F>1, the MIEI is more sensitive to the changes of data similarity than the changes of amplitude, and vice versus. In this study, we adopt F=2 that was proposed by Xu et al. (2017) and used in Han et al. (2022).

The MIEI combined the RMSL deviations and VSC into a single metric, taking the pattern similarity of multiple scalar fields and amplitude into account. Additionally, it fulfills the requirement that a model performance index should have the monotonic property with respect to the performance of ESMs. The MIEI can comprehensively reflect the overall performance of the simulated multiple scalar fields. A smaller value of the MIEI denotes a better performance of a model in simulating multiple scalar fields.

Third, Zhang et al. (2021) further improved the MVIE method by incorporating the area-weighted statistics and the combination of multiple scalar and vector fields. The area-weighted RMSL, VSC, and RMSVD are reintroduced as follows (Equations 7–9):

where wi is an area-weighting factor and the sum of wi is equal to 1. The area-weighting statistics are more accurate for a global evaluation of ESMs (Zhang et al., 2021). Since the first term on the RHS of Equation 5 can vary from 0 to + ∞, while the second term on the RHS of Equation 5 only ranges from 0 to 4, the MIEI may be somehow too sensitive to the RMS bias rather than the pattern bias. To further improve the representation of the evaluation index, the multivariable integrated skill score (MISS) was proposed as a normalized index (Zhang et al., 2021). The MISS is a flexible index that can adjust the relative importance of the pattern similarity and amplitude errors, which is defined as follows:

where Rm* in Equation 10 is a piecewise function determined by the value of LAm* (Equation 11). The MISS varies monotonically with the overall performance of ESMs, typically ranging from 0 to 1. If ESMs results are same as the reference values, the MISS is equal to 1.

Zhang et al. (2021) further proposed a centered mode of the above statistical metrics. Uncentered statistical metrics are calculated using the original field, while centered statistical metrics are calculated by the anomaly field generated by removing the spatial mean from each grid point of the original field. The statistical metrics in the centered mode provide insights into the evaluation of the multivariable anomalous fields. The centered RMSL, VSC, RMSVD, and the vector mean error (VME) are introduced as follows (Equations 12–15):

Indeed, the framework of the MVIE can be introduced to the evaluation of both scalar and vector fields (Han et al., 2022). In this study, we combine the multiple normalized scalar and vector fields into a single vector field for every selected ESM and reference data, respectively. Then, the VFE diagram is used to compare the combined vector fields from ESMs with reference data.

To show more details of the performance of the ESMs, we also use several statistical metrics to reveal the model ability in simulating a single scalar or vector field. For the first level of scalar and vector fields, we use the mean error (ME), root-mean-square difference (RMSD), correlation coefficient (CORR), and standard deviation (SD) to evaluate the model performance. For the second level of a multiple dimension vector field, which is composed of group scalar and vector fields, we use the VME, VSC, RMSL, and RMSVD to show the model performance. For the third level, we use the MIEI and MISS to provide a synthesized evaluation index for the model performance. Note that the SD and ME are normalized by dividing by the SD of the reference data.

In previous studies, this MVIE method has been employed in evaluations of ESMs. For example, Huang et al. (2019) used the VFE diagram and statistical quantities devised by Xu et al. (2016) to evaluate vector winds in the Asian-Australian monsoon region simulated by CMIP5 models. Building on the MVIE framework, Lv et al. (2020) assessed the overall performance of the Weather Research and Forecasting (WRF) model with various physics schemes in simulating precipitation and soil moisture over the central Tibetan Plateau. Han et al. (2022) evaluated the performance of CMIP6 models in simulating the large-scale environmental fields of tropical cyclones in the low and middle latitudes. Dai et al. (2021) diagnosed the influences of different parametrization schemes in the WRF model on the precipitation and temperature in northern China. Zhang et al. (2022) evaluated and ranked the ability of CMIP6 ESMs over coordinated regional downscaling experiment domains.

To our knowledge, the performance of the CMIP6 ESMs in the STZ has not been comprehensively evaluated in term of multiple variables, and thus we tend to assess the CMIP6 ESMs in the STZ in this study. The θ and S fields are divided into three layers, including the surface layer, the upper 1500 m layer, and the layer below 1500 m. Then, depth-averaged θ and S are calculated to derive 2D scalar fields. To assess the ability of ESMs on a uniform grid, we interpolate all the fields onto a 0.25° × 0.25° grid mesh with the bilinear interpolation method. Bilinear interpolation is a statistical method and widely used in model evaluation (e.g., Han et al., 2022; Qiu et al., 2024; Talukder et al., 2025). It is especially suitable for continuous variables such as θ, S, and winds. And this method preserves spatial gradients reasonably well while avoiding the excessive smoothing or artefacts from higher-order interpolation methods. Note that the in-situ temperature from WOA23 is converted to θ as the reference values.

Before conducting the multivariable assessment, we show the climatological annual mean of the reference fields in STZ (Figure 1). In the STZ, the u_10m is generally weak around 30°S (Figure 1A), mainly dominated by trade winds (Spiridonov and Ćurić, 2021), with the wind direction from east to west. Between 35°S and 55°S, the westerlies are prevailing and strong. The v_10m is dominated by north wind over higher latitudes, while strong south winds occur at the western boundaries of the mainland (Figure 1B). The t_2m in the STZ has a remarkable meridional gradient (Figure 1C), with a gradual decrease as the latitude increases. The climatological annual mean of t_2m typically ranges from 10 °C to 20 °C north of 50°S, whereas it cools significantly south of 50°S. The q_2m decreases with increasing latitude in the STZ, with a relatively large band to the east of the mainland (Figure 1D). The P is generally low in the Southern Ocean (Figure 1E), while the strong precipitation occurs at the western boundaries of New Zealand and South America, due to the influence of westerlies and terrain features (Garreaud et al., 2009). The rsds and rlds represent the downward solar radiation reaching the Earth surface and the downward longwave radiation, respectively. They influence the surface radiation budget directly (He et al., 2023; Wild et al., 2015). As latitude increases southward, the climatological annual means of both rsds and rlds gradually decrease (Figures 1F, G), with strong radiation at 30°S and weak radiation at higher latitudes.

The climatological mean of θ_surface is higher around 30°S, typically exceeding 20 °C, and gradually decreases southward (Figure 1H). In the 50°S, θ_surface approaches 0 °C. The climatological mean of θ_above1500 ranges from 0 °C to 15 °C in most regions in the STZ (Figure 1I), while θ_below1500 is colder (Figure 1J). The spatial distribution of S_surface and S_above1500 in the STZ also exhibits noticeable meridional gradients (Figures 1K, L). As latitude increases, the climatological annual mean of S_surface and S_above1500 decreases, whereas the spatial gradient of S_below1500 is much weaker (Figure 1M).

Figure 2 shows the difference between the CMIP6 MME climatological mean and the climatological mean of two reference datasets (ERA5 and WOA23). The CMIP6 ESMs show stronger biases in westerlies (Figures 2A, B), with intensified west and north winds in the ESMs. The CMIP6 MME mean overestimates u_10m by 0.5-1.5 m s⁻¹ over the STZ, particularly in the Indian Ocean sector and the southern Australian ocean. Meanwhile, the CMIP6 MME mean overestimates the negative v_10m by 0.5–1 m s⁻¹ in most areas of the STZ. In terms of t_2m, a cold bias of 1-2 °C dominates the simulated 2 m air temperature across most ocean surface at 30-45°S (Figure 2C). At 45-55°S, there is a warm bias of ~1 °C in the Indian and Atlantic Oceans (Figure 2C). In the STZ region, q_2m is generally overestimated by about 1-1.5 g kg^-1 in the CMIP6 simulations, while an underestimation occurs in the South American continent (Figure 2D). In the north of 40°S, the CMIP6 ESMs show a negative bias in the simulated climatological annual mean of P, whereas P generally presents a positive bias in the south of 40°S (Figure 2E). The biases of the rsds are generally opposite to that of the rlds bias (Figures 2F, G). The rsds shows an overestimation in the north of 40°S and an underestimation in the south (Figure 2F), while the rlds shows an opposite spatial distribution (Figure 2G). This result shows reasonable concordance with Xu et al. (2022). Xu et al. (2022) evaluated the simulated rlds by comparing the CMIP6 simulation results with CMIP5 results and ground measurements and identified the positive biases of rlds at 40-55°S.

The oceanic fields reveal pronounced thermal biases: the MME mean underestimates the θ_surface by 1-2 °C at 40-55°S but shows 1-2 °C warm anomalies at 30-40°S (Figure 2H). In terms of θ_above1500, the CMIP6 ESMs overestimate θ by 1-2 °C across most of the STZ ocean (Figure 2I). However, the simulated θ_below1500 is grossly overestimated by CMIP6 MME (Figure 2J). At 30°S, this warm deviation approaches 5 °C (Figure 2J). On the contrary, the simulated climatological annual mean of S_surface and S_above1500 shows a clearly fresh bias (Figures 2K, L). In terms of S_below1500, the fresh bias between the CMIP6 MME mean and reference data is still very large (Figure 2M).

These biases underscore the substantial discrepancies in the CMIP6 MME mean across variables. Indeed, the magnitude of biases between the CMIP6 MME mean and reference data may remarkably depend on the region analyzed.

The SD_ESM (Equation 1) can serve as a metric to quantify the dispersion of CMIP6 ESMs outputs relative to their MME mean. The inter-model spread, quantified by the SD_ESM across CMIP6 simulations (Figure 3), highlights significant uncertainty in the near-surface atmospheric fields and the oceanic θ and S. The u_10m shows maximum variability over the ACC region, reflecting divergent simulations of westerly jet intensity (Figure 3A). The v_10m and t_2m fields exhibit similar patterns with respect to the inter-model spread (Figures 3B, C), and their SD_ESM are generally smaller than 1 m s^-1 and 1.5 °C across most regions in the STZ, respectively. In terms of q_2m, CMIP6 ESMs show a greater inter-model spread at 30-40°S, especially in the coastal region (Figure 3D). Relatively larger inter-model spread of P is identified at the western Andes in South America and the western side of the New Zealand Island (Figure 3E). Compared to other variables, the rsds and rlds fields exhibit consistently higher inter-model SD_ESM values in CMIP6 ESMs (Figures 3F, G), underscoring systemic challenges in radiative fluxes. The inter-model spread of θ_surface is less than 2 °C (Figure 3H), except for the region of the Brazil Current on the eastern side of South America. In terms of the inter-model spread of θ_above1500 and θ_below1500, the CMIP6 ESMs exhibit similar spatial patterns, with a greater spread in 30-45°S, and the Southern Indian Ocean, and a relatively small spread in 45-55°S (Figures 3I, J). The distribution of the SD_ESM of S exhibits distinct spatial patterns in different layers: larger values of the SD_ESM are concentrated in subtropical regions (30-40°S) at the surface layer (Figure 3K), while there is a relatively larger inter-model spread over 45-55°S below 1500 m depths (Figure 3M). At the upper 1500 m depth, the SD_ESM of S_above1500 is relatively small (Figure 3L). The maximums of the SD_ESM of S_above1500 are still around the South American continent (Figures 3K-M).

Based on the inter-model spread analysis of CMIP6 ESMs, significant discrepancies still exist among simulations of these key fields. These inter-model spreads underscore persistent systematic challenges in representing the near-surface atmospheric fields, sea fields, and air-sea coupling interactions in the ESMs.

As discussed in sections 3.1 and 3.2, there are notable discrepancies between the CMIP6 ESMs and the reference data, and the inter-model spread is identified in the simulated multiple variables from the CMIP6. To further analyze the simulation ability of the CMIP6 ESMs, in this section, we employ the MVIE method to systematically evaluate the individual CMIP6 ESMs in simulating multiple variables. Our evaluation focuses on the nine variables specified in section 2.1, encompassing both near-surface atmospheric fields and three-dimensional oceanic fields (denoted by θ and S fields at varying depths).

The VFE diagram uses the RMSL, RMSVD, and VSC to offer comprehensive statistics on the abilities of the ESMs, including the differences between various ESMs and the differences between ESMs and reference data. Figure 4 provides a straightforward intercomparison of 17 CMIP6 ESMs by evaluating their simulated climatological annual mean vector fields that represent all variables.

In the VFE diagram, the model ability can be diagnosed by the RMSVD between the model and reference data, with a smaller distance from the REF (the black reference point in Figure 4) suggesting a good representation of an ESM. A smaller RMSVD is usually associated with a higher VSC and an RMSL approaching 1. The statistics, such as RMSVD, VSC, and RMSL, of 10 ESMs (the blue names in Figure 4) are close to the reference data, indicating that the differences between these 10 ESMs products and the reference data are relatively small over the STZ region. The VSCs of the 10 ESMs (the blue names in Figure 4) range from 0.94 to 0.99, suggesting that these ESMs can reproduce the spatial pattern of the near-surface atmospheric fields and the oceanic θ and S, with the best estimation from the ACCESS-ESM1-5. The normalized RMSLs of the 10 ESMs are generally larger than 1, indicating that these ESMs tend to overestimate the simulated climatological annual mean vector fields, with the best estimation from the INM-CM4-8.

In contrast, the VSCs of the rest 7 ESMs (the black names in Figure 4) range from 0.6 to 0.85, implying great differences in the spatial patterns of the climatological annual mean vector field between the ESMs products and the reference data. In addition, these 7 ESMs also greatly overestimate the amplitude of the climatological annual mean vector field. Furthermore, we also used the centered VFE diagram to evaluate the anomaly vector fields (Supplementary Figure S2), and the results of model evaluation are similar to those derived from the climatological annual mean vector fields.

The VFE diagram provides insights into the model abilities of ESMs by the VSC and RMSL. To show the performance of ESMs in detail, we further assess more statistical metrics (Figure 5). In addition, we calculated the cMISS to show the multivariable integrated skill, which takes the cVSC and SD into account together. The metrics table adopted the centered statistics that decompose the original fields into anomaly and mean fields. Evaluations of anomaly fields are conducted from three perspectives: the variance characteristics (SD, cRMSL), the spatial pattern consistency (CORR, cVSC), and the root-mean-square differences between the reference data and ESMs results (cRMSD, cRMSVD). The ME is also incorporated in Figure 5 to show the systematic biases of the simulated original fields of ESMs from the reference data.

In terms of the simulated spatial pattern of t_2m, q_2m, rsds, rlds, and θ_surface, the CMIP6 ESMs generally show relatively good performance, with CORRs higher than 0.9. In contrast, with CORRs ranging from 0.1 to 0.9, the CMIP6 ESMs cannot match closely with the spatial pattern of P, θ_below1500 and S from the reference data. It is not easy for these ESMs to adequately reproduce the spatial pattern of θ and S in the deep ocean. Nonetheless, some models still perform relatively well, including CanESM5, CanESM5-CanOE, GFDL-ESM4, and MRI-ESM2-0, with CORRs larger than 0.8. Most ESMs can capture the spatial patterns of the u_10m and v_10m, yet CORRs of winds from four models are lower than 0.9, including MIROC-ES2L, MRI-ESM2-0, NorESM2-LM, and INM-CM4-8. To derive a comprehensive estimation of the ESMs in simulating the spatial pattern of fields in the STZ, we calculated the cVSC to evaluate the overall performance. The GFDL-ESM4 model shows the highest cVSC (~0.939) among the 17 CMIP6 ESMs, indicating that this model is the most consistent with the reference data in simulating the spatial pattern of the near-surface atmospheric fields and the oceanic θ and S in the STZ region.

The CMIP6 ESMs exhibit considerable differences in simulating the spatial SD of the different fields. For instance, most ESMs tend to overestimate the spatial variability of 10 m vector wind, with 11 out of 17 CMIP6 ESMs overestimating both u_10m and v_10m over the STZ region. In terms of the spatial variability of q_2m, most ESMs tend to have an overestimation by 1%-35%, while the MIROC-ES2L, MPI-ESM1-2-LR, and INM-CM4–8 models tend to have an underestimation by 3%-9%. In terms of the spatial variability of θ_below1500 and S_below1500, characterized by an SD larger than 2, some CMIP6 ESMs have a significant overestimation, including the ACCESS-ESM1-5, CESM2, CESM2-WACCM, GFDL-CM4, IPSL-CM6A-LR, MPI-ESM1-2-LR, MPI-ESM1-2-HR, NorESM2-LM, NorESM2-MM, UKESM1-0-LL, and INM-CM4–8 models. In terms of t_2m, P, rlds, θ_surface, and S_surface, most CMIP6 ESMs overestimate the spatial variability, whereas the simulated rsds tends to be systematically underestimated. Yet, the SD values of these fields are smaller than 2. As the total SD across all selected fields, the value of the normalized cRMSL larger (smaller) than 1 denotes that the ESM overestimates (underestimates) the anomaly field’s amplitude error. Most ESMs have overestimations, whereas the CNRM-ESM2–1 and MIROC-ES2L models have underestimations. Although the INM-CM4–8 model overestimates the spatial variability of S_below1500, it is in most agreement with the reference data, with the cRMSL approaching 1.

The CMIP6 ESMs also show remarkable diversity in simulating the ME of different variables, with ME ranging from -9.2 to 19.9 (Figure 5). Apart from IPSL-CM6A-LR and MIROC-ES2L, with the ME ranging from -0.22 to -0.02, most CMIP6 ESMs overestimate u_10m over the STZ. These stronger biases are consistent with the differences between the MME mean and the reference data (Figure 2A). Conversely, apart from the IPSL-CM6A-LR model, most CMIP6 ESMs underestimate v_10m, which is in agreement with the broad negative values in Figure 2B. For t_2m, the ME of CESM2 has the minimal absolute value that approaches 0, whereas the MIROC-ES2L model has the largest maximum up to 0.297. The INM-CM4–8 model captures the minimum ME of P, while the MIROC-ES2L model tends to overestimate the ME of P mostly. In terms of q_2m, the CMIP6 ESMs all tend to have overestimation as indicated by the positive ME. Indeed, the ME values of the near-surface atmospheric fields are all less than 1, implying the relatively good agreement with the reference data. In contrast, the θ and S in the CMIP6 ESMs show relatively larger biases in the ME, especially in the θ_below1500 and S_below1500. The θ_surface and S_surface still match well with the reference data, except for the INM-CM4–8 model with the ME of S up to -1.188. In terms of θ_above1500 and S_above1500, the ME values of the CESM2, CESM2-WACCM, and INM-CM4–8 models are larger than 1. Twelve ESMs show strong biases in simulating the θ_below1500, whereas the ACCESS-ESM1-5, MRI-ESM2-0, NorESM2-LM, NorESM2-MM, and INM-CM4–8 models still can have the ME values of θ_below1500 less than 1. Similarly, thirteen ESMs show strong biases in simulating the S_below1500, whereas the ACCESS-ESM1-5, CanESM5, CanESM5-CanOE, and CNRM-ESM2–1 models can be close to the S_below1500 of the reference data, with the ME values less than 1. The VME measures the differences between two vector fields, and the ACCESS-ESM1-5, CanESM5, CanESM5-CanOE, CNRM-ESM2-1, GFDL-ESM4, NorESM2-LM, NorESM2-MM, and INM-CM4–8 models have relatively smaller VME, with the values less than 0.5.

The statistics of cRMSD are similar to those of ME (Figure 5). The cRMSD values of the near-surface atmospheric fields are also less than 1, indicating statistical agreement with the reference data. The θ_surface, θ_above1500, S_surface, and S_above1500 still match well with the reference data, with the cRMSD values all less than 1. Akin to the ME values, the CMIP6 ESMs exhibit larger values of the cRMSD for the oceanic θ_below1500 and S_below1500. Most ESMs show strong differences in anomaly fields between the simulations and the reference data in the abyssal ocean, whereas the cRMSD values of θ_below1500 and S_below1500 of the CanESM5, CanESM5-CanOE, CNRM-ESM2-1, and MRI-ESM2–0 models are still less than 1. The cRMSD of multiple fields is measured by the cRMSVD, which indicates the overall difference of the anomaly field in terms of the near-surface atmospheric fields and the oceanic θ and S. Among the CMIP6 ESMs, GFDL-ESM4 has the minimal cRMSVD (0.367), indicating the smallest overall error of multiple anomaly fields.

On the whole, there is no CMIP6 ESM that performs best in every simulated field. The cRMSVD provides an overall evaluation of the model performance, with a smaller RMSVD value corresponding to a better consistency between the CMIP6 ESM and the reference data. However, improvements in the model performance may not always be associated with a monotonically decreasing RMSVD (Huang et al., 2019; Xu et al., 2017). Therefore, the values of the cMIEI and cMISS are computed to provide an overall evaluation that is monotonically associated with the model performance (Figure 5).

In order to assess the overall simulation capability for the near-surface atmospheric fields and the oceanic θ and S, we further focus on the values of the centered MIEI and MISS of the 17 CMIP6 ESMs in the STZ region (Figure 6). A better model performance is indicated by a smaller value of the cMIEI and a larger value of the cMISS. We define two benchmark thresholds for a quantitative evaluation: (i) the ESMs with cMIEI < 1 exhibit better simulation skill with respect to the climatology annual mean vector field (including the ACCESS-ESM1-5, CanESM5, CanESM5-CanOE, CNRM-ESM2-1, GFDL-ESM4, MIROC-ES2L, MRI-ESM2-0, NorESM2-LM, NorESM2-MM, and INM-CM4-8); (ii) the ESMs with cMISS > 0.9 exhibit better simulation skill (including the ACCESS-ESM1-5, CanESM5, CanESM5-CanOE, CNRM-ESM2-1, GFDL-ESM4, MRI-ESM2-0, NorESM2-LM, and NorESM2-MM).

The cMISS is an advance over the cMIEI in the reduced sensitivity to amplitude errors, yet we still provide the values of cMIEI for a reference. Based on these two indices, the model performance of the CMIP6 ESMs evaluated is generally consistent, except for the MIROC-ES2L and INM-CM4-8. Finally, eight ESMs (including the ACCESS-ESM1-5, CanESM5, CanESM5-CanOE, CNRM-ESM2-1, GFDL-ESM4, MRI-ESM2-0, NorESM2-LM, and NorESM2-MM) satisfy both criteria, suggesting their better representation of the near-surface atmospheric fields and the oceanic θ and S in the STZ.

Notable discrepancies of climatological annual means between the CMIP6 ESMs and reference data have been discussed above. Since there is a strong seasonality in the near-surface atmospheric fields and the oceanic θ and S in the STZ (Supplementary Figure S3), a quantitative evaluation of the simulated seasonal climatology could also provide insights into the ESMs ability. Based on the classification of four seasons (DJF, MAM, JJA, SON) described in section 2.2, we further compare the CMIP6 ESMs with the reference data across seasons. Note that our assessments of the seasonality exclude the oceanic θ and S below 1500 m depth because the monthly climatology data of WOA23 only provides data in the upper 1500 m layers.

To compare the performance of various CMIP6 ESMs in reproducing multivariable fields in different seasons, Figure 7 illustrates the VFE diagram of the climatological seasonal means with multiple statistics. The VSCs of all ESMs are larger than 0.95, suggesting that these ESMs can properly replicate the spatial patterns for near-surface atmospheric fields and oceanic θ and S in different seasons. In contrast, the normalized RMSLs generally exceed 1 across seasons, indicating that most ESMs tend to overestimate the vector fields of climatological seasonal means. According to the differences of the RMSL between the ESMs and reference data, the ESMs have better representation in the austral summer (13 ESMs with blue names) and lowest model ability in the austral winter (4 ESMs with blue names).

In stark contrast to the VSCs of the climatological annual mean (Figure 4), the VSCs are mostly improved in the ESMs seasonal products (Figure 7), implying better representations of the spatial patterns of the ESMs seasonal products. The better representation of the seasonal evaluation should be attributed to the exclusion of the oceanic θ and S in the deeper layer, indicating that the oceanic θ and S in the deep layer have larger uncertainties than in the upper layer in the CMIP6 ESMs. We have used the centered VFE diagram to evaluate the anomaly vector fields of seasonal climatology (Supplementary Figure S4), with the outcomes closely resembling those obtained from the uncentered VFE diagram.

The VFE diagrams preliminarily show the model abilities of ESMs through the VSC and RMSL metrics across all four seasons. To delineate the capacity of these models in detail, the 17 CMIP6 ESMs are compared with the reference data. We conduct the comparison by evaluating the VME (ME), cRMSVD (cRMSD), cRMSL (SD), and cVSC (CORR) of the seasonal climatology of the near-surface atmospheric fields and the oceanic θ and S (Figures 8, 9; see Supplementary Figures S5-S8 in the Supplementary Material for the quantitative values in detail). Consistent with the annual assessment, these metrics also employ the centered statistics decomposing seasonal fields into anomaly and mean components.

Considerable differences have been identified among the CMIP6 ESMs in the simulated spatial SD of different fields and seasons (Figure 8). Across all seasons in the STZ region, a larger number of models, 10 of 17 CMIP6 ESMs, tend to overestimate the variability of 10 m winds. In terms of the spatial variability of q_2m, most ESMs also tend to have an overestimation by 2%-45% in the austral summer, autumn, and winter. The CanESM5 and CanESM5-CanOE models both exhibit overestimations up to 45% in winter. In contrast, most CMIP6 ESMs have a better representation of q_2m in spring. Most CMIP6 ESMs underestimate the spatial variability of t_2m, rlds, θ_surface, and S_surface fields across all seasons, with the lowest estimation of S_surface by the INM-CM4-8, while the simulation of rsds field tends to be systematically overestimated. Most CMIP6 ESMs can properly capture the spatial SD of P in different seasons. In terms of the spatial variability of θ_above1500 and S_above1500, the CESM2 and CESM2-WACCM models exhibit significant overestimation across all seasons. In contrast, the MPI-ESM1-2-LR and INM-CM4–8 models show a significant underestimation in simulating the S_above1500. With cRMSL larger than 1, most ESMs have overestimation across all seasons, whereas the IPSL-CM6A-LR models have underestimations throughout the year. Since the CESM2 and CESM2-WACCM models have relatively larger overestimations of the spatial variability of u_10m, v_10m, and θ_above1500, these two models have the largest cRMSL in the austral autumn, with values of 1.392 and 1.409, respectively.

For the simulated spatial pattern of t_2m, q_2m, rsds, rlds, θ_surface, and θ_above1500, the CMIP6 ESMs generally exhibit good performance across four seasons, with CORRs typically exceeding 0.9 (Figure 8). In contrast, in terms of the simulated spatial pattern of P, S_surface, and S_above1500, the CMIP6 ESMs show relatively poorer performance in reproducing the spatial pattern of the reference data, exhibiting CORRs between 0.5 and 0.9. While reproducing oceanic S spatial patterns remains challenging for the CMIP6 ESMs, several models, including the CanESM5, CanESM5-CanOE, GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, and UKESM1-0-LL show relatively good performance, with the CORRs of S exceeding 0.8 across all seasons. While most ESMs can capture the spatial patterns of the u_10m and v_10m during the austral summer, the MIROC-ES2L exhibits a poor wind pattern similarity with the winds of the reference data across all seasons, with the CORRs lower than 0.9. To comprehensively assess the simulated spatial patterns of the ESMs in the STZ, we calculate the cVSC metric for an overall performance evaluation. The MIROC-ES2L model shows the lowest cVSC among the 17 CMIP6 ESMs across four seasons. Most CMIP6 ESMs exhibit relatively poorer performance in simulating the spatial pattern of the near-surface atmospheric fields and the oceanic θ and S during the austral winter, with cVSC generally lower than 0.9.

Remarkable diversity exists among the CMIP6 ESMs in simulating the ME of different variables across different seasons (Figure 9). Most CMIP6 ESMs overestimate u_10m over the STZ across four seasons, except for IPSL-CM6A-LR and MIROC-ES2L. In contrast, most CMIP6 ESMs underestimate v_10m across four seasons, with the lowest value of -0.774 from the MRI-ESM2–0 in spring, except for the IPSL-CM6A-LR, MIROC-ES2L, and INM-CM4–8 models. Most CMIP6 ESMs have a good representation of t_2m and P in simulating the ME across seasons. Yet, the MIROC-ES2L model has a relatively large estimation of the ME of P across all seasons. The ME of q_2m and rlds, are generally overestimated in most CMIP6 ESMs, whereas the rsds in most CMIP6 ESMs shows negative biases in the ME. The ME of θ_surface of all the CMIP6 models shows negative biases in the austral winter. In terms of S_surface, the MPI-ESM1-2-LR, MPI-ESM1-2-HR, and NorESM2-MM models match well with the reference data, with the absolute values of ME less than 0.1 across all seasons. In terms of both S_surface and S_above1500, most CMIP6 ESMs exhibit underestimation as indicated by the negative ME values, while the INM-CM4–8 model shows relatively poor performance with the ME values less than -1 across all seasons. Measuring the differences between two vector fields, the VME metric shows that the MPI-ESM1-2-HR model has the minimum VME value in the austral summer, autumn, and winter, while the UKESM1-0-LL model has the minimum VME value in the austral spring.

The statistics of cRMSD exhibit similar results to the estimation of ME (Figure 9). During the austral summer, autumn, and winter, the values of cRMSD of the near-surface atmospheric fields and the oceanic θ and S are all less than 1, indicating relatively good agreement with the reference data. Yet, the MIROC-ES2L model shows relatively poor performance in simulating the u_10m, v_10m, and P, with the values of cRMSD larger than 0.95 in winter. Among the CMIP6 ESMs, the GFDL-ESM4 model has the minimal cRMSVD (0.242 and 0.308) in the austral summer and autumn. The MPI-ESM1-2-HR model has the minimal cRMSVD (0.361 and 0.472) in the austral spring and winter. These results indicate that these two models have the smallest overall error of multiple anomaly fields in the corresponding season.

Similar to the evaluation of the climatology annual mean vector field, we further calculate the values of the cMIEI and cMISS. These statistical metrics are used to assess the simulation capability of the 17 CMIP6 ESMs for the near-surface atmospheric fields and the oceanic θ and S in the STZ region in different seasons (Figures 9, 10). Most CMIP6 ESMs show good performance in the austral summer and relatively poor performance in the austral winter (Figure 10A). The MIROC-ES2L model has the maximum values of cMIEI across all seasons, implying a relatively poor performance. The evaluation of cMISS is largely consistent with the results of cMIEI (Figure 10B). Based on the cMIEI and cMISS metrics, GFDL-ESM4 shows the best performance during the austral summer and autumn. MPI-ESM1-2-HR and NorESM2-MM perform best in the austral winter, and MPI-ESM1-2-HR leads in the austral spring. In addition, the original vector fields are also analyzed with uncentered statistical metrics (Supplementary Figures S9-11).