An ocean modeling study to quantify wind forcing and oceanic mixing effects on the tropical North Pacific subsurface warm bias in CMIP and OMIP simulations

Sea surface temperature (SST) bias in the climate models has been a focus in the past, but subsurface temperature biases have not received much attention yet. In this study, subsurface temperature biases in the tropical North Pacific (TNP) are investigated by analyzing the CMIP6, CMIP5 and OMIP products, and by performing ocean model simulations. It is found that almost all the CMIP and OMIP simulations have a pronounced subsurface warm bias (SWB) in the northeastern tropical Pacific (NETP), and the model developments over the past decade do not indicate obvious improvements in bias pattern and magnitude from CMIP5 to the latest version CMIP6. This SWB is primarily caused by the model deficiencies in the simulated surface wind stress curl (WSC) in the NETP, which is too weak to produce a sufficient Ekman upwelling, a bias that also exists in OMIP simulations. The uncertainties in the parameterizations of the oceanic vertical mixing processes also make a great contribution, and it is demonstrated that the estimated oceanic vertical diffusivities are overestimated both in the upper boundary layer and the interior in the CMIP and OMIP simulations. The relationships between the SWB and the misrepresented oceanic vertical mixing processes are investigated by conducting several ocean-only experiments, in which the upper boundary layer mixing is modified by reducing the wind stirring effect in the Kraus-Turner type bulk mixed-layer scheme, and the interior mixing is constrained by using the Argo-derived diffusivity. By applying these modifications to oceanic vertical mixing schemes, the SWB is greatly reduced in the NETP. The consequences of this SWB are further analyzed. Because the thermal structure in the NETP can influence the simulations of oceanic circulations and equatorial upper-ocean thermal structure, the large SWB in the CMIP6 models tends to produce a weak equatorward water transport in the subsurface TNP, a weak equatorial upwelling and a warm equatorial upper ocean.


Introduction
Understanding climate model biases is critically important for the assessments of future climate change, and hence "What are the origins and consequences of systematic model biases?" is one of the three key scientific issues that have been addressed in CMIP6 (Eyring et al. 2016). Given the important role played by sea surface temperature (SST) in air-sea coupling, SST biases in climate models have received considerable attention in the past. For example, the too cold tongue bias along the equatorial Pacific can be attributed to the misrepresented oceanic vertical turbulent mixing processes along the equator (Jochum 2009;Sasaki et al. 2013;Furue et al. 2015;Jia et al. 2015;Zhang 2018b, 2019), overestimated cloud albedo in the subtropics (Burls et al. 2017;Thomas and Fedorov 2017), uncertainties in the atmospheric convection parameterizations (Woelfle et al. 2018), and so on. The warm SST bias in the southeastern tropical Pacific and Atlantic can be attributed to coarse resolutions of climate models (Seo et al. 2006;Small et al. 2014), uncertainties in the parameterizations of atmospheric convection and cloud physics (Ma et al. 1996;Gordon et al. 2000), and eastward propagation of equatorial subsurface warm bias (Xu et al. 2014). SST biases can severely degrade the credibility of climate predictions and projections. For example, too cold tongue bias in the equatorial Pacific promotes a La Niña-like warming pattern in the tropical Pacific under the increasing concentrations of greenhouse gases, whereas an El Niño-like warming pattern is produced when the too cold tongue bias is removed from model projections (Li et al. 2016;Ying et al. 2019).
Thermal structure and variability in the subsurface oceans are known to play an important role in the climate change. For example, the acceleration of ocean subsurface warming and the slowdown of surface warming in the early decade of the 21st century indicate that subsurface oceans play a critical role in regulating the global warming (Chen and Tung 2014;Wang et al. 2018). However, compared with that at sea surface, temperature bias in the subsurface oceans is poorly understood. Therefore, it is essential to investigate the oceanic subsurface temperature bias for understanding the simulations and predictions of global energy and heat redistribution. Particularly, subsurface temperature simulation in the tropics is very important to the climate simulations. For instance, Xu et al. (2014) find that the subsurface warm bias within the Atlantic equatorial thermocline can be transported to the Benguela coast by horizontal currents and subsequently surfaces by coastal upwelling, contributing to the longstanding warm SST bias in the eastern tropical Atlantic. In addition, subsurface warm temperature bias along the Pacific equator generates a too diffuse equatorial thermocline, leading to a weak SST-thermocline feedback in the tropical air-sea coupling (Guilyardi et al. 2009;Gao and Zhang, 2017;Zhang et al. 2020). Off the equator, subsurface temperature bias in the climate models is of the opposite sign in the northern and southern tropical Pacific: a large warm bias between 5°-10° N and a large cold bias near 10° S (Wittenberg et al. 2006;Zhu et al. 2020). Though being substantial in its magnitude, the subsurface temperature bias off the equator has not received much attention yet. Importantly, owing to the local Ekman pumping in the Tropical North Pacific (TNP), climatological pycnocline tends to arise to form a potential vorticity barrier, acting to block the local water exchange between the subtropics and tropics (Lu and McCreary 1995;Rothstein et al. 1998;Johnson and McPhaden 1999). In this way, the subsurface warm bias in the TNP acts to weaken the potential vorticity barrier, and would affect the width of subtropical-tropical water exchange window.
Subsurface temperature bias also reflects an erroneous representation of oceanic vertical heat distribution in climate models. Thus, in addition to the poorly simulated atmospheric forcing, flaws in ocean models might be an important source of subsurface temperature bias. One flaw is associated with the oceanic mesoscale eddies, which act to transport heat upward and offset the downward heat transport; these processes are unsolved in coarse resolution models. By increasing the ocean model resolution to the eddy resolving level, subsurface warm bias in the subtropical gyres can be reduced (Griffies et al. 2015;Rackow et al. 2019). Besides, ocean simulations are very sensitive to the parameterizations representing the oceanic vertical turbulent mixing. For example, an overestimated vertical mixing by strong vertical shear is responsible for the subsurface warm bias in the tropical Indian Ocean (Chowdary et al. 2016), and a salty and warm bias in the Antarctic Intermediate Water simulations is largely caused by the misrepresentation of oceanic mixing processes (Zhu et al. 2018).
Based on the previous studies, we will continue to investigate the subsurface temperature bias in the TNP, focusing mainly on the bias contributions from the atmospheric wind forcing fields and the oceanic vertical mixing parameterization by performing CMIP-based analyses and ocean model-based experiments. This paper is organized as follows. Section 2 describes the datasets from CMIP simulations, reanalysis and observational products for model evaluation, and model configurations for MOM5-based numerical experiments. Section 3 describes the characteristics of subsurface temperature biases in the TNP. The influences of the temperature biases on the upper-ocean thermal and current structures are discussed in Sect. 4. Atmospheric and oceanic origins of the subsurface temperature biases are investigated in Sects. 5 and 6, respectively. Finally, summaries and discussions are given in Sect. 7.

Datasets and ocean model used
This study is primarily based on the historical simulations from 53 CMIP6 models (Eyring et al. 2016), which are available online at https:// esgf-node. llnl. gov/ proje cts/ cmip6/. The historical simulations are forced by the observed greenhouse gases, solar forcing, and volcanic aerosols from 1850-2014, providing an opportunity to evaluate model abilities to simulate the past climate. In our analysis, the last 35 years  of historical simulations are selected for model evaluations. All the selected CMIP6 outputs are interpolated onto a 1° horizontal grid, and the potential temperature fields are further interpolated to 87 standard levels with a vertical resolution of 10 m near the sea surface. As we are interested in the long-term mean biases in climate simulations, the mean states for all fields are calculated by averaging the entire selected period. In addition, 41 CMIP5 models and 16 Ocean Model Intercomparison Project (OMIP) (Griffies et al. 2016) models are also used in our study. All the models are listed in Table 1.
In order to evaluate the subsurface temperature biases in CMIP6 simulations, the EN4 objective analyses of subsurface temperature from the Met Office Hadley Centre are used (Good et al. 2013). The EN4 product consists of temperature and salinity fields from 1900 to the present, with a 1° horizontal resolution and 42 vertical levels. In addition, wind stress fields are taken from the fifth generation of ECMWF atmospheric reanalyses (ERA5; Copernicus Climate Change Service 2017) and the Scatterometer Climatology of Ocean Winds (SCOW); the latter is estimated from the 122-month record of the QuikSCAT wind measurements (Risien and Chelton 2008).
Bias quantification is made by conducting the MOM5based ocean-only experiments (Griffies et al. 2009). This ocean model has a 1° horizontal resolution with the meridional resolution increased to 1/3 near the equator, and 50 vertical levels with 10 m resolution in the upper 220 m. In order to investigate the relationship between the subsurface temperature bias in the TNP and the poorly prescribed atmospheric forcing fields, two ocean-only simulations are conducted. In the control run (the LY09 run), the Normal Year atmospheric forcing fields from Large and Yeager (2009), which are the required forcing fields by OMIP (Eyring et al. 2016), is applied. In the sensitivity run (the ERA5 run), forcing field of wind is replaced by the climatological 6-hourly wind forcing based on 1980-2014 ERA5 Reanalysis. Furthermore, in order to investigate the relationship between the subsurface temperature bias in the TNP and the parameterized oceanic vertical mixing processes, several ocean-only experiments are conducted, in which mixing strength in the upper ocean and the interior ocean is reduced. Specifically, the upper boundary layer mixing is modified by reducing the wind stirring effect in a Kraus-Turner type mixing scheme (Chen et al. 1994), and the interior mixing is constrained by using the Argo-derived diffusivity (Zhu and Zhang 2018b), whose details are given in Sect. 6.

Subsurface temperature bias in the TNP
In the northeastern tropical Pacific (NETP), a subsurface warm bias (SWB) emerges in most of CMIP6 simulations, with its bias center located near (10° N, 130° W) at the depth of 100 m (Fig. 1a). The multimodel mean bias is about 4 °C with the maximum greater than 8 °C (Fig. 2). Although great efforts have been made to improve climate model performances over the past decade, the SWB in the NETP has not been reduced significantly from CMIP5 to CMIP6 (Fig. 1b). Figure 3 shows the seasonal variation of the SWB. The warm bias is deeper during April-June, and is relatively shallower during boreal winter, accompanied by a deep thermocline bias throughout the year.
Although this SWB occupies a small domain of the tropical Pacific basin, its influences might be widespread as the upper-ocean thermal structures in the TNP control the subtropical-tropical water exchanges. In the next section, some local and remote consequences of the SWB are further discussed.

The consequences of the SWB
In order to investigate the relationships between the SWB and the simulated oceanic circulation pathways in CMIP6 outputs, we define two model groups according to the magnitude of SWB (Fig. 4a): 14 models are seen to have large SWB (Group1), and 10 models are seen to have trivial or little SWB (Group2). As shown in Fig. 4b-d, the CMIP6 models with a large SWB in the NETP tend to produce much flatter isopycnals over the central TNP, which acts to reduce interior water transport. As a consequence, the equatorward water transport at 50-150 m in the westerncentral TNP is weaker than that in the models with a small SWB. Furthermore, the equatorward water transport in the subsurface TNP regulates the supply of cold water that upwells along the equator, which can be associated with the decadal changes in SST over the central and eastern equatorial Pacific (McPhaden and Zhang 2002;Capotondi et al. 2005). Figure 5a and b contrast the equatorial Pacific upwelling simulated in Group1 and Group2. Consistent with our understanding of subtropical-tropical water exchange, models with the large SWB tend to produce a weak equatorial upwelling (Fig. 5c). Thus, the cooling effect by the equatorial upwelling in Group1 is weak, and the upper-ocean temperature difference between Group1 and Group2 is positive ( Fig. 6a and b). Figure 6c shows the scatterplots of the SWB versus the equatorial Pacific upper-ocean temperature bias among the 53 CMIP6 models. Obviously, these two quantities show a positive correlation (R = 0.75), indicating that the CMIP6 models with the larger SWB in the NETP tend to produce a warmer equatorial upper-ocean.
This SWB is highly related to the simulations in the upper equatorial Pacific. Thus it is critical to understand the mechanisms that lead to this bias. This SWB is located in an important region where a positive wind stress curl (WSC) produces a positive Ekman pumping near the Intertropical Convergence Zone (ITCZ). As a consequence, a ridge-like climatological thermocline arises to form a potential vorticity barrier blocking the local water exchange between the subtropics and tropics GISS-E2-1-H (Lu and McCreary 1995;Johnson and McPhaden 1999). Therefore, WSC over the NETP might be underestimated by CMIP6 models, and the consequent Ekman upwelling is too weak to maintain the shallow thermocline as observed.
In the next section, the relationship between the SWB and the wind stress simulations is examined. Figure 7a shows the linear regression of the intermodel WSC onto the normalized SWB series (Fig. 2). The NETP happens to be the region where negative regression  1 Spatial distributions of multimodel annual mean temperature bias greater than 3 °C in CMIP6 and CMIP5 multimodel ensemble coefficients arise. As the weak WSC tends to produce a weak Ekman upwelling, intermodel differences in SWB are largely explained by the differences in the simulations of local WSC. However, whether the WSC intensity in the NETP is underestimated by CMIP6 models is quite uncertain. Figure 7b shows the annual-mean WSC difference between the CMIP6 multimodel ensemble (MME) and the ERA5. It seems plausible that the WSC intensity in the NETP is indeed underestimated by CMIP6 models when the wind fields in ERA5 are considered to be realistic. But it is widely accepted that wind measurements by Quik-SCAT are more reliable, and so the QuikSCAT winds are widely used to correct the atmospheric reanalyses winds (Large and Yeager 2009;Tsujino et al. 2018). Figure 7c  shows the WSC difference between the CMIP6 MME and the SCOW. It is obvious that the positive WSC difference extends farther south, almost occupying the entire 10° N. Hence, the WSC bias in the NETP is positive rather than negative when the QuikSCAT measurements are considered to be more realistic than ERA5. However, the reliability of the QuikSCAT measurements has been questioned. Many previous studies have found that the accuracy of the QuikSCAT measurements can be degraded by rain (Weissman et al. 2002;Draper and Long 2004;Sun et al. 2019). Therefore, the WSC in the ITCZ region can be poorly measured by the QuikSCAT, and correcting the reanalysis winds towards the QuikSCAT winds might introduce a large error. Nevertheless, the forcing fields of wind in OMIP experiments just come from the NCEP and JRA55 wind fields corrected by the QuikSCAT measurements (Griffies et al. 2016). Thus, it is not surprising that this SWB also exists in the OMIP simulations (Fig. 8).

The contribution of surface wind stress to SWB
To further investigate the relationship between the SWB and the wind forcing fields, two MOM5-based ocean-only experiments are conducted. In the LY09 run, the prescribed wind forcing fields are the same as those in the OMIP experiments. In the ERA5 run, ERA5 reanalysis winds are used to drive the ocean model. Both runs are integrated for 30 years using a Kraus-Turner-type vertical mixing scheme (Chen et al. 1994), and the outputs for the last 5-year are chosen to calculate the mean state. Figure 9 shows the temperature differences between the ERA5 run and the LY09 run. Compared with that in the LY09 run (Fig. 9a), the SWB is greatly reduced by ~ 6 °C in the ERA5 run (Fig. 9c), implying that the failure of reproducing the subsurface temperature distribution in the Tropical North Pacific is due to the poorly prescribed wind forcing fields in OMIP experiments.
Although the above analyses confirm a significant contribution of the wind stress simulations to the SWB in the NETP, its quantitative contribution is still unknown as it is difficult to obtain the accurate wind measurements and to realistically describe the WSC distributions in the tropics for ocean modeling. Here, we further discuss the possible mechanisms impacting the wind simulations in the NETP. As shown in Fig. 10, the SWB is closely related to the northeasterly wind in the NETP, which helps to produce a negative WSC on its right flank (Fig. 7a). Song and Zhang (2020) have found that the climate models fail to reproduce the seasonal wind reversal of the North American monsoon, leading to a year-round northeasterly wind bias in the NETP. Consistent with their study, Fig. 10 shows positive regression coefficients between the sea level pressure and the SWB over the North America. Thus, the overestimated sea level pressure over the North America acts to produce a year-round northeasterly wind bias in the NETP, which further causes a negative WSC bias on its right flank. This negative WSC bias acts to suppress the local Ekman upwelling, leading to the emergence of the SWB in the NETP as represented in the CMIP simulations.

The contribution of oceanic vertical mixing scheme to SWB
Apart from the erroneous simulations in atmospheric state, the deficiencies in ocean models can also contribute to the SWB in the NETP. The overly strong penetration of the ocean boundary layer mixing tends to warm the subsurface layer. Thus, sensitivity  (Fig. 2). The dots denote the region where P-value is smaller than 0.05. b, c WSC bias relative to ERA5 and SCOW, respectively experiments are conducted to investigate the relationship between the SWB in NETP and the overly strong mixing effect. In this subsection, a Kraus-Turner-type vertical mixing scheme (Niiler 1977;Chen et al. 1994) is applied to determine the mixed layer depth as follows: where Δb is the buoyancy jump across the base of ocean surface boundary layer, w e = h∕ t is the entrainment velocity, h is the depth of ocean surface boundary layer, Fig. 9 a Temperature bias at the depth of 100 m and b the vertical-season section of temperature bias horizontally averaged over the NETP in the LY09 run. c, d The differences between the ERA5 run and the LY09 run  (Acreman and Jeffery 2007;Zhu and Zhang 2018a). This prognostic equation of h has been employed into MOM5 to depict the evolution of ocean surface boundary layer. When the h is determined by the Eq. (1), vertical mixing coefficients for the model layers within the h are assigned to be a constant value of 5 × 10 -3 m 2 s −1 , the maximum observed diffusivity reported by Peters et al. (1988). Thus, the overly deep penetration of boundary layer mixing demonstrated in Fig. 12 can be relieved by reducing the m 0 in the Eq. (1). Two MOM5 based ocean-only simulations are thus conducted: m 0 = 0.4 is taken in the control run; in the reduced m 0 run, the m 0 is prescribed as in which "lon" and "lat" are the longitude and latitude of a model grid point. Below the h, shear instability mixing for the ocean interior is parameterized as a function of Richardson number (Peters et al. 1988), and background diffusivity is prescribed with a constant value (10 -5 m 2 s −1 ). Both runs are integrated for 30 years using the atmospheric climatological forcing fields from LY09, and the outputs for the last 5 years are selected for comparisons. Figure 13a shows the annual-mean difference in ocean surface boundary layer depth between the reduced m 0 run and the control run. In accord with Eq. (2), surface boundary layer shoals by ~ 15 m in the NETP when the m 0 is reduced, greatly removing the too deep surface boundary layer depth bias in Fig. 12b. As a consequence, subsurface layer cools by ~ 1 °C in the NETP (Fig. 13b). Moreover, the SWB is found to be more pronounced during the first half of a year (Fig. 9b). By reducing the penetration depth of the boundary layer mixing, improvements in the subsurface temperature simulations are much substantial in the reduced m 0 run during the first half of a year (Fig. 13c). Besides, surface boundary layer in the reduced m 0 run is shallower than that in the control run, and their difference is large during the boreal winter when the subsurface cooling effect is the most pronounced, revealing that the overly strong penetration of the mixing effect and the ocean surface boundary layer deepening indeed contribute to the SWB in the NETP.
Besides the misrepresented boundary layer mixing, the overestimated vertical mixing in the ocean interior may also contribute to the SWB by inhibiting the downward heat transport to the deeper ocean. Indeed, our previous study (Zhu and Zhang 2018b) found that the overestimated vertical diffusivity along the equatorial Pacific contributes to a Fig. 12 a Microstructure profiles from the LADDER project (kindly provided by Professor Andreas M. Thurnherr at Lamont-Doherty Earth Observatory). These data were collected in the eastern tropical Pacific near the crest of the East Pacific Rise (9°30′-10° N, 103°45′-105° W) using Vertical Microstructure Profiler (VMP) during November to December, 2007. b The corresponding vertical diffusivity in MOM5 based ocean-only simulation warm bias within the equatorial thermocline, and employing the Argo-derived background diffusivity produces a reduction in subsurface warm bias and a more realistic equatorial thermocline. Off the equator, the spatial pattern of diapycnal mixing in the tropical Pacific inferred from the strain-based finescale parameterization (Kunze et al. 2006) also confirms the weak mixing strength in the NETP, coinciding with the location of the SWB (Fig. 14a). Thus similar to our previous studies (Zhu and Zhang 2018b;Zhu et al. 2020Zhu et al. , 2021, diapycnal diffusivity over the NETP are calculated using the strain-based finescale parameterization: where ⟨ 2 z ⟩ is the observed strain variance which can be calculated from the Argo profiles (available online at ftp:// ftp. ifrem er. fr/ ifrem er/ argo/) with 2-10 m vertical resolution, is the strain variance from the Garrett-Munk model spectrum, h R is a function of the ratio between shear and strain variance, and j(f ∕N) is a latitudinal correction given the influence of the Coriolis force on the internal wave breaking. Diffusivity estimates derived from Argo profiles are grouped into 2° × 2° bins, and the median diffusivity averaged between 250-500 m is used to present the diapycnal diffusivity in each bin. More details can be found in Zhu and Zhang (2018b). Two ocean-only experiments are conducted to investigate the relationship between the SWB and diapycnal mixing intensity in the NETP. In the control run, background diffusivity is taken as the commonly used value (10 -5 m 2 s −1 ), a value that is considered to be too large compared to what is observed. In the reduced background diffusivity (RBD) run considering the observed weak diapycnal mixing, background diffusivity over the NETP is replaced by the Argo-derived one (Fig. 14a). It is worth noting that the finescale parameterization has a relatively high uncertainty. This is the reason why the VMP observations from the LADDER project are also used in this study (Fig. 12a). In fact, the diffusivities derived from both the LADDER observations and the finescale parameterization are ~ O(10 -6 m 2 s −1 ) in the interior ocean, giving confidence that the Argo-derived diffusivity is able to realistically represent the vertical mixing strength in the NETP. Two runs are integrated for 30 years using the LY09 forcing fields, and the vertical mixing scheme for the surface boundary layer is also the Kraus-Turner-type vertical mixing scheme. The temperature and diffusivity differences averaged over the last 5 years are shown in Fig. 14b-d. It is shown that the total diffusivity is reduced by nearly an order of magnitude below the surface boundary layer when the Argo-derived diffusivity is applied (Fig. 14c). The reduced diffusivity inhibits the heat transfer downward, so that the SWB can be reduced when constraining the background diffusivity to match the observations. The above discussions support the fact that the poorly represented vertical diffusivity in the climate models contributes to the SWB in the NETP. But it is still unclear whether the change in vertical diffusivity magnitude can fully explain the subsurface cooling in the reduced m 0 run and the RBD run. Particularly, Jia et al. (2015) find that the ocean temperature response to the change in the vertical structure of the vertical diffusivity is often stronger than the response to the change in the magnitude of the vertical diffusivity. Thus, following Furue et al. (2015), Jia et al. (2015) and Jia et al. (2021), the 1-D equation of temperature change resulting from a change in the vertical diffusivity is performed to quantify the individual contributions to the SWB reduction: where ΔT and Δk v are the changes in temperature and vertical diffusivity when the reduced m 0 and the Argo-derived background diffusivity are applied. The first term on the right-hand side represents the temperature response to a change in the vertical diffusivity, whereas the second term is associated with the vertical gradient of the change in vertical diffusivity. Figure 15a and b show these two terms in the reduced m 0 run. Considering the subsurface cooling mainly appears between 100-150 m in the reduced m 0 run (Fig. 13c), the warm bias reduction is primarily caused by the change in the vertical structure of the vertical diffusivity. But in the RBD run, subsurface cooling appears between 100-125 m (Fig. 14d), where the first term appears to have a cooling effect (Fig. 15c and d). Thus, the warm bias reduction in the RBD run is caused by the weakened magnitude of the vertical diffusivity when the Argo-derived diffusivity is applied.
Our modeling experiments clearly demonstrate that the uncertainties in the oceanic vertical mixing parameterization can produce a ~ 3 °C warm bias, and the poorly simulated wind stress fields (Fig. 9c) can cause a ~ 6 °C warm bias. Putting together, these two effects could account for the magnitude of the SWB in coupled simulations (Fig. 16).

Summary and discussion
Thermal structure and variability in the subsurface oceans are known to play an important role in the climate change, and the studies of oceanic subsurface temperature biases are important for understanding the simulations and predictions of global energy and heat redistribution. Though being substantial in climate simulations, subsurface biases have not received much attention yet. In this study, subsurface temperature biases in the TNP are investigated using the newly released CMIP6 historical simulations. It is found that almost all the CMIP6 simulations have a pronounced SWB that is persistent throughout the year in the NETP, and the model updates from CMIP5 to CMIP6 do not provide obvious bias alleviations for temperature biases. The SWB is primarily caused by the model deficiencies in simulating the surface WSC. For a similar reason, the poorly prescribed wind forcing also causes a pronounced SWB in the NETP in OMIP simulations. Besides, this SWB can also be partly attributed to the uncertainties in the representations of the oceanic vertical mixing processes. Compared with the observations, the estimated vertical diffusivity in ocean models is misrepresented in both the upper boundary layer and the ocean interior. By constraining the diffusivity to match observations, the SWB in the NETP is reduced by ~ 3 °C. The consequence of the SWB in the NETP is further examined using the CMIP6 products. The SWB can influence the simulations of oceanic circulations and equatorial upper-ocean thermal structure. The CMIP6 simulations with a large SWB tend to produce a weak equatorward water transport in the subsurface TNP, a weak equatorial upwelling and thus a warm upper ocean layer along the equator.
Substantial subsurface temperature biases are also seen in other regions of the world ocean (Fig. 17a). Unfortunately, it seems that no substantial improvements are achieved from CMIP5 to CMIP6 (Fig. 17b), and much more efforts are needed to understand the sources of subsurface biases and ultimately to remove them. It is interesting to note that the subsurface temperature bias is often accompanied by a bias in salinity simulation, and the biases in temperature and salinity are compensated for by each other in terms of their effects on density. For example, the large subsurface cold bias centered at the Pacific sector (8° S, 200 m) is accompanied by a fresh bias. Because these two biases tend to be compensated for by each other in terms of their effects on density, density bias is reduced. Compensation of temperature and salinity is a property commonly observed in water masses (Lilly et al. 1999;Rudnick and Ferrari 1999). Thus, subsurface thermohaline biases might be closely related to the water mass formation and transformation in the global oceans, which needs to be analyzed in the future. Figure 6 shows a positive correlation between the SWB and the equatorial upper-ocean temperature in CMIP6 models. The positive correlation may be not related to the vertical diffusivity parameterized in the NETP, because the background diffusivity in the NETP is ~ 10 -5 m 2 s −1 among the CMIP6 models and the SWB produced by the errors in the vertical diffusivity should not be so scattered. It is interesting to note that the intermodel correlation between The basinwide zonally averaged biases for temperature (colors) and salinity (contours) biases in a CMIP6 and b CMIP5 the precipitation and the SWB index (Fig. 18) has a similar spatial pattern to the second EOF mode of the tropical precipitation bias (Fig. 2b in Li and Xie 2014). Thus, the positive correlation reflects the links among the double-ITCZ bias, the SWB and the equatorial upper-ocean temperature. Further quantitative assessment is needed.
Climate models suffer from a significant cold bias in the Pacific cold tongue region (Fig. 6c). Thus, it seems likely that increasing the SWB helps to reduce the too cold tongue bias in the equatorial Pacific, and the mechanisms for reducing one bias may have the effect of increasing another. However, the origins of Pacific cold tongue bias are still unclear. Atmospheric model deficiencies can make a great contribution because the ocean models participating the OMIP tend to produce a warm bias in the Pacific cold tongue region (Tsujino et al. 2020). If the atmospheric model produces a background cold bias in the cold tongue region, the warm SST bias produced by the SWB can be erased. Thus further modeling studies are needed to understand the origins of Pacific cold tongue bias.
Previous studies have found that the cold bias in the upper tropical Pacific Ocean affects the ENSO simulations in climate models, leading to an intermodel diversity of ENSO representations in CMIP5 simulations (Vannière et al. 2013;Kim et al. 2014;Zhang et al. 2020). The upper ocean temperature bias in the equatorial Pacific is a possible cause of forecast errors in ENSO amplitude (Kim et al. 2017). Thus, the SWB might be linked to the ENSO simulations in climate models, and further research efforts are needed to test this hypothesis.