**3. Results**

#### *3.1. Dominant Modes of Variation in Summer SH over the TP*

Figure 1a,b shows the spatial pattern of summer SH over the TP during the period 2001–2018 obtained from ERA-Interim reanalysis (Figure 1a) and the remote sensing product (Figure 1b). In general, ERA-Interim and remote sensing-based SH presented similar spatial patterns, with a pattern correlation of 0.8, exceeding the 99% confidence test. The values of SH were positive over the whole TP, acting as a gigantic SH air pump and having crucial climatic effects [9]. SH increased from south to north, with the maximum value in the Qaidam Basin over the northern TP. Figure 1c,d shows the spatial distributions of the standard deviation in SH derived from ERA-Interim reanalysis (Figure 1c) and the remote sensing product (Figure 1d). The spatial pattern of ERA-Interim SH standard deviation was close to the remote sensing-based SH standard deviation. In the two datasets, values of the standard deviation in SH over the northern and southern TP were also larger than those in the eastern TP, which were similar to the spatial distributions of summer seasonal-mean SH.

To identify the dominant modes of summer SH anomalies, EOF analysis was performed. The spatial distribution of the first dominant mode (EOF1) was marked by the consistent variations of summer SH over the whole TP, with an explained variance of 20.1% (Figure 2a). EOF1 was set apart from the other modes based on the North test. This means that EOF1 is considered statistically distinguishable and significant. Figure 2b shows the spatial distribution of EOF2. EOF2 accounted for 14.2% of the total variance and passed the North test [32]. A clear zonal seesaw pattern could be observed over the TP. The same analysis based on the remote sensing product derived from Han et al. [28,29] also indicated

that the first leading mode of summer SH showed significant consistent variations, and the second leading mode of summer SH showed a zonal dipolar pattern (Figure S1).

**Figure 1.** The spatial pattern of the TP summer SH climatology (W m<sup>−</sup>2) from ERA-Interim (**a**) and remote sensing product derived from Han et al. [28,29] (**b**), standard deviation (W m<sup>−</sup>2) from ERA-Interim reanalysis (**c**), and remote sensing (**d**) data.

Figure 2c shows the PC1 and its interdecadal component (PC1-ID). PC1-ID presented an interdecadal shift from the positive phase to the negative phase from around 1996. From the perspective of long-term trends, SH demonstrated reductions from the 1980s and a slight recovery from the 2000s, which is in agreemen<sup>t</sup> with previous studies [24,25]. The spatial distribution of the summer SH interdecadal component anomalies onto PC1-ID is shown in Figure 2e. Significant positive anomalies appeared in almost the entire TP, while weak negative anomalies occurred in the southeastern TP, which was similar to the first leading mode of SH driving from the EOF method, with a pattern correlation of up to 0.96 (Figure 2a). This indicates that the PC1-ID can effectively represent the interdecadal variation of summer SH over the TP.

Figure 2d shows the PC2 and its interannual component (PC2-IA). A remarkable feature of this was that the relation between PC2-IA and PC2 had a correlation coefficient of up to 0.95, which indicates that EOF2 presents the interannual variation of summer SH. The spatial distribution of SH anomalies associated with PC2-IA is shown in Figure 2f. A clear zonal seesaw pattern could be observed over the TP, with negative SH in the east of the TP and positive SH over the western TP. On the interannual time scale, when the PC2-IA was in a positive phase, SH over the TP exhibited a zonal asymmetric spatial mode with anomalous strengthening in SH over the western TP and anomalous weakening in SH over the eastern TP. Anomalies in the interannual component of SH associated with PC2-IA were the same as the second leading mode of SH derived from the EOF method (Figure 2b), which indicates that the EOF2 pattern can effectively represent the dominant mode of interannual variation in summer SH over the TP.

**Figure 2.** The first two EOFs of summer SH over the TP for the period 1981–2018 (**<sup>a</sup>**,**b**). The explained variance of EOF is in the upper right of (**<sup>a</sup>**,**b**). Normalized time series of PC1 (column) and its interdecadal component (PC1-ID) (line) (**c**). Normalized time series of PC2 (column) and its interannual component (PC2-IA) (line) (**d**). Regression of the interdecadal component of summer SH (shading, W m<sup>−</sup>2) on PC1-ID (**e**). Regression of the interannual component of summer SH (shading, W m<sup>−</sup>2) on PC2-IA (**f**). Grid points with statistically significant anomalies passing the 95% confidence level are denoted by an oblique line.

### *3.2. Physical Mechanisms of Variations in Summer SH*

#### 3.2.1. First Dominant Mode of Variation in Summer SH

The first dominant mode of variation in summer SH presented an interdecadal weakening. As SH was determined by the local surface energy budget, the local surface energy budget related to anomalies in SH on the interdecadal time scale is investigated below. Considering SH experienced a decadal phase shift around 1996, we used the two periods, 1981–1995 and 1997–2018, to represent the decadal change in SH accordingly. Figure 3 presents the difference in surface fluxes between the two periods. Anomalies in SWR bore a close resemblance to the spatial pattern of SH with negative anomalies over the majority of the TP (Figure 3a,c). The anomalous decrease in SWR corresponded to significantly reduced DSWR over the TP (Figure 3e). This indicates that the DSWR has a positive contribution to the interdecadal variation of SH. The spatial distribution of anomalies in LWR presented strengthening over the TP (Figure 3d), which was dominated by the change in DLWR (Figure 3g). Moreover, the magnitude of anomalies in LH, USWR, ULWR, and G was too small to impact SH (Figure 3b,f,h,i). The contributions of the components to changes in SH are depicted in Figure 4. It was found that the first dominant mode of SH was mainly due to anomalous SWR induced by DSWR, while the LWR anomalies induced by DLWR

helped to offset the decadal weakening of SH. The contribution of the other components was insignificant.

**Figure 3.** Spatial distribution of difference in sensible heating (SH) (**a**), latent heating (LH) (**b**), surface net shortwave radiation (SWR) (**c**), surface net longwave radiation (LWR) (**d**), downwelling surface shortwave radiation (DSWR) (**e**), upwelling surface shortwave radiation (USWR) (**f**), downwelling surface longwave radiation (DLWR) (**g**), upwelling surface longwave radiation (ULWR) (**h**), and ground heating (G) (**i**) between 1997–2018 and 1981–1995 (shading, W m<sup>−</sup>2). Grid points with statistically significant anomalies passing the 95% confidence level are denoted by an oblique line.

**Figure 4.** The difference in SH, LH, SWR, LWR, DSWR, LSWR, DLWR, ULWR, and G between 1997–2018 and 1981–1995 over the TP (W m<sup>−</sup>2).

According to previous studies, cloud cover can regulate the DSWR and DLWR [33]. Therefore, the opposite variations of DSWR and DLWR may have been induced by the cloud cover. The total cloud cover variations associated with the interdecadal variation of SH were investigated next. Figure 5a displays the difference in the total cloud cover between 1997–2018 and 1981–1995. The TP was dominated by enhanced cloud cover, which induced more reflection of downwelling solar radiation and resulted in reduced DSWR. A decrease in DSWR can enhance SH owing to the surface energy budget. Figure 5b depicts the difference in water vapor flux and divergence between 1997–2018 and 1981–1995. A large westward water vapor flux emerged over the TP. This can enhance the convergence of water vapor flux over the TP by means of decreasing the water vapor exported from the eastern boundary of the TP. Associated with this increase in water vapor convergence, more summer cloud cover appeared over the TP (Figure 5a), inducing a decadal weakening of SH (Figure 3a).

**Figure 5.** Spatial distribution of difference in total cloud cover (**a**), water vapor flux exceeding the 95% confidence level (vector, kg m<sup>−</sup><sup>1</sup> s<sup>−</sup>1), and divergence (shading, 10−<sup>4</sup> kg m<sup>−</sup><sup>2</sup> s<sup>−</sup>1) (**b**) between 1997–2018 and 1981–1995. Grid points with statistically significant anomalies passing the 95% confidence level are denoted by oblique lines.

Large-scale circulation anomalies can usually contribute to regional climatic change. Studying the impact of atmospheric circulation on TP SH can improve our understanding of climatic change over the TP, which is helpful for climate prediction. The atmospheric circulation anomalies linked to the interdecadal weakening of summer SH over the TP were examined further. Figure 6a displays the 500 hPa geopotential height and horizontal wind anomalies between 1997–2018 and 1981–1995. There was a wave train in the Eurasian continent. A considerably positive geopotential height center with an anticyclonic anomaly was located to the northeast of the TP, corresponding to this anomalous wave-like train. This led to a weakened subtropical westerly jet over the TP and less water vapor export from the eastern boundary of the TP. Figure 6b shows the 200 hPa geopotential height anomalies during the two periods. A zonal wave train occurred from the North Atlantic to the Eurasian continent, with two negative-height anomalies over the North Atlantic and West Asia and two positive-height anomalies over eastern Europe and Lake Baikal. The wave activity flux indicated that the wave train's propagation direction was expected to be from the North Atlantic to the Eurasian continent. Such a circulation anomaly pattern is favorable for a reduced westerly jet and led to anomalous water vapor convergence of the TP (Figure 5b). This conclusion coincides with Zhou et al. [34], who proposed that the interdecadal variation in summer water vapor over the TP is related to a similar teleconnection pattern. Furthermore, the atmospheric circulation anomalies linked to PC1-ID were similar to those related to the interdecadal variation of SRP [35], implying that the SRP pattern may be crucial in the interdecadal variability of summer SH over the TP.

**Figure 6.** Spatial distribution of difference in 500 hPa geopotential height anomalies (shading, gpm) and the wind exceeding the 95% confidence level (vector, m s<sup>−</sup>1) (**a**), the 200 hPa geopotential height anomalies (shading, gpm), and the wave activity flux (vector, m<sup>2</sup> s<sup>−</sup>2) (**b**) between 1997–2018 and 1981–1995. Grid points with statistically significant anomalies passing the 95% confidence level are denoted by oblique lines.

### 3.2.2. Second Dominant Mode of Variation in Summer SH

We next investigated the surface flux anomalies related to the second dominant mode of variation in TP summer SH. Figure 7 displays the surface flux anomalies related to the PC2-IA of summer SH. The upwelling surface longwave and shortwave radiations ranged from the surface to the atmosphere. In general, a large decrease in surface latent heating (LH) mainly appeared in the southern TP, and an insignificant increase occurred in the northern TP. Corresponding to the EOF2 of SH, a significant increase in the surface net shortwave radiation (SWR) was observed over the western TP and a remarkable reduction appeared in the eastern TP (Figure 7c). Anomalies in SWR were dominated by the zonally asymmetric pattern of downwelling surface shortwave radiation (DSWR) (Figure 7e). However, the contribution of upwelling surface shortwave radiation (USWR) to anomalies

in SWR was insignificant (Figure 7f). This means that the DSWR has a positive contribution to anomalies in SH over the TP. On the contrary, anomalies of surface longwave radiation (LWR) (Figure 7d) mainly induced by anomalous downwelling surface longwave radiation (DLWR) (Figure 7g) negatively affected SH due to the opposite spatial pattern anomalies. Moreover, anomalies in upwelling surface longwave radiation (USWR) and the ground heat flux (G) were insignificant and had little impact on the change in SH.

**Figure 7.** Anomalies in the interannual component of SH (**a**), LH (**b**), SWR (**c**), LWR (**d**), DSWR (**e**), LSWR (**f**), DLWR (**g**), ULWR (**h**), and G (**i**) regressed on PC2-IA (shading, W m<sup>−</sup>2). Grid points with statistically significant anomalies passing the 95% confidence level are denoted by oblique lines. The two black rectangles denote the western TP (28–36◦ N, 78–88◦ E) and eastern TP (28–36◦ N, 92–102◦ E).

To further examine the contribution of each surface flux to SH anomalies, the average regional anomalies of SH, LH, SWR, LWR, DSWR, USWR, DLWR, ULWR, and G were calculated for the eastern TP (28–36◦ N, 78–88◦ E) (Figure 8a) and western TP (28–36◦ N, 92–102◦ E) (Figure 8b), respectively. For the eastern TP, a large decrease in SWR of approximately −9Wm−<sup>2</sup> contributed to the decrease in SH and LH. The decrease in SWR is mainly attributed to the weakening of DSWR. For the western TP, enhanced SH was consistent with the increase in DSWR and decrease in LH, which means that the increase in SH was mainly due to DSWR and LH. The LWR presented a significant reduction of approximately −2.2 W m<sup>−</sup>2. Among the contributions to the SH change, the role of DSWR was predominant. Therefore, we investigated the role of DSWR in changing SH on the interannual time scale.

**Figure 8.** Anomalies in the interannual component of SH, LH, SWR, LWR, DSWR, LSWR, DLWR, ULWR, and G regressed on PC2-IA (W m<sup>−</sup>2) in the region of eastern TP (**a**) and western TP (**b**).

As investigated by previous studies [33], cloud cover has an important influence on the shortwave radiation reaching the surface. Figure 9a shows anomalies in total cloud cover related to the interannual variation of the TP summer SH. In contrast, less cloud cover appeared in the western TP, and more cloud cover appeared in the eastern TP. This indicates that the anomalous DSWR leading to the zonally asymmetric pattern of SH may be partly related to the variations of total cloud cover over the TP. Figure 9b presents the regression of water vapor flux and divergence onto PC2-IA. There was a significant southwesterly vapor supply from the Indian monsoon region to the main body of the TP, which was induced by the southwesterly anomalies. The anomalous water vapor flux divergence also presented a dipole mode with divergence in the west of the TP and convergence in the east of the TP. Sufficient moisture convergence (divergence) is favorable for the increase (decrease) in total cloud cover, and therefore induced weakened (enhanced) SH over the eastern (western) TP. According to previous studies, the variations of SH significantly respond to the changes in surface wind speed [17]. Figure 9c also displays the regression of the surface wind speed on PC2-IA. Remarkable southwesterly anomalies and positive surface wind speed anomalies were located in the western TP, in accordance with the anomalous enhancement in SH of the western TP (Figure 7a). This means that the positive SH of the western TP was also partly due to the strengthened surface wind speed, owing to the southwesterly wind anomalies.

**Figure 9.** Anomalies in the interannual component of total cloud cover (**a**), water vapor flux exceeding the 95% confidence level (vector, kg m<sup>−</sup><sup>1</sup> s<sup>−</sup>1) and divergence (shading, 10−<sup>4</sup> kg m<sup>−</sup><sup>2</sup> s<sup>−</sup>1) (**b**), surface wind vectors passing the 95% confidence level (vector, m s<sup>−</sup>1), and wind speed (shading, m s<sup>−</sup>1) (**c**) regressed on PC2-IA. Grid points with statistically significant anomalies passing the 95% confidence level are denoted by an oblique line.

The above result indicates that interannual variation in SH can be mainly explained by the DSWR change, which is associated with moisture convergence and divergence. In addition, the southwesterly wind anomaly can also partly explain the SH anomalies in the western TP. These analyses sugges<sup>t</sup> that atmospheric circulation changes may have an important impact on the SH anomalies by impacting the water vapor convergence and divergence and wind anomalies over the TP. Atmosphere circulation anomalies related to the anomalies in SH were further examined. Figure 10a displays 500 hPa geopotential height and horizontal wind anomalies related to the interannual variation of SH. An evident feature is that there was a significant, negative, high-pressure center to the west of the TP corresponding to the cyclonic anomaly, which led to enhanced southwest wind anomalies at the southwestern boundary of the TP and water vapor convergence (divergence) over the eastern (western) TP. Furthermore, 200 hPa geopotential height regressed against normalized PC2-IA presented a similar spatial pattern of 500 hPa, showing an equivalent barotropic structure in the vertical direction (Figure 10b). According to the horizontal wave activity flux, the teleconnection wave associated with SH anomalies originated in the North Atlantic and traveled eastward across the Eurasian continent. This suggests that atmospheric wave-like patterns originating in the North Atlantic may be important in the interannual variation of SH.

**Figure 10.** Anomalies in the interannual component of 500 hPa geopotential height (shading, gpm) and wind passing the 95% confidence level (vector, m/s) (**a**), the 200 hPa geopotential height (shading, gpm), and the wave activity flux (vector, m<sup>2</sup> s<sup>−</sup>2) (**b**) regressed on PC2-IA. Grid points with statistically significant anomalies passing the 95% confidence level are denoted by oblique lines.

3.2.3. Association between Variations in Summer SH and Atmospheric Wave Trains

The above analysis indicated that the dominant modes of variations in summer SH over the TP had a close connection to atmospheric wave trains. Associations between variations in summer SH and atmospheric wave trains are investigated below.

We first investigated the impact of interdecadal variation in SRP on summer SH. According to previous studies [35–37], the SRP pattern is defined as the EOF1 of the summer mean 200 hPa meridional wind within the region 20–60◦ N, 0–150◦ E (Figure 11a), and the normalized PC1 is referred to as the SRP index (SRPI) (Figure 11b). A 9-year Lanczos low-pass filter was used to calculate the interdecadal component of SRPI (SRPI-ID), which can be seen in Figure 11c. The temporal correlation coefficient between PC1-ID and the SRPI-ID was up to 0.84, with a significance level of 99% (Figure 11c). Figure 11d shows the spatial distribution of SH anomalies related to the interdecadal variation of SRP. Significantly decreased anomalies occurred in almost the whole TP, while weak negative anomalies appeared in the southeastern TP. The spatial pattern was in good agreemen<sup>t</sup> with the distribution of interdecadal weakening of SH (Figure 3a). This suggests

that the SRP is closely related to the interdecadal variability of SH. Figure 11e,f depict atmospheric circulation anomalies linked to SRP interdecadal variation. The 500 hPa geopotential height and horizontal wind anomalies regressed upon the SRPI-ID multiplied by −1 are presented in Figure 11e. To the northeast of the TP, significant positive height anomalies arose, which were related to the decadal weakening of SH. The anomalous positive height accompanied the anticyclonic anomaly to the northeast of the TP. In relation to the anticyclonic circulation at 500 hPa, easterly wind anomalies to the south of the anticyclonic center flowed into the TP from the eastern boundary of the TP, resulting in a diminished subtropical westerly jet. Corresponding to the easterly anomalies there was significant water vapor convergence and enhanced cloud cover, which could reduce DSWR and SH in the TP. The above analysis indicated that the interdecadal variation of SRP can significantly induce anomalous weakening in SH over the TP by impacting the circulation anomalies and, therefore, cloud cover and DSWR anomalies.

**Figure 11.** (**a**) The first EOF of summer means 200 hPa meridional wind over the region 20–60◦ N, 0–150◦ E for the period 1981–2018. (**b**) Normalized time series of SRPindex (column) and its interdecadal component (line). (**c**) Normalized time series of interdecadal component of SRP index (column) and PC1-ID (line). (\*\*: significant at 99% confidience level). Anomalies in the interdecadal component of sensible heating (**d**), 500 hPa geopotential height (shading, gpm) and the wind exceeding the 95% confidence level (vector, m/s) (**e**), the 200 hPa geopotential height (shading, gpm) (**f**) regressed on SRPI-ID multiplied by −1. Grid points with statistically significant anomalies passing the 95% confidence level are denoted by oblique lines.

We also found that the second leading mode of the upper-tropospheric 200 hPa meridional wind anomalies in 20–60◦ N, 0–150◦ E during summer, which is defined as the North Atlantic-East and North Asia pattern (NAENA), had significant impacts on the interannual variation of SH (Figure 12a). The NAENA index (NAENAI) is defined as the PC2 corresponding to EOF2 of the summer 200 hPa meridional wind over the Eurasian continent (Figure 12b). Figure 12c displays the normalized time series of the interannual component of NAENAI (NAENAI-IA) (column) and PC2-IA (black line). The temporal correlation coefficient between PC2-IA and NAENAI-IA during the period of 1981–2018 was up to 0.63, which indicated that the interannual variation in NAENA could explain approximately 40% of the anomalies in summer SH on the interannual time scale. Figure 12d shows anomalies in SH obtained by regression upon NAENAI. Positive SH anomalies predominated in the west of the TP, whereas weakened SH anomalies arose in the east of the TP, resembling the spatial pattern of EOF2 (Figure 2b). Moreover, the spatial pattern of NAENA showed close similarity to the atmospheric circulation anomalies associated with the interannual variation of SH (Figure 12e,f). The above analysis suggested that NAENA can impact the interannual variation of summer SH over the TP by means of changing the wind speed and water vapor convergence and divergence over the TP.

**Figure 12.** (**a**) The second EOF of summer mean 200 hPa meridional wind in 20–60◦ N, 0–150◦ E for the period 1981–2018. (**b**) Normalized time series of NAENA index (column) and its interannual component (line). (**c**) Normalized time series of interannual components of NAENAI (column) and PC2-IA (line), (\*\*: significant at 99% confidience level). Anomalies in the interannual component of sensible heating (**d**), 500 hPa geopotential height (shading, gpm) and the wind exceeding the 95% confidence level (vector, m/s) (**e**), the 200 hPa geopotential height (shading, gpm), and the wave activity flux (vector, m<sup>2</sup> s<sup>−</sup>2) (**f**) regressed on PC2-IA. Grid points with statistically significant anomalies passing the 95% confidence level are denoted by oblique lines.
