Remote effect of model systematic bias in tropical SST on the cold bias over the Tibetan Plateau

Most state-of-the-art climate models substantially underestimate the near-surface air temperature (SAT) over the Tibetan Plateau (TP), especially for the cold season. While previous studies have attributed this cold bias to local factors such as the elevation difference, coarse resolution, and excessive snow cover, this investigation addresses the potential contributions of the systematic bias of tropical sea surface temperature (SST) to the TP cold bias. Experiments with the NCAR Community Atmosphere Model demonstrate that tropical SST bias results in an apparent cold bias over the TP, especially in boreal winter, and explains about 40% of the deviation in multi-model mean SAT over the TP relative to station observations. Forced by the tropical SST bias, heat flux exhibits an anomalous divergence over the plateau, causing a cooling center in the mid- and lower-troposphere over the TP. This atmospheric cooling in turn leads to a reduction of the downward longwave radiative fluxes reaching the surface, less energy supply, and thus a cold bias over the TP.


Introduction
The Tibetan Plateau (TP) exerts a profound impact on atmospheric circulation due to its extraordinary altitude and horizontal extent, and thereby influences the global climate, especially the Asian monsoon climate (Flohn 1957;Yanai et al. 1992;Duan and Wu 2005;Lu et al. 2018Lu et al. , 2019. Its elevated heating, which is closely associated with the difference between surface skin temperature and surface air temperature (SAT) and related to surface wind speed, is a dominant impacting factor during most of the year (Zhao et al. 2018). As an indicator of the thermal condition The original online version of this article was revised: "The figures have been placed correctly and affiliation details of author has been updated in the original article." Xiaoming Hu huxm6@mail.sysu.edu.cn Song Yang yangsong3@mail.sysu.edu.cn 1 (GCMs) from the Coupled Model Intercomparison Project Phase 5 (CMIP5) exhibit a substantially cold bias over the plateau, especially during the cold season (Su et al. 2013;Duan et al. 2014;Chen and Frauenfeld 2014). The mean cold bias in CMIP6 models is only slightly reduced in the summertime compared with that in CMIP5 models (Xie & Wang 2021;Lun et al. 2021). Massive effort has been devoted to understanding the sources of the TP cold bias. It was found that the bias of interpolated temperature was in proportion to the increase in local elevation and topographical complexity, and that the improvement of interpolated surface temperature could be achieved via "topographic correction" (Zhao et al. 2008;Hu et al. 2014;Wang et al. 2015). Moreover, the unrealistic local parameterization schemes in models influence surface heat transfer and resistance, which is partly responsible for the TP SAT bias (Yang et al. 2007;Zhuo et al. 2016). Lee et al. (2019) also reported that the lack of ice radiative effect (i.e. the miss of emission from falling snow) in most CMIP5 models, causing more solar radiation to penetrate through the optically thin atmosphere and less downward longwave (LW) radiation, could also explain a portion of the model biases of surface radiative fluxes and surface temperature over the TP. Several studies focused on the land processes in model simulations, especially the snow processes, further found that overestimation of TP snow could lead to larger solar energy reflection and then a colder surface and that the simulated cold bias would be improved with more realistic albedo over the TP (Rangwala et al. 2010;Ghatak et al. 2014;You et al. 2016;Meng et al. 2018).
Apparently, the aforementioned studies have been mainly focused on the local processes that may induce cold bias in models. However, the TP SAT is simultaneously regulated by tropical sea surface temperature (SST). The spring SST and air-sea coupling of the Indian Ocean significantly affect the TP heat source in summer, and the Indian Ocean SST and TP SAT exhibit a possibly teleconnected mode (Zhang et al. 2006;Ji et al. 2018;Zhao et al. 2018;He et al. 2019;Wang et al. 2019). Besides, El Niño is closely associated with an increase in wintertime snowfall over the plateau, causing a colder TP (Shaman and Tziperman 2005); Wang and Xu 2018) . The remote SST in the warm pool could also govern TP warming by modulating the coherent variations of tropical easterly jet stream and subtropical westerly jet stream (Wang et al. 2012). Ma et al. (2017) suggested that the tropical and extratropical oceanic internal variability exerted a relatively uniform warming effect on the whole TP during 1980-2012. While SST especially the tropical SST shows a profound influence on the TP SAT as discussed above, the SST itself displays a remarkable systematic bias in the state-of-the-art models. SST errors are comparable or even larger than the observed interannual variability and the projected change of SST in the 21st century (Li and Xie 2012). The equatorial cold tongue penetrates too far westward in the Pacific with a double intertropical convergence zone (Mechoso et al. 1995;Li and Xie 2014;Zhou et al. 2020), and there exist warm SST biases over the southeastern Pacific off South America (Xie et al. 2007;de Szoeke et al. 2010). Also, the warm SST bias in the southeastern tropical Atlantic is a common problem in many previous and current climate models (Xu et al. 2014;Zhang et al. 2014;Exarchou et al. 2018;Kurian et al. 2021), and coupled atmosphere-ocean GCMs show substantial difficulties in correctly simulating the sign of the nearequatorial zonal gradient of the Atlantic Ocean SST (Davey. et al. 2002;DeWitt 2005;Song et al. 2012).
Given that tropical SST influences the TP SAT extensively and that climate models often suffer from large biases especially over the tropical Pacific and Atlantic oceans, the unrealistic climatological SST in model simulations may also contribute to the TP SAT bias. The current study is aimed to explore the possible role of tropical SST biases in the modeled TP cold biases. In Sect. 2, datasets, model design, and analysis methods are introduced. Section 3 describes the TP cold bias and the tropical SST bias in climate models. Section 4 presents the influence of tropical SST bias forcing on the TP SAT bias and the responsible mechanisms. Summary and discussion are presented in Sect. 5.

Data
The Extended Reconstructed SST version 5 (ERSSTv5; Huang et al. 2017) dataset from the National Oceanic and Atmospheric Administration (NOAA) and the SAT dataset from the National Meteorological Information Center (China Meteorological Administration, http://data.cma.cn/ data/cdcdetail/dataCode/SURF_CLI_CHN_TEM_MON_ GRID_0.5.html) are used as the observed references to calculate models' systematic biases. The monthly gridded SAT data during 1961-2018 with a high spatial resolution of 0.5˚ × 0.5˚ are interpolated from the observation records of more than 2400 Chinese meteorological stations. The SST segment from the ERSSTv5 dataset is from 1950 to 1999 on 2˚ × 2˚ resolution grids. The domain of the TP is limited to the area higher than 1500 m over 25˚-40˚N/75˚-105˚E.
Model biases are defined as the deviations of model simulations in historical experiments from the observed monthly SAT and surface skin temperature. Note that the values of surface skin temperature over the oceans are the same as the SST in climate models. The systematic bias is the mean of model deviations obtained from 28 CMIP5 models and 20 CMIP6 models. Tables S1 and S2 in the supplemental materials provide a summary of the details of all model simulations. To facilitate a comparison with observations, all the simulated SAT and SST data by the GCMs are interpolated to 0.5˚ × 0.5˚ and 2˚ × 2˚ grids, respectively. Since previous studies have demonstrated a robust annual cycle of the strength of TP cold bias (Su et al. 2013;Jia et al. 2019), all model biases analyzed in this study are measured by monthly mean values.

Model and experiment design
The model used in this study is the Community Atmosphere Model version 4 (CAM4), the atmospheric component of the Community Earth System Model version 1.2.2 (CESM 1.2.2) from the National Center for Atmospheric Research (NCAR, Neale et al. 2013), with a horizontal resolution of 1.25˚ longitude by 0.9˚ latitude and 26 vertical levels. The CESM as a fully-coupled global climate model can provide state-of-the-art simulations of the climate states and be used for a broad class of scientific research (Hurrell et al. 2013).
This study is mainly based on six model experiments, whose design is briefly described in Table 1. The simulations are implemented by F_AMIP component setting, in which both the atmosphere model and land surface model are active. The control experiment, referred to as CTRL, is forced by the observed monthly-mean global SST from 1979 to 2005. The surface boundary forcing dataset for uncoupled simulations with CAM4 is a merged product based on the monthly-mean Hadley Center Sea Ice and Sea Surface Temperature dataset (HadISST) and the NOAA Optimum Interpolation Sea Surface Temperature dataset (OISST). The details about the default optional product have been described in Hurrell et al. 2008). It should be noted that the TP SAT bias exists in the CTRL experiment which is forced by observed SST. The TP bias can be attributed to not only remote factors (e.g. tropical SST bias) but also local factors such as surface heat fluxes, snow cover, and so on. Therefore, the TP SAT bias can still exist even though the model is forced by observed SST. We can remove this inherent bias by using the difference between the perturbation and control simulations to represent the contribution of tropical SST bias. The second to fifth sensitivity experiments are similar to the CTRL, except that the annual cycle of tropical (30˚S-30˚N) SST bias is imposed on the observed monthlymean global SSTs over all the three basins (referred to as SEN-ALL), the Pacific Ocean basin only (referred to as SEN-PAC), the Atlantic Ocean only (referred to as SEN-ATL), and the Indian Ocean only (referred to as SEN-IND), respectively. It is noted that the tropical SST bias is given by the climatological monthly SST difference between the CMIP5 multi-model ensemble (MME) and the ERSSTv5 during 1950-1999. The last sensitivity experiment (referred to as SEN-ALL6) is similar to SEN-ALL, except that the tropical SST bias is replaced by the difference between the CMIP6 MME and the ERSSTv5. All experiments involve model runs for 27 years and the results from the last 25 years are analyzed.

Downward LW radiation decomposition
We apply the radiative transfer model by Liou (1992, 1993) to decompose the total downward LW radiative surface energy flux perturbations (∆R s ↓ ) between CTRL and SEN-ALL into the sum of partial downward LW surface energy flux perturbations due to the changes in CO 2 (∆R CO2 s ↓ ), ozone (∆R O3 s ↓ ), cloud (∆R Cloud s ↓ ), water vapor (∆R W V s ↓ ), and air temperature (∆R T A s ↓ ). The decomposition can be expressed mathematically as: where the subscript s denotes the near-surface layer and R represents LW radiation. In fact, each term in the right-hand side of Eq.
(1) is calculated as the difference between two runs of the radiative transfer model: one with the fields of CO 2 , ozone, clouds, water vapor, and air temperature derived from CTRL and the other with identical inputs except the fields denoted by the superscripts, which are taken from experiment SEN-ALL. For the left-hand side ∆R s ↓ , the inputs of the first radiative transfer model calculation are all taken from CTRL and the fields for the second calculation are all from SEN-ALL. Considering the same external

SEN-IND
Observed SST forcing + CMIP5 Indian Ocean basin tropical (30°S-30°N) SST bias with annual cycle only CMIP6 all-basin tropical SST bias run

SEN-ALL6
Observed SST forcing + CMIP6 tropical (30°S-30°N) SST bias with annual cycle only cloud ∆T Cloud , surface albedo (∆T α ), surface dynamic processes ∆T NRS , and atmospheric dynamic processes ∆T NRA . Then, the DLWRF anomalies at the surface associated with the partial air temperature anomalies related to ∆T X can be obtained according to (after the linearization approximation), where − ∂Rs ∂T j corresponds to the DLWRF perturbations at the surface layer due to a 1-K warming at the jth layer (for model outputs, J = 26). By substituting the atmospheric component of all seven ∆T X terms obtained by Eq. (2) Sejas and Cai (2016) and Hu et al. (2017), they represent the effect of air temperature response to the energy flux perturbations at the surface whereas other three ∆R T A_X s ↓ terms include the effect of air temperature response to the energy flux perturbations both at the surface and in the atmosphere. In summary, we obtain, to further decompose the DLWRF anomalies induced by air temperatures into the partial DLWRF anomalies due to cloud, water vapor, and atmospheric dynamic processes. .

Offline error correction
According to Eq. (1), offline errors come from the difference between ∆R s ↓ and the sum of ∆R X s ↓ . They are attributed to (i) using the climatological mean fields as inputs of the radiative transfer model in our offline calculations instead of the climatological mean of taking the radiative heating rates calculated by the instantaneous fields, especially for the cloud fields (Kang et al. 2020), (ii) adopting different radiative transfer models, and (iii) excluding the changes in other processes such as aerosols.
Following Hu et al. (2017), we identify the dominant terms that are responsible for ∆R err s ↓ by comparing the spatial distribution of each term ∆R X s ↓ in the right-hand side of Eq. (1) with error (∆R err s ↓). We evaluate the pattern correlation ( r X ) between each individual process ∆R X s ↓ and error ∆R err s ↓ over region A, namely, Furthermore, we apply the coupled feedback response analysis method (CFRAM) by Cai and Lu (2009) and Lu and Cai (2009) to decompose the total air temperature changes into partial air temperature changes due to individual processes including water vapor and others. This approach allows us to further evaluate the downward LW radiative flux (DLWRF) anomalies at the surface associated with the partial air temperature anomalies due to individual processes. To be specific, we use to evaluate the vertical profile of partial temperature anomalies ( ∆T X in units of K) from the vertical profile of partial radiative or non-radiative heating anomalies due to a specific process "X" (∆Q X in units of W/m 2 ). In Eq.
(2), (∂R/∂T ) is the Planck feedback matrix whose jth column corresponds to the vertical profile of the radiative energy perturbation due to 1-K warming at the jth layer. There are five radiative processes of "X", namely, CO 2 concentration ∆Q CO2 , ozone ∆Q O3 , water vapor ∆Q W V , cloud ∆Q Cloud , and surface albedo (∆Q α ). All the five terms can be obtained from the same calculation as for the right-side terms in Eq. (1). Besides, non-radiative processes ∆Q NR such as surface turbulent sensible and latent heat fluxes, atmospheric large-scale advective and convective processes, and heat storage terms in the surface and atmospheric layers can also induce energy flux convergence anomalies. Estimation of the vertical profile of total nonradiative heating perturbations is based on the energy balance between radiative and non-radiative processes. That is, where ∆Q R corresponds to the vertical profile of radiative energy flux convergence due to the changes in all the radiative processes. ∆Q R can be derived from the same-way calculation of the left-side term in Eq. (1). We further divide the ∆Q NR into the surface part ∆Q NRS and the atmospheric part ∆Q NRA , where ∆Q NRS represents the surface turbulent heat flux process and heat storage term, with its atmospheric component set to zero, and ∆Q NRA represents the atmospheric large-scale advective and convective processes, with its surface component set to zero. Finally, according to Eq.
(2), we obtain the total seven partial temperature anomaly fields induced by CO 2 concen- where the offline error mainly comes from the atmospheric dynamic process instead of the cloud field. Figure 1 shows the climatological monthly mean SAT over the TP in observation, CMIP5 models, and CMIP6 models. Both the CMIP5 and CMIP6 MMEs can reproduce the annual cycle of the TP SAT, but with a cold bias, especially in the cold season. Compared with the observation (black solid line), almost all CMIP5 models (dashed lines) significantly underestimate the TP SAT throughout the year, with the annual mean cold biases ranging from − 5.54 K to 0.04 K. The CMIP5 MME (purple solid line) presents a larger cold bias (-3.62 K) in boreal winter and a smaller one (-1.38 K) in boreal summer. Compared with CMIP5, the CMIP6 MME (green solid line) shows a progress in boreal summer, with a bias of about − 0.33 K, but a comparable cold bias in boreal winter, with a magnitude of -3.43 K. Figure 2a shows the correlation between the annualmean TP SAT and global SST in observations. The TP SAT is closely related to tropical SST, especially that in the tropical Pacific and Atlantic oceans. The multi-model mean of correlation coefficients derived from the individual CMIP5 models displays a remarkably high correlation coefficient

TP cold bias and tropical SST bias
where φ and λ are latitude and longitude, respectively. A indicates the TP region in this study, and a is the mean radius of the earth. It is found that r Cloud , the spatial pattern correlation between ∆R err s ↓ and ∆R Cloud s ↓ is the largest among the three processes, which is up to -0.5. The strong negative correlations suggest that the input of time-mean cloud fields tends to overestimate the DLWRF of the cloud feedback.
To address the above issue about overestimation, we correct the DLWRF perturbations of the offline error fields by removing the counterparts of the partial DLWRF due to the changes in clouds, namely, The corrected DLWRF perturbations due to the changes in cloud are shown in this paper while the uncorrected results are omitted.. The same method is used in the further decomposition of the downward LW flux perturbations due to the changes in atmospheric temperatures as given by Eq. (5), where (SST MME − SST OBS ) is the climatological annualmean SST bias derived from the CMIP5 MME (i.e. Figure 3m), SST j is the climatological annual mean SST of the CMIP5 model j, and SST OBS corresponds to the climatological annual mean SST in observations. Besides, other parameters have the same meanings as those in Eq. 6, except that A represents the tropical ocean domain (30˚S-30˚N). Figure 4 is the scatter plot of the two indices. The index E SST ranges from 0.61 to 1.22 K and E TP ranges from − 5.54 K to 0.04 K. The non-negative values of E SST indicate that all CMIP5 models have similar spatial patterns of SST bias with the same polarity as the MME SST bias shown in Fig. 3m. By definition, the models with smaller (larger) values of index E SST would perform better (worse) in reproducing the observed tropical SST or suffer from a smaller (greater) SST bias. The negative correlation (about − 0.4) between E TP and E SST of individual models indicates that smaller E SST corresponds to a larger E TP (i.e. a smaller cold bias), which indeed confirms that the models with a better (worse) ability of reproducing tropical SST would perform better (worse) in reproducing the SAT over the TP. Such a negative correlation is consistent with the observed feature shown in Fig. 2, namely cold SST bias over the tropical Pacific (which corresponds to positive projection onto Fig. 3m) would partially account for the cold SAT anomalies over the TP. Moreover, the non-perfect correlation between these two indices may be due to that the cold TP response to the tropical SST bias is nonlinear and that other factors (such as local effects) besides tropical SST bias also contribute to the cold TP bias.
To further substantiate the above statistical evidence that the systematic cold bias over the TP could be partially originated from the tropical SST bias, we perform a series of model experiments forced by tropical SST biases using the NCAR CAM4 model. The modeling experiments would help shed light on the underlying mechanisms for the cold response over the TP to tropical SST biases. Figure 5a -l show the spatial distributions of SAT difference between experiments SEN-ALL and CTRL from January to December. Under the CMIP5 MME tropical SST bias forcing, the TP SAT experiences an obvious cooling with strong signals in the cold season and relatively weaker signals in the warm season, which is similar to the feature shown in Fig. 1. The TP SAT anomaly forced by tropical SST bias can also reproduce the spatial distribution of MME cold bias in CMIP5/6 to some degree, which is larger over between the SAT and the SST (Fig. 2b), and the large values in tropical oceans indicate that the models reach a consensus on the significant correlation. The difference between Fig. 2a and b mainly rises from two aspects: the lack of observational data over the Southern Ocean produces potentially unrealistic correlation (Sumner et al. 2003;Rapizo et al. 2015) and climate models vary widely in their simulations of the Southern Ocean and its climate impact (Russell et al. 2018). Figure 3 displays the spatial pattern of monthly-mean and annual-mean tropical SST biases in CMIP5 MME relative to observations. The biases exhibit a small annual cycle, with a large region of cold bias and a smaller region of warm bias. The main cold bias appears in the tropical Pacific and Atlantic oceans, while the main warm bias emerges from the east coast of the mainland like Peru and western Africa. The tropical Indian Ocean shows a relatively large cold bias in the cold season and a smaller bias in the warm season. The tropical SST bias of the CMIP6 MME illustrates similar patterns with the CMIP5 MME, but with a relatively smaller magnitude (Fig. S1 in the supplemental materials).

Change in TP SAT forced by tropical SST bias
Given the close relationship between the TP SAT and the tropical SST, one may naturally ask whether the systematic cold bias over the TP presented in Fig. 1 can be partially explained by the tropical SST bias shown in Fig. 3. To gain an insight onto this question, we construct two indices: E SST and E TP , with the former for the severity of SST bias and the latter for the severity of TP SAT bias of individual models. Specifically, the index E TP is defined as the areal-mean and annual-mean TP SAT bias of the individual CMIP5 models, whereas the index E SST is the projection coefficients of tropical SST bias from individual models onto the MME SST bias. That is Apparently, the tropical Pacific shows the largest contribution to the total cold bias forced by the tropical SST biases. Nevertheless, it should be noted that the contribution of each basin is not linearly addable due to the complex interaction among the three ocean basins. Thus, the systematic tropical SST bias forcing of each ocean basin, especially the forcing of tropical Pacific bias, could explain a portion of the systematic cold bias over the TP. Here, we focus on the collective contributions of SST bias in all the three tropical basins to the TP cold bias in boreal winter, and the following analysis will be based on the western TP and smaller over the central and eastern TP. Figure 5m presents the area-averaged change in TP SAT in different tropical SST bias runs with respect to the CTRL, which displays that the tropical SST bias forcing by both CMIP5 MME and CMIP6 MME would lead to a cooler surface over the TP with an analogous annual cycle. In boreal winter, SEN-ALL (SEN-ALL6) shows a cold bias of about 1.49 K (1.31 K), which explains about 41% (38%) of the total winter cold bias in Fig. 1. Moreover, under the individual forcing of the tropical Pacific, Atlantic, and Indian Ocean SST biases, the TP exhibits an anomalous surface cooling of about − 0.65 K, -0.64 K, and − 0.31 K in winter SEN-ALL. The decreased upward LW radiation is balanced by (1) the obviously decreased downward LW radiation (Fig. 6b), (2) the increased downward and upward SW radiations ( Fig. 6d and e), and (3) the small changes in SH, LH, and G0 (Fig. 6g, h, and i). The changes in surface SH, surface LH, and G0 are relatively small compared with those in surface radiative heat fluxes and thus are neglectable. Also, the sum of these three terms seems positive, indicating that the net effect of the biases in SH, LH, and G0 would lead to more energy going to the surface layer (or less energy leaving the surface layer). Moreover, the total effect of LW radiation change exhibits a decreased net LW radiation (Fig. 6c, positive value indicating downward) over the nearly entire TP, suggesting that the decreased downward LW is dominant so that the surface loses energy via LW radiation. However, the total effect of SW radiation change, i.e. the net SW radiation (Fig. 6f, positive value indicating downward), shows a small increase over the central and western TP and a small decrease over the eastern TP, which indicates more solar radiation supply. These results of surface energy budget over the TP suggest that the plateau gains more energy via SW radiation and surface non-radiative processes but loses more energy via LW radiation. The decrease in surface downward LW radiation from the atmosphere is the dominant contributor to the decreased SAT over the TP.
Then, which factors are responsible for the decreased downward LW radiation under the forcing of tropical SST bias? To answer this question, we adopt the Liou (1992, 1993) radiative transfer model (see Sect. 2) to decompose the reduction of DLWRF.

Decomposition of winter downward LW radiative fluxes
The decreased surface DLWRF in Fig. 6b has been decomposed into the LW perturbations due to the changes in cloud (Fig. 7a), water vapor (Fig. 7b), and air temperature (Fig. 7c). The sum of the partial surface DLWRF (Fig. 7a, b, and c) is shown as the shading in Fig. 7d, which is close to the model output (the contour of Fig. 7d). The error term, namely the difference between the shading and the contour of Fig. 7d, is small, indicating that the decomposition is reasonable and valid (figure not shown). It is also found that the change in clouds mainly leads to a decreased DLWRF over the central TP, while the change in water vapor causes a generally uniformly-reduced DLWRF over the entire TP. Moreover, the change in air temperature is the largest contributor to the total downward LW perturbation, causing larger DLWRF over the western TP and smaller DLWRF over the eastern TP. The decreased DLWRF due to the change in air temperature is mainly caused by the cold atmosphere over the TP (Fig. S2 in the supplemental materials). the outputs from CTRL and SEN-ALL for the December-January-February (DJF) means.

Surface energy budget analysis of winter TP cooling
The sources of winter surface cooling over the TP can be revealed via a surface energy budget analysis. According to the Stefan-Boltzmann law, the upward LW radiation is directly proportional to surface temperature. The surface cooling over the TP is equivalent to the decrease in upward LW radiation. Based on the surface energy budget equation, this decrease in upward LW radiation is coupled with the changes in downward LW radiation, upward and downward shortwave (SW) radiation, latent heat flux (LH), sensible heat flux (SH), and ground heat flux (G0) at the surface. That is, where the subscript s denotes the surface layer, and R and S represent LW and SW radiations, respectively. Figure 6 illustrates the spatial distributions of the changes in winter surface energy fluxes between SEN-ALL and CTRL. The upward LW radiation (Fig. 6a) is much smaller in SEN-ALL than in CTRL, which is consistent with the cooler TP in Since the decreased DLWRF is mainly affected by atmospheric dynamic processes, it is worth exploring the possible link between the change in atmospheric circulation and the decreased DLWRF induced by anomalous atmospheric cooling.

Changes in winter heat flux convergence and large-scale atmospheric circulation
As discussed above, atmospheric dynamic processes play an important role in the anomalous atmospheric cooling that leads to a decreased DLWRF over the TP. Heat flux divergence is a main atmospheric dynamic process that may be responsible for the change in air temperature. Figure 9a, which is the same as Fig. 8c but with a different color bar for a better comparison in the following analysis, shows the DLWRF perturbations due to the anomalous atmospheric temperature induced by atmospheric dynamic processes, To reveal the origins of the anomalous atmospheric cooling over the TP, we further decompose the surface DLWRF perturbations induced by air temperature anomalies into three terms, i.e. partial air temperature anomalies associated with the changes in cloud, water vapor, and large-scale atmospheric heat transport. As shown in Fig. 8a, the reduction of DLWRF due to the cloud-induced air temperature change is mainly located over the southeastern TP and the western TP, and the reduction due to the water vapor induced change in air temperature is approximately uniform over the TP with a relatively larger magnitude over the western TP (Fig. 8b). Moreover, the reduction due to air temperature anomalies is dominated by the non-radiative heating contribution from the atmospheric advective and convective processes over the central TP (Fig. 8c). The shading in Fig. 8d, which is the sum of Fig. 8a, b, and c, shows the anomalously decreased DLWRF over the western and southeastern TP. The originally-decomposed DLWRF perturbations induced by anomalous atmospheric temperature are shown in Fig. 8d (contour), with a decreased magnitude of DLWRF reduction along the western TP to the eastern TP. The difference between the shading and the contour in Fig. 8d presents that TP. This linkage between tropical SST bias and TP SAT in experiment SEN-ALL is also applicable for experiments SEN-PAC and SEN-ALL6.
To further explore which component of the large-scale circulation is responsible for the energy divergence over the TP, 300-hPa heat flux and its divergence are analyzed. In Fig. 10a, b,  , and the sum of the two, respectively. For climatological mean, the meridional component is mainly convergent and the zonal component is divergent due to an accelerated jet stream over the TP and its east ( Fig. 10a and  b). The sum of the heat flux divergence shows a convergence over the eastern TP and a divergence over the western TP (Fig. 10c). Compared with CTRL, the heat flux change exhibits a cyclonic circulation over the western TP, accompanied by a stronger zonal heat transport along 30˚N and a and Fig. 9b, c, and d illustrate the changes in 300-hPa heat flux divergence of SEN-ALL, SEN-PAC, and SEN-ALL6 with respect to CTRL. Under the forcing of tropical SST bias, heat flux exhibits an anomalous divergence over the main body of the plateau, with a small convergence area over the southeastern TP (Fig. 9b). The analogous pattern between Fig. 9a and b strongly suggests that the atmospheric dynamic processes related to heat transport are responsible for the atmospheric cooling leading to a decreased DLWRF and then a colder surface over the TP. Moreover, giving the obvious cold bias over the TP in all SEN-ALL, SEN-PAC, and SEN-ALL6 (Fig. 5m), the linkage between the tropical SST bias and the TP SAT in SEN-ALL may also exist in the other two cases. Figure 9c and d reveal an apparent heat flux divergence over the main body of the TP, which indicates that, when considering the forcing of the tropical SST bias over the Pacific Ocean or the tropical SST bias derived from the CMIP6 MME, there would be an analogous heat flux divergence over the TP, leading to a cooling atmosphere and then a reduced DLWRF and a colder surface over the In short, our physical analysis provides a mechanical understanding of the link of TP cold bias to tropical SST biases. In response to the tropical SST biases, the westerly jet stream tends to intensify above the TP, which leads to more energy transported out of the plateau, causing a cold air bias above the TP. The cold air results in a reduction of downward thermal radiative energy fluxes at the surface, manifesting the TP cold bias.
weaker poleward and equatorward heat transport over the southern TP and the northern TP, respectively (vectors in Fig. 10d, e, and f, for the difference between SEN-ALL and CTRL). Therefore, the change in the meridional component of heat flux divergence shows an apparent divergence over the TP along with its southern flank and a convergence over the western TP (Fig. 10d), while the change in the zonal component presents an anomalous convergence over the southeastern TP and a remarkable divergence over the western TP (Fig. 10e). Overall, the sum of the change in heat flux divergence (Fig. 10f) depicts an obvious energy divergence over the entire TP, except for a small energy convergence over the southeastern TP. Therefore, due to the forcing of tropical SST bias, the change in heat flux shows an anomalous divergence over the TP and the atmospheric cools down consequently, resulting in a decrease in DLWRF at the surface. The changes in heat flux at other levels present Fig. 7 Spatial pattern of partial perturbations of the surface downward LW radiative flux (units: W m − 2 ) in DJF-mean due to a cloud change alone, b water vapor change alone, c air temperature change alone, and d the shading denotes the sum of a-c, the contour denotes total perturbation of the surface DLWRF in SEN-ALL with respect to CTRL. Note that the cloud induced partial perturbation of downward LW radiative surface energy flux is corrected by using the pattern-amplitude projection method (see the text for details). The grey solid lines denote the topographic height of 1500 m annual cycle, with a negative bias over a large portion of the oceans and a small positive bias in the east coast of the mainland. The significant tropical SST bias largely overlaps the high correlation between the TP SAT and the tropical SST, indicating the potential impact of the SST bias on the TP SAT in the models. A series of model experiments have been conducted to understand the influence of remote tropical SST bias forcing on the TP SAT. A surface energy budget analysis shows that, under the forcing of tropical SST bias, the DLWRF decreases significantly, and thus the surface gains less energy and a colder surface appears. The decomposition of the decreased DLWRF indicates that the cooling atmosphere over the TP is the main contributor and that the anomalous heat flux divergence induced mainly by weaker meridional heat flux convergence is responsible for the cooling of the

Conclusions
Compared with observational gridded data, the CMIP5 and CMIP6 models underestimate the SAT over the TP. While previous studies have been focused on local land surface factors such as the overestimation of surface albedo, in this study we explore the role of tropical SST bias as a critical remote external forcing in determining the characteristics of surface energy budget and its influence on the bias of the TP SAT. First, we display the systematic cold bias over the TP and the bias over the tropical oceans in both CMIP5 and CMIP6 models. The cold bias over the TP is larger in the cold season than in the warm season. Compared with the CMIP5, the CMIP6 MME shows an obvious improvement in reducing the bias in summer but with a comparable bias in winter. The systematic tropical SST bias exhibits a small Fig. 8 Spatial pattern of partial surface downward LW radiative fluxes (units: W m − 2 ) in DJF-mean due to the partial air temperature changes whose thermal radiative cooling rates are in balance with a radiative heating rate perturbations induced by cloud change, b radiative heating rate perturbations induced by water vapor change, c nonradiative heating rate perturbations associated with atmospheric advective/convective processes, and d the shading denotes the sum of a-c, the contour denotes perturbation of the surface DLWRF due to air temperature change alone. Note that the atmospheric dynamic processes induced partial perturbation of DLWRF is corrected by using the pattern-amplitude projection method (see the text for details). The grey solid lines denote the topographic height of 1500 m and TP SAT has not been fully clarified. For example, how the tropical SST bias causes the change in the subtropical westerly jet stream needs to be better explained. One may also notice that there is a strong annual cycle of the TP SAT bias response to the tropical SST bias forcing which exhibits a neglectable annual cycle. The meridional migration of the climatological subtropical westerly jet stream may be responsible for the distinct TP responses between summer and winter. The climatological westerly jet stream is more or less located over the TP in winter but moves to the north of the TP in summer. The westerly anomalies above the TP in response to the tropical SST anomalies may exert a stronger impact in winter because they are directly superimposed on the westerly jet core. Furthermore, while the collective contributions of SST biases in all three tropical basins to atmosphere. Importantly, the tropical Pacific SST bias forcing from the CMIP5 MME and the all-basin tropical SST bias forcing derived from the CMIP6 MME can both generate the linkage discussed above.
In short, climate models suffer from systematic tropical SST bias that influences heat transfer, leading to a cooled atmosphere and a reduced DLWRF and thus a colder surface over the TP. The tropical SST bias in the CMIP5 and CMIP6 models explains about 40% of the total cold bias over the TP.
It should be pointed out that although we have established a strong relationship between tropical SST and TP SAT and explained the response of TP SAT bias to the SST bias via a series of model experiments and diagnostic analyses, the physical process connecting the biases of SST Fig. 9 a Spatial pattern of partial downward LW radiative surface energy fluxes (units: W m − 2 ) in DJF-mean due to the partial air temperature changes whose thermal radiative cooling rates are in balance with the nonradiative heating rate perturbations associated with atmospheric advective/convective processes, which is the same as in Fig. 8c with a different color bar. b Changes in DJF-mean horizontal heat flux divergence (units: K day − 1 ) at 300 hPa in SEN-ALL with respect to CTRL. c Changes in DJF-mean horizontal heat flux divergence (units: K day − 1 ) at 300 hPa in SEN-PAC with respect to CTRL. d Changes in DJF-mean horizontal heat flux divergence (units: K day − 1 ) at 300 hPa in SEN-ALL6 with respect to CTRL. The grey solid lines denote the topographic height of 1500 m indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless