*2.1. Data*

To study the variations of summer SH over the TP, ERA-interim surface sensible heat flux datasets were used in this paper, which were available at a horizontal resolution of 0.5◦ × 0.5◦ from 1981 to 2018. To improve the handling of data bias and background error constraints, ERA-interim uses a 12-h 4D-Varassimilation system. Compared to ERA-40, ERA-Interim has additional input remote data, including Meteosat-2 clear-sky radiances, Global Ozone Monitoring Experiment (GOME) ozone profiles, and radio occultation measurements from the Challenging Mini Satellite Payload (CHAMP), Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC), and Gravity Recovery and Climate Experiment (GRACE). Reprocessed ocean wave height data from ERS-1 and ERS-2, as well as upper-level winds from Meteosat-2, were also included. Previous investigations have shown that ERA-interim SH data can accurately reflect the surface heat flux [27]. Element fields such as the surface wind speed, total cloud cover, surface latent heating, and surface radiation data are also provided by ERA-interim.

A remote-sensing product derived from Han et al. [28,29] was also used to analyze summer SH on the TP. A detailed retrieval algorithm of SH can be found in Han et al. [28,29].

To quantify the association between SH and atmospheric circulation, we used the 38-year (1981–2018) monthly mean geopotential height, zonal wind, and meridional wind use from the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP/NCAR), with a horizontal grid spacing of 2.5◦ × 2.5◦. The data we used were for the summer season (JJA). Table 1 shows the details of each variables used in this study.

**Table 1.** Summary of the datasets used in the study.


*2.2. Methods*

2.2.1. Surface Energy Balance Analysis

Surface energy balance analysis was employed to identify the cause of change in SH. It is expressed as follows [30]:

$$\rm{SW}\_{\downarrow} - \rm{SW}\_{\uparrow} + \rm{LW}\_{\downarrow} - \rm{LW}\_{\uparrow} = \rm{SH} + \rm{LH} + \rm{G}\_{\prime} \tag{1}$$

where SW↓ and SW↑ indicate the downwelling and upwelling shortwave solar radiation, respectively; LW↓ and LW↑ indicate the downwelling and upwelling longwave thermal infrared radiation, respectively; SH and LH indicate the upward surface sensible and latent heating; and G is the ground heating. The data were obtained from ERA-Interim. Therefore, the change in SH was determined by the change in the other terms in the equation, which can be expressed as follows:

$$
\Delta \text{SH} = \Delta \text{SW}\_{\downarrow} - \Delta \text{SW}\_{\uparrow} + \Delta \text{L} \text{W}\_{\downarrow} - \Delta \text{L} \text{W}\_{\uparrow} - \Delta \text{L} \text{H} - \Delta \text{G}\_{\prime} \tag{2}
$$

### 2.2.2. Wave Activity Flux

The wave activity flux formulated by Takaya and Nakamura [31] was used to determine the features of stationary wave propagation, which is given as follows:

$$\mathcal{W} = \frac{1}{2|\overline{\mathcal{U}}|} \left\{ \begin{array}{l} \overline{\pi} (\boldsymbol{\psi}^{\prime}\_{\boldsymbol{x}} - \boldsymbol{\psi}^{\prime} \boldsymbol{\psi}^{\prime}\_{\boldsymbol{x}\boldsymbol{x}}) + \overline{\pi} \Big(\boldsymbol{\psi}^{\prime}\_{\boldsymbol{x}} \boldsymbol{\psi}^{\prime}\_{\boldsymbol{y}} - \boldsymbol{\psi}^{\prime} \boldsymbol{\psi}^{\prime}\_{\boldsymbol{x}\boldsymbol{y}}\Big) \\ \overline{\pi} (\boldsymbol{\psi}^{\prime}\_{\boldsymbol{x}} \boldsymbol{\psi}^{\prime}\_{\boldsymbol{y}} - \boldsymbol{\psi}^{\prime} \boldsymbol{\psi}^{\prime}\_{\boldsymbol{x}\boldsymbol{y}}) + \overline{\pi} (\boldsymbol{\psi}^{\prime 2}\_{\boldsymbol{y}} - \boldsymbol{\psi}^{\prime} \boldsymbol{\psi}^{\prime}\_{\boldsymbol{y}\boldsymbol{y}}) \end{array} \right\} \tag{3}$$

Here, *ψ* is the stre-m function, *u* is the zonal wind v-locity, and *v* denotes the meridio-al wind velocity. The data were obtained from NCEP/NCAR. The overbars represent the basic states and primes represent perturbations.

### 2.2.3. Empirical Orthogonal Functions (EOF) Analysis

EOF analysis was used to study the features of interannual and interdecadal variations of summer SH. In addition, the EOF analysis was also used to extract the teleconnection wave trains related to variations of SH. The EOF method is frequently used to investigate potential spatial patterns of climatic variability and how they develop over time. In EOF analysis, the original climatic data were also projected on an orthogonal basis [32]. Moreover, the orthogonal basis was determined by calculating the eigenvector of the spatially weighted anomaly covariance matrix, with the corresponding eigenvalues indicating the percentage variance explained by each pattern. Therefore, the EOFs of spatiotemporal physical processes could represent mutually orthogonal spatial patterns in the data change set, in which the first pattern accounts for most of the variance, the second pattern accounts for most of the residual variance, etc. As the principal component (PC) of an EOF mode shows how the spatial pattern of this mode oscillates over time, we used the corresponding PC of the dominant mode as the reference time series for each summer of SH.

### 2.2.4. Linear Regression Analysis and Composite Analysis

To investigate the anomalies in surface energy flux and atmospheric circulation associated with SH, linear regression and composite analysis were employed in this study. Using the statistical analysis method of linear regression analysis, the linear relationship between two or more variables can be determined quantitatively. A composite analysis is frequently used in climate change research to explore the salient characteristics of different periods. The positive phase period of SH and the negative phase period were selected first. Second, we investigated the atmospheric circulation differences between the two periods. The significance of the regression coefficient and difference in the composite analysis was evaluated using Student's *t*-test.
