Observed Multi-Timescale Di ﬀ erences between Summertime Near-Surface Equivalent Temperature and Temperature for China and Their Linkage with Global Sea Surface Temperatures

: Based on the ensemble empirical mode decomposition method, this study explores the di ﬀ erences and similarities in multiple time-scale characteristics of summer air temperature ( T ) and equivalent temperature ( T e ) over China during 1961–2017, using daily meteorological observations collected at 412 stations in China. Their relationships to global sea surface temperature variations is also discussed. Results show that both T and T e can be decomposed into ﬁve components, which includes multiple timescales, from interannual to long-term trends. The spatial patterns of each timescale’s leading mode show that the variations of T e are generally larger than that of T . Meanwhile, both T and T e are dominated by their inter-annual, multi-decadal variations and the non-linear trend. High correlations of T and T e can also be found in these major scales. The related sea surface temperature variations in these major scales also show consistent patterns, which correspond to El Niño–Southern Oscillation, Atlantic Multidecadal Oscillation and the global warming trend in the sea, respectively. In other scales, both spatial patterns of T and T e and the corresponding correlation patterns with sea surface temperature are distinguishable. The current results explore the compound changes of surface temperature-humidity during the past ﬁve decades from a new perspective, which provides some insights for a better understanding of the possible causes of climate change over China.


Introduction
Global warming during the past century is an unequivocal fact [1,2]. Based on the assessment of the latest report by the Intergovernmental Panel on Climate Change (IPCC) on the global warming by 1.5 • C, human-induced warming has exceeded 1 • C in 2017 above the preindustrial period (1750-1850) at a warming rate of approximately 0.2 • C/10a [3]. Many studies have indicated that the acceleration of global warming will lead to more hot extremes and fewer cold extremes [4,5], presenting human society and economic development with unprecedented challenges.
Atmosphere 2019, 10 Currently, most studies use a single climatic factor to characterize or define climate change, global warming or extreme events. For instance, surface air temperature (T) is the most common-used straightforward temperature metric that represents climate change, variability and trend. However, climate change or the occurrence of extreme events may not be appropriately presented by a change in T alone, and it is likely to be a compound change of multiple variables (e.g., the joint change of temperature and humidity) [6][7][8][9][10][11][12][13]. Furthermore, the changes and causes of joint changes in multiple variables are more complicated than a single variable, and the compound extreme changes of multiple variables can amplify the impact of independent events, resulting in more severe consequences [14,15]. In particular, for heatwaves (or heat extremes), a growing number of studies in recent years have highlighted the advantages and importance of using the multi-variable-based temperature metric (e.g., equivalent temperature) [6,8,12]. For example, Pal and Elthair [12] stated that the warmings of wet-bulb temperature (an index that represents the joint temperature-humidity variations) over southwest Asia in future scenarios (i.e., RCP 4.5 and RCP 8.5) are likely to exceed the critical threshold of the human capacity to adapt to climate change, which will possibly bring deadly threats to human health. Pielke et al. [6] found that there are large differences between T and moist enthalpy (ME)-based equivalent temperature (T e ) for 2002, especially when the absolute humidity is high. Hence, he noted that a more comprehensive multivariable heat metric (i.e., ME-based T e ) should be used to comprehensively measure global surface heat content change. Fall et al. [8] revealed that the moisture component of T e induces larger trends and variability relative to T during 1979-2005 over the United States. Furthermore, their work also shows that T e can be used as a complementary variable to assess the land-atmosphere interaction processes. In recent years, T e is becoming more and more widely used in climate research due to it being easy-to-compute, having clear physical meanings and the significance of its representation (e.g., Pielke et al. [6], Fall et al. [8], Kang and Eltahir [16] and Younger et al. [17]). However, the characteristics of T e during recent decades, especially the differences between T and T e , remain to be investigated and fully understood.
In recent decades, heat waves in summer have become increasingly frequent in China [18]. Numerous studies have investigated the characteristics and causes of variations in air temperature (T) during summer. For example, Sun et al. [19] indicate that since the early 1950s, greenhouse gas emissions caused by human activities have led to an increase in extreme high-temperature events during summer in China. Taking the 2013 extreme summer heatwave as an example, the possibility of that has increased more than 60-fold due to anthropogenic influences [19]. Some researchers have also found that urbanization could contribute to over 20% of the total warming trend in summer air temperature and lead to significant increases (decreases) in the related hot (cold) extremes during the past half-century [20,21]. On multidecadal scales, Qian et al. [22] noted that both the frequency of extreme heat events and the evolution of T in summer over Shanghai had clear multidecadal variability for a 60-to 80-year period, which was probably linked with the Atlantic Multidecadal Oscillation (AMO). On the interdecadal scale, Wu et al. [23] found that the first empirical orthogonal function (EOF) mode of the heat waves in summer in China is closely related to the snow cover in the western part of the Tibetan Plateau. Hu et al. [24] demonstrated that on an interannual scale, the sea surface temperature (SST) anomalies in the Indian Ocean Basin affected the summer T in southern and northeastern China. Wu et al. [25] illustrated that the variation in T in northeastern China is related to the El Niño-Southern Oscillation (ENSO). Despite these pre-existing studies, for the compound temperature-humidity changes, their multiple timescales characteristics and differences from those of the mean air temperature remain unclear. Therefore, based on the above considerations, the primary scientific questions discussed in this paper are as follows: (1) What are the similarities and differences in the characteristics between T and T e at multiple timescales in China in recent decades? (2) What is the relationship between the multiple timescale changes of the two temperatures (i.e., T and T e ) and the global SSTs?

Observational and Reanalysis Datasets
Daily surface air temperature, relative humidity and pressure data from 824 meteorological stations of China's surface climatological dataset (Version 3.0) (http://data.cma.cn/data/cdcdetail/dataCode/ SURF_CLI_CHN_MUL_DAY_V3.0.html) were used to calculate T e . This dataset has been released by the National Meteorological Information Center (NMIC) after performing complete quality control processes (e.g., unusually large values check, spatial and temporal consistency check). The quality and completeness of the data for each variable were significantly improved compared to the previously published surface data of the same type. To minimize any inhomogeneity of the data caused by nonclimatic factors (e.g., relocation of stations, replacement of equipment) and missing data, screening of the 824 stations was conducted based on two criteria: (1) Following Chen and Zhai [26], we chose 20 km and 50 m as the relocation distance thresholds in the horizontal and vertical direction, respectively. The station relocation distance that exceeded either threshold was marked as a significantly relocated station, and then we excluded those marked stations to avoid the inhomogeneity due to the significant station relocations.
(2) For each variable, the exclusion of stations that had missing data for more than 5% of summer days (i.e., June-August, 92 days in total) in any year during 1961-2017.
Ultimately, 412 meteorological stations were retained in the study.
To further reduce the possible uncertainties that were introduced by the inhomogeneity of the stations' geographic distribution, we followed Zhang et al. [27] and Chai et al. [28] to conduct a redistribution algorithm after the screening of the stations. In detail, this process redistributed the station data onto a 0.5 • × 0.5 • equal longitude-latitude grid, and each grid box only contained one station (if any grid contained more than one station, the station data that was last processed would represent the final value of that grid). The advantage of this method is to make sure that the density (or weight) of each grid cell is homogeneous. Finally, 408 equal-weighted grid boxes were retained for further analyses (Figure 1b).

Observational and Reanalysis Datasets
Daily surface air temperature, relative humidity and pressure data from 824 meteorological stations of China's surface climatological dataset (Version 3.0) (http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html) were used to calculate Te. This dataset has been released by the National Meteorological Information Center (NMIC) after performing complete quality control processes (e.g., unusually large values check, spatial and temporal consistency check). The quality and completeness of the data for each variable were significantly improved compared to the previously published surface data of the same type. To minimize any inhomogeneity of the data caused by nonclimatic factors (e.g., relocation of stations, replacement of equipment) and missing data, screening of the 824 stations was conducted based on two criteria: (1) Following Chen and Zhai [26], we chose 20 km and 50 m as the relocation distance thresholds in the horizontal and vertical direction, respectively. The station relocation distance that exceeded either threshold was marked as a significantly relocated station, and then we excluded those marked stations to avoid the inhomogeneity due to the significant station relocations.
(2) For each variable, the exclusion of stations that had missing data for more than 5% of summer days (i.e., June-August, 92 days in total) in any year during 1961-2017.
Ultimately, 412 meteorological stations were retained in the study.
To further reduce the possible uncertainties that were introduced by the inhomogeneity of the stations' geographic distribution, we followed Zhang et al. [27] and Chai et al. [28] to conduct a redistribution algorithm after the screening of the stations. In detail, this process redistributed the station data onto a 0.5° × 0.5° equal longitude-latitude grid, and each grid box only contained one station (if any grid contained more than one station, the station data that was last processed would represent the final value of that grid). The advantage of this method is to make sure that the density (or weight) of each grid cell is homogeneous. Finally, 408 equal-weighted grid boxes were retained for further analyses (Figure 1b). The monthly mean global SST data from the Hadley Center sea ice and SST dataset (HadISST, can be downloaded from https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html) were used to analyze the key SST regions that are possibly linked to T and Te. The horizontal resolution of the data was 1° × 1°. SST indices, including the AMO, Niño 3.4 index and the Pacific Decadal Oscillation (PDO), were obtained from the United States' National Oceanic and Atmospheric Administration (NOAA), which can be downloaded from https://www.esrl.noaa.gov/psd/data/climateindices/list. North Atlantic SST tripole (NAT) index can be downloaded from https://cmdp.ncc-cma.net/download/precipitation/diagnosis/NAT/NAT-monnor.tms. The monthly mean global SST data from the Hadley Center sea ice and SST dataset (HadISST, can be downloaded from https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html) were used to analyze the key SST regions that are possibly linked to T and T e . The horizontal resolution of the data was 1 • × 1 • . SST indices, including the AMO, Niño 3.4 index and the Pacific Decadal Oscillation (PDO), were obtained from the United States' National Oceanic and Atmospheric Administration (NOAA), which can be downloaded from https://www.esrl.noaa.gov/psd/data/climateindices/list. North Atlantic SST tripole (NAT) index can be downloaded from https://cmdp.ncc-cma.net/download/ precipitation/diagnosis/NAT/NAT-mon-nor.tms.
In this study, we focused on the period from 1961 to 2017, and all the above-mentioned data covers this study period. For the seasonality analyses, SST data were averaged into seasonal values: Winter (January, February, and December of the preceding year), spring (March-May), and summer (June-August). The two temperatures (i.e., T and T e ) were averaged into summer (June-August) values.

Calculation of Equivalent Temperature
Pielke et al. [6] proposed that because a single temperature variable does not contain humidity-related heat, it cannot fully represent the near-surface heat content. For example, the change in heat characterized by an increase of 1 • C in dew-point temperature is comparable to that characterized by a direct increase of 2.5 • C in air temperature [29]. Therefore, T e , defined based on the ME of moist air, can comprehensively characterize the joint variations of near-surface temperature-humidity. The detailed definition is as follows: In Equation (1), C p represents the specific heat capacity of dry air, which is approximately equal to 1.1013 × 10 −3 MJ kg −1 • C −1 ; T represents the air temperature (unit: • C); L v represents the latent heat of vaporization, which is a function of T (calculated using the equation provided by Henderson-Sellers [30]); and q represents the specific humidity (unit: g/kg), obtained from the station air temperature, relative humidity and pressure. Equation (1) can be further transformed to calculate the T e : Equation (2) indicates that T e is composed of two parts: One part is air temperature (T) and the other part is the "temperature" component related to latent heat ( L v q C p , short for "LH term" hereafter). In particular, for the LH term, it represents the coupled nonlinear variations of T and q. Therefore, using T e can comprehensively characterize the magnitude of surface heat content (or enthalpy) when joint temperature-humidity change is considered. Physically, it represents the potential "air temperature" if all of the latent heat of moist air is isobarically converted into sensible heat. In addition, T e is similar to the heat index, humiture and apparent temperature in physical terms, and all of those indices can characterize the comprehensive variations of surface meteorological conditions (such as temperature, humidity and wind).

Ensemble Empirical Mode Decomposition
To comprehensively understand the similarities and differences of T and T e from the perspective of multiple timescales, we introduce a widely used time-frequency decomposition method for studying climate change, namely, the ensemble empirical mode decomposition (EEMD) method [31]. The EEMD method was developed from the empirical mode decomposition (EMD) method [32] to overcome the mode mixing problem of the original method by adding random white noise into the original time series in a large number of trials. The EEMD method has been widely applied in the research of meteorology, hydrology and climatology (e.g., Qian et al. [33], Franzke [34] and Kim et al. [35]). Detailed information on the EEMD can be found in the study of Wu and Huang [31]. To minimize the impact of noise, 10,000 simulations were selected, and the commonly-used number of decomposed modes was determined using the algorithm log 2 (N)−1 (detailed information can be found in Wu and Huang [36]), where N was the sample length. In our case, considering the situation of all stations, five modes (calculated by log 2 (N)-1) were finally selected, i.e., all the stations were decomposed into four (1-4) intrinsic mode functions (IMFs) and a residual. The original MATLAB code of EEMD and the corresponding introductions and tutorial can be found on https: //in.ncu.edu.tw/~{}ncu34951/research1.htm.
In addition to the aforementioned methods, some common climatological statistical methods were also used (i.e., rotated EOF analysis, correlation analysis and power spectrum analysis). For studying the multi-timescale differences between T and T e , we first performed the EEMD method on all the grids of T and T e to get their five IMFs components (i.e., IMFs 1-4 and the residual component), respectively. EOF analysis was then conducted for each IMF from the EEMD results, and the first ten modes were obtained. Moreover, according to Wilks [37], we further performed rotation to the original 10 EOF modes to enhance the physical meanings of the leading modes of each IMF. Ultimately, the first mode of the rotated EOF (REOF) result represented the most important variability of each IMF. In addition, the mean periods of REOF time series of each IMF scale were calculated using the fast Fourier transform method. Moreover, we had also calculated the proportion of variance explained of each IMF modes to the total variance of the original temperatures (i.e., T and T e ). The results show that the variations of both T and T e are predominantly controlled by the interannual variability (i.e., IMF1; Figure S1 and S2).
It should be noted that to better interpret and present the physical meanings of the REOF loadings of each timescale, we had conducted the "scaling conventions" following Wilks [37] beyond the original EOF/REOF results. For the following REOF results, the eigenvector (principal component) was multiplied (divided) by (λ) 1/2 , where λ is the eigenvalue of the corresponding REOF mode. Wilks [37] revealed that the rescaled eigenvector elements are more directly interpretable in terms of the relationship between the principal components and the original data when using this rescaling. The advantage of this rescaling procedure is that the rescaled principal component has been standardized and nondimensionalized (because (λ) 1/2 is equal to the standard deviation of the corresponding principal component), which in turn means the rescaled eigenvector has the same unit (i.e., • C) and physical meanings (i.e., temperature variation) as the original data has. Therefore, the REOF loadings in our results show not only the strength of the connection of the respective station to each EOF, but also the amplitude of the changes in related temperatures (i.e., T and T e ) when the corresponding rescaled principal component increased by 1.
In addition, the decomposed IMFs using EEMD will result in changes in the effective degree of freedom of the data. Hence, to avoid the statistical artificial significance from the above-mentioned circumstance, we used the method in reference [38] to calculate the effective degree of freedom before testing the correlation coefficient and then used the two-tailed Student's t-test to perform the significance test.

Leading Modes of T and T e at Different Timescales
To discuss the differences and similarities between the characteristics of T and T e , we first show the EEMD-REOF results for T and T e at multiple timescales in China. Generally, both T and T e can be decomposed into five frequency bands: IMFs 1-4 (3-to 4-year period, quasi-6-year period, quasi-19-year period and quasi-57-year period, respectively) (Figures 2-5) and the residual (nonlinear trend, Figure 6). First inspections show that both the spatial distribution and magnitude of T (T e ) are different among timescales; at a specific scale, the spatial distribution and magnitude between T and T e are also different. Despite these differences, there was also a high correlation between temporal variations in T and T e at certain frequencies. Since the signs of REOF spatial-temporal results are arbitrary, and hence we take one condition as an example (i.e., shown in Figures 2-6) in the following text to describe the related EEMD-REOF results. Based on the rescaling procedure (which has been introduced in the previous section), the magnitudes and signs of REOF loadings hereafter refer to temperature changes when the corresponding principal component (PC; or time series) increased by 1. It should be noted that the REOF spatial-temporal results under the exact opposite condition are also true.
Specifically, the first REOF mode of the IMF 1 of T primarily shows a high-frequency time variation at the interannual scale (period of 3.2-3.6 years; Figure 2c), and the spatial pattern exhibits consistent changes throughout China (except for parts of South China and regions to the north of 45 • N in China; Figure 2a). Evident variations are detected in the Yangtze River basin, with an increase in T of more than 0.4 • C. For the first REOF mode of IMF 1 of T e , the spatial feature of T e also exhibits a consistent distribution pattern across the country, the PC of T and T e are significantly positively correlated (the correlation coefficient is 0.90, p < 0.01). However, unlike the T distribution, where the large values are distributed solely in the Yangtze River basin, the large values of T e are located in both the Yangtze River basin and most regions in northern China. In addition, the amplitude of change in T e is significantly greater than that in T, specifically in northern China, where the amplitude of variation in T e is approximately twice than that in T (Figure 2a vs. Figure 2b). Specifically, the first REOF mode of the IMF 1 of T primarily shows a high-frequency time variation at the interannual scale (period of 3.2-3.6 years; Figure 2c), and the spatial pattern exhibits consistent changes throughout China (except for parts of South China and regions to the north of 45 °N in China; Figure 2a). Evident variations are detected in the Yangtze River basin, with an increase in T of more than 0.4 °C. For the first REOF mode of IMF 1 of Te, the spatial feature of Te also exhibits a consistent distribution pattern across the country, the PC of T and Te are significantly positively correlated (the correlation coefficient is 0.90, p < 0.01). However, unlike the T distribution, where the large values are distributed solely in the Yangtze River basin, the large values of Te are located in both the Yangtze River basin and most regions in northern China. In addition, the amplitude of change in Te is significantly greater than that in T, specifically in northern China, where the amplitude of variation in Te is approximately twice than that in T (Figure 2a vs. Figure 2b). In the first REOF mode of IMF 2 of T with a quasi-6-year period (Figure 3c), the spatial pattern shows a "negative-positive" anti-phase relationship between northern and southern China, with a boundary of approximately 35° N, i.e., the temperature variations between southern and northern China are opposites (Figure 3a). The amplitude of the temperature variation of IMF 2 in the Yangtze-Huaihe region (approximately 0.4 °C) is significantly larger than that in the North China (smaller than −0.15 °C), which suggest that IMF 2 is probably dominated by the temperature variations in the quasi-6-year period in the Yangtze-Huaihe region. On the IMF 2 scale of Te, in terms of both the spatial distribution and time series, there are considerable differences between T and Te (Figure 3a vs. Figure  3b). Specifically, although Te also exhibits opposite distributions in the south and north regions, the variations of Te in northern China (especially in North China and Northeast China) are much stronger than that of T, and the largest warming amplitude is greater than 0.6 °C. Different from the IMF 2 leading mode of T that reflects the variations in the Yangtze-Huaihe region, the IMF 2 leading mode of Te primarily reflects the significant anti-variations between the middle and low reaches of the Yangtze River (MLRYR) and northern China. In addition, there is no significant correlation between the time series corresponding to T and Te (the correlation coefficient is only 0.22).
On the decadal scale (IMF 3 of a quasi-19-year period), the first REOF mode of IMF 3 of T also exhibits an anti-phase distribution in both northern and southern China (but with the Yangtze River as the boundary; Figure 4a). In contrast to IMF 2, the IMF 3′s REOF leading mode of T generally represents the variation in the Bohai Rim region (with the largest cooling of more than −0.6 °C), and In the first REOF mode of IMF 2 of T with a quasi-6-year period (Figure 3c), the spatial pattern shows a "negative-positive" anti-phase relationship between northern and southern China, with a boundary of approximately 35 • N, i.e., the temperature variations between southern and northern China are opposites (Figure 3a). The amplitude of the temperature variation of IMF 2 in the Yangtze-Huaihe region (approximately 0.4 • C) is significantly larger than that in the North China (smaller than −0.15 • C), which suggest that IMF 2 is probably dominated by the temperature variations in the quasi-6-year period in the Yangtze-Huaihe region. On the IMF 2 scale of T e , in terms of both the spatial distribution and time series, there are considerable differences between T and T e (Figure 3a vs. Figure 3b). Specifically, although T e also exhibits opposite distributions in the south and north regions, the variations of T e in northern China (especially in North China and Northeast China) are much stronger than that of T, and the largest warming amplitude is greater than 0.6 • C. Different from the IMF 2 leading mode of T that reflects the variations in the Yangtze-Huaihe region, the IMF 2 leading mode of T e primarily reflects the significant anti-variations between the middle and low reaches of the Yangtze River (MLRYR) and northern China. In addition, there is no significant correlation between the time series corresponding to T and T e (the correlation coefficient is only 0.22).
On the decadal scale (IMF 3 of a quasi-19-year period), the first REOF mode of IMF 3 of T also exhibits an anti-phase distribution in both northern and southern China (but with the Yangtze River as the boundary; Figure 4a). In contrast to IMF 2, the IMF 3 s REOF leading mode of T generally represents the variation in the Bohai Rim region (with the largest cooling of more than −0.6 • C), and the T variations in southern China are very slight (slight warming of no more than 0.2 • C). For T e , although the correlation of the PCs corresponding to T and T e is high (0.76, p < 0.01; Figure 4c), their spatial fields are distinguishable. In contrast to the spatial distribution of T exhibiting anti-variations in southern and northern China, the spatial feature of T e exhibits cooling in most regions of China (except for slight warming occurs in some areas of MLRYR; Figure 4b), and the amplitude of T e variation is significantly larger than that of T (more than −0.4 • C; Figure 4a vs. Figure 4b).   In the quasi-57-year cycle (IMF 4), overall, the variation in T is also generally consistent nationwide (Figure 5a); the large temperature variations primarily located in North China, the eastern   In the quasi-57-year cycle (IMF 4), overall, the variation in T is also generally consistent nationwide (Figure 5a); the large temperature variations primarily located in North China, the eastern  In the quasi-57-year cycle (IMF 4), overall, the variation in T is also generally consistent nationwide ( Figure 5a); the large temperature variations primarily located in North China, the eastern part of Northwest China and East China (where the temperature increases approximately 0.3-0.4 • C). In other regions, the temperature variations of T are weak (approximately ± 0.15 • C; Figure 5a). Similar to the spatial distribution of T, T e also shows consistent warming in most regions of China (Figure 5b), and their PCs are highly positively correlated (0.95, p < 0.01; Figure 5c). However, the primary difference between them is the amplitude of variations. In detail, the variation amplitude of T e in North China, Northwest China and East China is significantly larger than that of T (Figure 5b).  Figure 5a). Similar to the spatial distribution of T, Te also shows consistent warming in most regions of China (Figure 5b), and their PCs are highly positively correlated (0.95, p < 0.01; Figure 5c). However, the primary difference between them is the amplitude of variations. In detail, the variation amplitude of Te in North China, Northwest China and East China is significantly larger than that of T (Figure 5b). For the residual component (which exhibited a nonlinear trend), the temperature variation is dominated by warming in the northern region, and the amplitude of such warming is higher than 0.5 °C in most northern regions. In the southern region, the temperature variation is relatively weak, and there is a slight cooling in the Yangtze River basin, where the temperature variation is approximately within ± 0.2 °C (Figure 6a). Furthermore, T and Te are consistent in that warming occurs in northern regions, and they are different in that the cooling is more prominent for Te than for T in some southern regions (less than −0.4 °C in most regions; Figure 6b). It is noted that the REOF spatial-temporal variations of the residual mode are consistent with the warming trend in China during the recent decades [39]. For the residual component (which exhibited a nonlinear trend), the temperature variation is dominated by warming in the northern region, and the amplitude of such warming is higher than 0.5 • C in most northern regions. In the southern region, the temperature variation is relatively weak, and there is a slight cooling in the Yangtze River basin, where the temperature variation is approximately within ± 0.2 • C ( Figure 6a). Furthermore, T and T e are consistent in that warming occurs in northern regions, and they are different in that the cooling is more prominent for T e than for T in some southern regions (less than −0.4 • C in most regions; Figure 6b). It is noted that the REOF spatial-temporal variations of the residual mode are consistent with the warming trend in China during the recent decades [39].
In addition, based on the percent variances of all the leading mode of T and T e at multiple timescales, REOF leading mode of both T and T e at interannual scale (IMF 1), multidecadal scale (IMF 4) and especially the nonlinear trend scale (the residual component), can all well represent the actual temperature variation of the corresponding scale. In conclusion, similarities of T and T e in both spatial and temporal variations are detected in three leading modes at IMF 1, IMF 4 and residual component scales, while differences between two temperatures (i.e., T and T e ) are detected in the other two scales (i.e., IMF 2 and IMF 3). It should be noted that although we have used a station redistribution algorithm to unify the station density in each grid, the inhomogeneity of the spatial distribution in still remains, which might introduce some uncertainties to the current results. For enhancing the confidence of the REOF results, we have tested the sensitivity of the EOF/REOF results to the station selection from the perspective of bootstrap resampling, details can be found in the supporting information (Text S1 and Figure S3). The results of random sampling show that the REOF1 spatial-temporal patterns at each timescale shown in our manuscript are generally not sensitive to the number and distribution of the stations. In addition, based on the percent variances of all the leading mode of T and Te at multiple timescales, REOF leading mode of both T and Te at interannual scale (IMF 1), multidecadal scale (IMF 4) and especially the nonlinear trend scale (the residual component), can all well represent the actual temperature variation of the corresponding scale. In conclusion, similarities of T and Te in both spatial and temporal variations are detected in three leading modes at IMF 1, IMF 4 and residual component scales, while differences between two temperatures (i.e., T and Te) are detected in the other two scales (i.e., IMF 2 and IMF 3). It should be noted that although we have used a station redistribution algorithm to unify the station density in each grid, the inhomogeneity of the spatial distribution in still remains, which might introduce some uncertainties to the current results. For enhancing the confidence of the REOF results, we have tested the sensitivity of the EOF/REOF results to the station selection from the perspective of bootstrap resampling, details can be found in the supporting information (Text S1 and Figure S3). The results of random sampling show that the REOF1 spatialtemporal patterns at each timescale shown in our manuscript are generally not sensitive to the number and distribution of the stations.

The Role of T and LH Terms in Determining the Variations of Te
In the previous section, we have comprehensively demonstrated the major differences between T and Te at different timescales. The differences between these two temperatures theoretically suggest the role of introducing a latent heat term on the variations of Te. Hence, in this section, further exploration and discussion will be conducted to show the relative role of T and LH terms in determining the variations of Te.
Generally, according to the formula of Te, firstly, if the variations of T and Te are similar, there are two possibilities: (1) The variability of T (VT) is far greater than that of LH term (VLH), i.e., VT >> VLH, which makes Te dominating by T; and (2) VT and VLH are comparable or VLH >> VT, and the variations of T and LH term must be similar at the same time (i.e., the correlation between T and LH term are positive and high), which would result in the enhancement of Te by the superposition effect of T and LH. Secondly, if the variations of T and Te are distinguishable, then VT and VLH are comparable or VLH >> VT, and the variations between T and LH term should be negatively or lowly correlated. Thus, Te would be dominated by the variation of LH and/or the joint variation of T and

The Role of T and LH Terms in Determining the Variations of T e
In the previous section, we have comprehensively demonstrated the major differences between T and T e at different timescales. The differences between these two temperatures theoretically suggest the role of introducing a latent heat term on the variations of T e . Hence, in this section, further exploration and discussion will be conducted to show the relative role of T and LH terms in determining the variations of T e .
Generally, according to the formula of T e , firstly, if the variations of T and T e are similar, there are two possibilities: (1) The variability of T (V T ) is far greater than that of LH term (V LH ), i.e., V T >> V LH , which makes T e dominating by T; and (2) V T and V LH are comparable or V LH >> V T , and the variations of T and LH term must be similar at the same time (i.e., the correlation between T and LH term are positive and high), which would result in the enhancement of T e by the superposition effect of T and LH. Secondly, if the variations of T and T e are distinguishable, then V T and V LH are comparable or V LH >> V T , and the variations between T and LH term should be negatively or lowly correlated. Thus, T e would be dominated by the variation of LH and/or the joint variation of T and LH. In order to verify the above-mentioned conjectures and further our understanding of the possible causes of similarities and differences between T and T e , we first additionally calculated the EEMD results of LH term (L v q/C p ) to get the five EEMD modes of each grid. Furthermore, we calculate the correlation coefficients between EEMD-decomposed T and LH terms and the ratio of the standard deviations (STDR) between those two terms at different timescales. The correlation and STDR results are shown in Figure 7. Here, we define that STDR < 0.5, 0.5 ≤ STDR ≤ 2 and STDR > 2 represent that V LH dominated the variability of T e (V Te ), V T and V LH are comparable and T dominated V Te , respectively. causes of similarities and differences between T and Te, we first additionally calculated the EEMD results of LH term (Lvq/Cp) to get the five EEMD modes of each grid. Furthermore, we calculate the correlation coefficients between EEMD-decomposed T and LH terms and the ratio of the standard deviations (STDR) between those two terms at different timescales. The correlation and STDR results are shown in Figure 7. Here, we define that STDR < 0.5, 0.5 ≤ STDR ≤ 2 and STDR > 2 represent that VLH dominated the variability of Te (VTe), VT and VLH are comparable and T dominated VTe, respectively. In detail, for IMF 1, the correlations between T and LH terms exhibit high values (more than 0.6) over the Yangtze river basin ((a1) in Figure 7a), meanwhile, STDR between T and LH terms are comparable over the same region ((b1) in Figure 7b). However, the correlations of T and LH terms are weak in northern China (no more than ±0.2), and the STDR results in northern China indicate that In detail, for IMF 1, the correlations between T and LH terms exhibit high values (more than 0.6) over the Yangtze river basin ((a1) in Figure 7a), meanwhile, STDR between T and LH terms are comparable over the same region ((b1) in Figure 7b). However, the correlations of T and LH terms are weak in northern China (no more than ±0.2), and the STDR results in northern China indicate that V LH is much larger than V T . Those pieces of evidence suggest that the warming/cooling of T e in northern China are dominated by V LH , while T and LH probably increase or decrease simultaneously in the Yangtze River basin, which finally leads to a warm-humid or a cold-dry state in the summer climate of Yangtze River basin. For IMF 2, both the spatial and temporal results of REOF 1 are distinguishable between T and T e (Figure 3). Unlike IMF 1, there is no evident correlation between V T and V LH (no more than ±0.4 over most of China; (a2) in Figure 7a), and the STDR results show that the anti-variation pattern in the REOF leading mode of T e is probably influenced by the joint V T -V LH ((b2) in Figure 7b). Those results highlight that the variation of T e is not a simple composition of T and LH, but can represent a unique complex nonlinear variation of coupled temperature-humidity. At the IMF 3 scale, the correlation pattern between T and LH resemble that of IMF 2, which also show sporadic correlations nationwide ((a3) in Figure 7a). From the STDR results of IMF 3, it is noted that V LH is much greater than V T over the regions with the strong variations of T e (Figure 4b and (b3) in Figure 7b), which indicates that V LH at the decadal scale is more important for surface heat content than V T . For the multi-decadal scale (IMF 4), the REOF 1 spatial and temporal patterns between T and T e are very similar, despite T e with larger variations. The correlations between T and LH terms exhibit high values (more than 0.6) for regions with strong variations in both T and T e (i.e., north of Yangtze River basin), which suggests that the multi-decadal variabilities of T and LH are similar ((a4) in Figure 7a). It is worth noting that STDR between T and LH terms in the previous three modes display either V LH -dominated V Te or a comparable situation, while at the multi-decadal timescale, V T seems to be stronger than previous modes, specifically for some regions, dominant features of V T are detected (i.e., some regions over MLRYR and northern China; (b4) in Figure 7b). Moreover, for the residual component, there are also high correlations (more than 0.6) all over the country, despite some evident negative correlations in North China ((a5) in Figure 7a). Interestingly, the role of V T in determining V Te further grows at the long-term trend scale from the STDR results. In detail, the STDR of T and LH show that V T and V LH are comparable over most of China, while V T is much larger than V LH in northern China, indicating that the nonlinear warming trend of T e over northern China is dominated by V T ((b5) in Figure 7b). In addition, from the REOF 1, correlations and STDR results of IMF 4 and the residual component, the relationships between T and LH on these two scales generally resemble that for IMF 1, i.e., the superposition effect of highly correlated T and LH (i.e., warm-humid or cold-dry) dominates the variation of T e .
In addition, the above-mentioned results are completely consistent with our conjectures about the role of V T and V LH in determining the similarities and differences between T and T e . In conclusion, V LH plays an equivalent or more important role than V T on V Te at interannual to decadal scales (i.e., IMF 1-IMF 3), while V T is being pronounced at multi-decadal and the long-term scales (i.e., IMF 4 and the residual component). These results also imply that the influencing factors and physical mechanisms leading to the changes in these two elements are likely to be different among timescales.

Similarities and Differences in the Linkage between SST and the Leading Mode of T and T e
Previous studies have demonstrated that there is a close relationship between temperature variation in summer in China and SST variation in many key regions around the world (Pacific, North Atlantic and Indian Oceans) [24,25,40], especially on interdecadal and multidecadal scales (e.g., Wang et al. [41]). In the previous section, we have determined that except for the residual time series (which exhibits a nonlinear trend), the four main IMFs of T and T e are interannual-to-multidecadal oscillations, and their similarities and differences are distinct at different timescales. Hence, what are the potential factors that cause a high correlation between the two temperatures (i.e., T and T e ) at the interannual, interdecadal and multidecadal scales (IMFs 1, 3 and 4)? What are the possible reasons for the low correlation between the two temperatures (i.e., T and T e ) at the interannual-to-interdecadal transition scale (IMF 2)? In this section, the possible linkage between the multi-timescale changes of the two temperatures (i.e., T and T e ) and the global SST is discussed. Note that we mainly focus on the first four IMF modes of T and T e in this section; therefore, the detrending process is firstly done on the SST before calculating its correlation with the two temperatures. Figure 8 presents the distribution of the correlation coefficients between the PCs of the first REOF modes of T and T e and the SST in summer on each timescale, meanwhile, Table 1 lists the correlation coefficients between the PCs and multiple SST indices. In general, for the IMF 1 (Figure 8a,b) and IMF 4 scales (Figure 8g,h), the corresponding SST patterns that are significantly correlated with T resemble that with T e , indicating that for these scales, the key SST patterns are associated with both V T and V LH . This also implies that the changes in SST in these key regions may not only affect the temperature variation in summer in China but also cause the coupled T-LH to increase or decrease simultaneously, i.e., warming (cooling) corresponds to high humidity (low humidity), which further causes the spatial distribution of T e to be similar to that of T. However, in IMF 2 (Figure 7c,d) and IMF 3 (Figure 7e,f) modes, the corresponding SST patterns are different, indicating that the dominant key SST regions affecting T and T e as well as the possible mechanisms are distinguishable for these two scales. This further suggests that for these two timescales, the differences in the corresponding SST patterns between T and T e potentially highlight the possible relationship between SST and V LH .
Atmosphere 2019, 10, x FOR PEER REVIEW 12 of 17 Figure 8 presents the distribution of the correlation coefficients between the PCs of the first REOF modes of T and Te and the SST in summer on each timescale, meanwhile, Table 1 lists the correlation coefficients between the PCs and multiple SST indices. In general, for the IMF 1 (Figure 8a,b) and IMF 4 scales (Figure 8g,h), the corresponding SST patterns that are significantly correlated with T resemble that with Te, indicating that for these scales, the key SST patterns are associated with both VT and VLH. This also implies that the changes in SST in these key regions may not only affect the temperature variation in summer in China but also cause the coupled T-LH to increase or decrease simultaneously, i.e., warming (cooling) corresponds to high humidity (low humidity), which further causes the spatial distribution of Te to be similar to that of T. However, in IMF 2 (Figure 7c,d) and IMF 3 (Figure 7e,f) modes, the corresponding SST patterns are different, indicating that the dominant key SST regions affecting T and Te as well as the possible mechanisms are distinguishable for these two scales. This further suggests that for these two timescales, the differences in the corresponding SST patterns between T and Te potentially highlight the possible relationship between SST and VLH. The SST fields associated with T and Te of IMF 1 all exhibit significant reverse-phase changes in eastern and the western Pacific (i.e., positive and negative correlations over the western and eastern Pacific, respectively), which is obviously the typical El Niño-Southern Oscillation (ENSO)-type pattern (Figure 8a, 8b). The correlation coefficient between the PCs of T and Te and the Niño 3.4 index is −0.30 (p < 0.05) and −0.43 (p < 0.05), respectively (Table 1), indicating that when the central and eastern equatorial Pacific SST is abnormally low in summer, T and Te is high over the whole country. In addition, T and Te are not only closely associated with ENSO-like SST anomalies (SSTA) in summer, but also positively correlated with the central and eastern equatorial Pacific SST in the preceding winter ( Figure S4). This implies that the variations of T and Te in summer may be The SST fields associated with T and T e of IMF 1 all exhibit significant reverse-phase changes in eastern and the western Pacific (i.e., positive and negative correlations over the western and eastern Pacific, respectively), which is obviously the typical El Niño-Southern Oscillation (ENSO)-type pattern (Figure 8a, 8b). The correlation coefficient between the PCs of T and T e and the Niño 3.4 index is −0.30 (p < 0.05) and −0.43 (p < 0.05), respectively (Table 1), indicating that when the central and eastern equatorial Pacific SST is abnormally low in summer, T and T e is high over the whole country. In addition, T and T e are not only closely associated with ENSO-like SST anomalies (SSTA) in summer, but also positively correlated with the central and eastern equatorial Pacific SST in the preceding winter ( Figure S4). This implies that the variations of T and T e in summer may be modulated by the phase transition from El Nño (La Niña) in the preceding winter to La Niña (El Nño) in summer. Numerous studies have revealed that such transition, e.g., from El Nño to La Niña, typically results in a significant strengthening of the western Pacific subtropical high in summer, leading to clear skies (i.e., fewer clouds) with high temperatures in China (e.g., Wang et al. [42]). Moreover, it is noted that although the correlation coefficient patterns of T and T e with SST is relatively consistent, T e is detected to have a higher correlation with more significant signals than T. This indicates that the central and eastern equatorial Pacific SST not only affects the temperature but also causes the humidity to change at the same time, leading to the warm and wet (or cold and dry) climate in summer in China. In IMF 4, the SST field is positively correlated with T and T e in the whole region of the North Atlantic (especially in the eastern part with significant positive correlations; Figure 8g,h). This correlation pattern is similar to the AMO. The correlation coefficient between the AMO index and T and T e reached 0.67 (p < 0.05) and 0.59 (p < 0.05), respectively (Table 1), indicating that the variations of T and T e are probably related to AMO in this mode. In addition, significant but sporadic correlations of T and T e with SST are also detected in some regions of Indian and Pacific Oceans (Figure 8g,h). On the remaining two scales that T and T e are evidently different (see Figures 3 and 4), their correlations with SST are also distinguishable (Figure 8c-f). In detail, at the IMF 2 scale, the SSTA that are significantly correlated with T manifest only as the low-latitude and mid-latitude dipole-type distributions in the North Atlantic (Figure 8c). But for T e , it is shown that broadly significant positive correlations are detected in the central and eastern equatorial Pacific, the northeast Pacific and Indian Oceans, whereas the correlations over North Atlantic show an insignificant NAT pattern. The correlations of T and T e with the SST indices are approximately consistent with the aforementioned spatial distribution (Table 1). These results also imply that although the relationship between T and the equatorial Pacific, North Pacific and Indian Oceans' SSTs is relatively weak in the quasi-6-year period, the abnormal change in SST in these key regions may be an important factor for the variations of surface heat content. In addition, based on the SST signals from the preceding winter and spring ( Figure S5), T is shown to correspond to the La Niña (El Niño) event, which has a weak development from the preceding winter to spring and decays to a non-significant signal during summer, whereas T e corresponds to the La Niña (El Niño) event that continuously developed from the preceding winter into the current summer. This indicates that whether the development of ENSO events can be maintained into the current summer is a possible cause of the differences in T and T e . From another perspective, continuously intensifying La Niña (El Niño) events may significantly affect the variation of humidity (precipitation) in summer in China, which in turn reflects in the variations of surface heat content. On the IMF 3 scale (Figure 8e,f), the T-related SST field does not feature a significantly correlated signal; however, the T e -related SST field has negative correlations in the central, east and north equatorial Pacific Ocean, indicating that this type of SST is more closely related with V LH . Furthermore, although the PCs of T and T e have a high correlation, the corresponding SST patterns are inconsistent. This may be because the correlation between the PCs of T and T e is relatively low before the 1990s (0.60 in the period 1961-1985, compared with 0.90 in the period 1986-2010).
Except for the above four modes, the PCs of the residual component of T and T e are consistent with the trend of global warming (Figure 6c), and both correspond to the global-scale increase of SST ( Figure S6). The change of the two temperatures is probably due to increased CO 2 concentrations caused by human activities [2].

Conclusions and Discussion
Based on the EEMD method and the observational data (i.e., surface air temperature, relative humidity, pressure and global SST), this study analyzed the similarities and differences of the surface air temperature (T) and equivalent temperature (T e ) defined by the ME in summer at multiple timescales for the period 1961-2017, and their relationship with global SST. The conclusions are as follows: (1) The similarity between T and T e is that both can be decomposed into five timescales: IMFs 1-4 (3.2-3.6-year period, quasi-6-year period, quasi-19-year period and quasi-57-year period, respectively) and a residual component (nonlinear trend). In terms of the results of all timescales, the variation of the leading modes (i.e., the first REOF modes) of T and T e is dominated by the high-frequency interannual (3.2-3.6-year) and multidecadal (quasi-57-year) variations and the nonlinear warming trend. Moreover, except for the interannual-to-interdecadal transitional scale (quasi-6-year), there are significant positive correlations between the PCs of the first REOF modes of T and T e .
(2) The differences between T and T e can be mainly concluded as follows: Seen from all the REOF leading modes of T and T e at multiple timescales, the magnitude of T e variation is always greater than that of T, indicating that the variation of multi-variable-defined surface heat content is larger than that of single temperature-defined. For IMF 1, IMF 4 and the residual component scales, the REOF 1 spatial patterns are generally similar, except for the large discrepancies in the amplitude of between T and T e variations. Further explorations suggest that V T and V LH are generally positively correlated and both changes of them tend to enhance the variability of surface heat content simultaneously, and in turn, leading to a warm-humid (or cold-dry) state in summer. However, at IMF 2 and IMF 3 scales, the differences between T and T e are evident, not only in the amplitude of variability but in the spatial pattern as well. For IMF 2 scale, although V T and V LH are comparable, their correlations are relatively weak, which in turn exert a unique coupled effect on the surface heat content. For IMF 3 scale, V LH are much larger than V T , and there are low correlations between them, which ultimately make V LH dominate the variation of the surface heat content (T e ).
(3) From the correlations of SST with T and T e , similarities are also found in IMF 1, IMF 4 and the residual component scales, which both T and T e show the ENSO phase transitions (e.g., El Niño phase to La Niña phase; IMF 1), AMO-like (IMF 4) and the global warming-like SSTA patterns (the residual component), respectively. In terms of quasi-6-year (IMF 2) and quasi-19-year periods (IMF 3), the SST fields associated with T and T e are different: The ENSO events that developed during the preceding winter and persisted until summer for the quasi-6-year period possibly caused the difference between T and T e . For the quasi-19-year period, the SSTA in the key regions of the central, east, and the north equatorial Pacific Ocean was closely related to T e but not significantly correlated with T.
From the current results, this paper primarily emphasizes the similarities and differences of T and T e at multiple timescales over China and shows some linkage between SST and T/T e . Additionally, there are some limitations in this study, primarily in the following aspects. First, in terms of SST, this article only explored its relationship with T and T e from the perspective of correlation, and the corresponding physical processes and mechanisms require further study. Second, in addition to the SST, there are many factors affecting T and T e . For example, the results of Findell et al. [13] revealed that the land-use change in the historical period made parts of northern China significantly warmer and drier and made some parts of southern China warmer and more humid. The effects of other factors on T and T e at different timescales require comprehensive investigation. Third, based on the main results of this paper, for the interannual-to-interdecadal scale (6 to 19-year period), there are very robust differences between T and T e . The specific physical mechanism for the difference also deserves comprehensive discussion. In addition, we only show the results of the differences between T and T e over China, however, the differences between T and T e may be a global phenomenon (e.g., Li et al. [43]), since the reason for their differences is possibly due to the large-scale external forcings (e.g., anthropogenic greenhouse gases) or the internal variability (e.g., Rossby wave trains). Many other studies have also found that there are significant differences between T and T e over other countries or globally (e.g., Pielke et al. [6], Younger et al. [17] and Li et al. [43]), and their differences at a multiscale over the rest of the world, as well as whether there is a potential linkage with our findings will require further investigations.
Despite these shortcomings, this article comprehensively discusses the characteristics of variations of T e at multiple timescales in recent decades and their similarities and differences with those of T from a new perspective of joint temperature-humidity-defined surface heat content. The relevant results are intended to provide clues for the variations in the T e -related total surface heat content and provide insights and bases for the comprehensive understanding of climate change as well as the related compound extreme events (e.g., compound temperature-humidity extremes) in recent decades.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/10/8/447/s1, Figure S1: The dominant IMF component for T and T e from the EEMD results, Figure S2: The percent contribution of the variance of each IMF to the total variance, Figure S3 and Text S1: The sensitivity of the number and distribution of stations to the REOF results, Figure S4: The spatial distributions of correlation coefficients between preceding winter SST and PCs of the first REOF mode of IMF 1−4, Figure S5: The spatial distributions of correlation coefficients between spring SST and PCs of the first REOF mode of IMF 1−4, Figure S6: The spatial distributions of correlation coefficients between summer SST and PCs of the first REOF mode at the residual scale.