Variability of the Early Summer Temperature in the Southeastern Tibetan Plateau in Recent Centuries and the Linkage to the Indian Ocean Basin Mode

As one of the areas’ most sensitive to climate change, the Tibetan Plateau (TP) not only has a significant climate change signal but the TP also significantly affects climate change in eastern China and globally [1–5]. The southeastern TP (SETP) contains the sources of the Yangtze, Yellow, and Lancang Rivers and is therefore referred to as the “water tower of China” [6, 7]. It is the most important water source and one of the largest and highest natural wetlands and biodiversity areas in China [8, 9]. However, owing to the limited length of observed climate records, knowledge regarding the characteristic variations in climate change in this area may be incomplete [10]. To understand long-term climate change characteristics, high-resolution climate proxy data, such as tree rings, ice cores, and lake sediments, have been used to reconstruct past climatic changes on the TP [11–16]. On the TP, tree-ring data represent the most commonly used data for the reconstruction of climate change during the past centuries, even millennia. To date, dozens of tree chronologies over the millennium have been obtained in this region [17–20]. Among these, Yang et al. [21] represented the longest chronology in Asia, which was approximately 6700 years. Because of the wide and long-lived natural forests (such as forests comprising Picea balfouriana and Juniperus tibetica) on the SETP, recently many dendroclimatologists have performed paleoclimate research using multiple tree-ring parameters. The most commonly reconstructed parameters include temperature [22–24], drought [25–27], streamflow [28–30], and precipitation [31–33] over periods of several centuries to millennia. These reconstructions help us improve our understanding of past climate change over the SETP. Recently, with the accumulation of tree chronologies, dendroclimatic studies have focused on large-scale or regional climate field reconstruction [25, 34–39]. However, most studies in the SETP used one or more sample sites to reconstruct the regional mean climate variation. To date, the only STEP climate field reconstructions include the Monsoon Asia Drought Atlas (MADA) [35], standardized precipitation evapotranspiration index (SPEI) for the eastern TP [25], and summer temperature anomalies for East Asia [36]. Cook et al. [35] and Deng et al. [25] focused on drought reconstruction, while Cook et al. [36] principally discussed large-area temperature changes, with few details and a lack of skill for the SETP [25]. Therefore, our objective was to reconstruct the SETP temperature field and elucidate its variations in detail.

The objectives of the present study were to (1) develop a tree-ring network that includes 53 chronologies (one δ¹³C chronology, one regional maximum latewood density (MXD) chronology, and 51 width chronologies); (2) reconstruct the May-June temperature field for the SETP using the modified point-by-point reconstruction (PPR) method with the discussion of the resultant spatial and temporal characteristics; and (3) elucidate the possible forces driving temperature variability in the SETP.

The study area (SETP) is located between 25–35° N and 90–104° E. This region is one of the most important areas for dendroclimatic research. The tree-ring network included 53 chronologies, the main data sources of which were the International Tree Ring Data Bank (ITRDB, https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring), our contributed chronologies, and several previously published chronologies (e.g., Yin et al. [23]; Xu et al. [40]; Li et al. [41, 42]; Ye et al. [43]). Thirty-seven chronologies were obtained from the ITRDB, which has been used to depict temperature changes in Asia over the last 2000 years [39]. The MXD chronology is a regional chronology developed from four sample sites [23]. The locations of these chronologies are shown in Figure 1.

For the chronologies developed by our group, the raw widths were measured using a Lintab measuring system with a precision of 0.01 mm. Subsequently, the COFFCHA program was used to ensure accurate dating results [44]. Nonclimatic information was removed or reduced from the original tree-ring width data, while retained climatic information was maximized using the ARSTAN_XP software [45] (http://www.ldeo.columbia.edu/tree-ring-laboratory/resources/software). Most tree-ring series were detrended by a spline function using a time step of 100 y. Three different types of chronologies were obtained, including a standard chronology (STD), residual chronology (RES), and ARSTAN chronology (ARS). To retain more of the characteristics of low-frequency variations, the STD chronology was used to produce the reconstructions in this study. Because the earlier years of tree-ring chronologies were not reliable [25], the value of the expressed population signal (EPS) was used to determine the reliable years; an EPS >0.85 was used to determine the reliability threshold.

A tree-ring network with 53 chronologies (one δ¹³C chronology, one regional maximum latewood density (MXD) chronology, and fifty-one width chronologies) was developed. All chronologies were used to represent the temperature [23, 39–43].

Climate data from 28 meteorological stations in the SETP were used in this study (Figure 1 and Table 1). Climate data were collected from the Chinese Homogeneous Temperature Data Set (CHTDS), which was detected and adjusted for temporal inhomogeneity [46, 47]. Quality control was performed by the National Meteorological Information Center of the China Meteorological Administration. The dataset was downloaded from http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_TEM_MUT_HOMO.html. The climate data covered the period of 1961–2004. CRU TS v. 3.20 data with a spatial resolution of $0.5°×0.5°$ (https://crudata.uea.ac.uk/cru/data/hrg/), supported by the Climatic Research Unit (CRU) of the University of East Anglia, were used to analyze the spatial correlation between the recorded series and reconstructed series.

The point-by-point reconstruction (PPR) method [34] was used to reconstruct the temperature field, which has been successfully used to reconstruct past climate fields in some regions [25, 34–37, 48]. The PPR method was described in detail by Cook et al. [34].

In general, a tree-ring network is developed in four steps [34–37]. In the first step, several chronologies from the candidate tree-ring chronologies were selected as “a pool,” which are likely to be closely related to the temperature at a given point or station. Generally, a “search radius,” depending on the correlation decay distance (CDD, e.g., 150 km), is selected to locate candidate tree-ring chronologies [34–36]; however, it may be less reliable in regions with complex topography [37]. In this study, a “modified search radius” (i.e., search spatial correlation) was used to locate candidate tree-ring chronologies [37]. Instead of the search radius, the spatial correlation was calculated between the instrument temperature records and the surrounding gridded temperature records from the CRU TS v. 3.20 dataset. Similar to the search radius method [34], a correlation coefficient of 0.8 was selected as the initial search spatial correlation. If there were insufficient tree-ring chronologies (at least 20 chronologies) within the area, the search area was expanded by reducing the spatial correlation coefficient by intervals of 0.05 until at least 20 chronologies were available. This method is more appropriate for the SETP owing to its variable climate and complex topography. As most chronologies provided warm-season temperatures [23, 39–43], in this study, the target temperature for reconstruction was the mean temperature from May to June. Second, we calculated the correlation between the temperature at the target meteorological stations and the candidate tree-ring chronologies, as the candidate tree-ring chronologies may not be significantly correlated with the temperature at a meteorological station. If the correlation coefficient exceeded 0.32 (⁠$P<0.05$⁠), it was retained for the next stage of analysis (Table 1). Third, to reduce the dimensionality of the field and retain the principal component signal of the retained chronologies, principal component regression analysis (PCRA) [49, 50] was used to construct a regression model. Principal components with eigenvalues greater than unity were extracted. Finally, the stability and reliability of the reconstructed temperature field were verified. Therefore, the correlation coefficients, explained variance, reduction of error (RE), and root mean square error (RMSE) were used to verify the reliability of each station. The RE is a highly sensitive signal, ranging from -∞ to 1, and a value greater than 0 indicates that the model is acceptable. Principal component analysis (PCA) and rotated principal component analysis (RPCA) were then used to spatially evaluate the reconstruction [27, 32].

Table S1 shows the mean correlation coefficients and the mean number of chronologies with a significant correlation (⁠$P<0.05$⁠) between the observed temperature and chronologies at every meteorological station in the SETP. These results are similar to those of previous studies [23, 39–43]. The current results suggest that 15 and 12 chronologies were significantly correlated with the mean correlation coefficients being 0.43 and 0.41 in May and June, respectively. These were the two were maximum observed during the year. This indicated that the temperature in early summer (May and June) was the most common climatic factor, which was represented by tree growth in the SETP. Therefore, early summer (May and June) temperatures were selected as the target of the reconstruction in this study.

The number of candidate tree-ring chronologies used in the regression model differed for each target station (Table 1). Generally, most of the search spatial correlation coefficients were not less than 0.5 (excluding those for the LaS and XiangGL stations, with search spatial correlation coefficients of 0.4 and 0.35, respectively) when sufficient candidate tree-ring chronologies were used. If the correlation coefficient between the candidate tree-ring chronology and temperature at the target meteorological station was not less than 0.32 (⁠$P<0.05$⁠), the chronology was retained for the PCRA analysis. Three to fifteen chronologies were generally retained (Table 1); however, two and one chronologies were retained for the XiaoJ and XiangGL stations, respectively. PCR was used to conduct the regression model for the XiaoJ station, while linear regression was performed instead of PCR for the XiangGL station. Generally, this work reconstructed the SETP temperature field with a common period of 1730–1998 and the longest single period of 1480–2002 (Table 1 and Figure 2).

The calibration and verification results are presented in Table 1. The RMSE ranged from 0.45 to 0.91, and all of them were less than the standard deviation accordingly, while the correlation coefficients ranged from 0.43 (DaoF and DaoC stations) to 0.73 (GanZ station), and all satisfied the 0.05 significance test. The explained variance ranged from 18.1% (DaoF and DaoC stations) to 53.5% (GanZ station). Although the RE values all exceeded 0, the RE values at stations LaS, DaoF, XinL, and DaoC were relatively low (0.01, 0.08, 0.07, and 0.09, respectively). This may suggest that although the reconstructions were acceptable, the resultant accuracies were poor and should be carefully considered.

Spatial features have received more attention for field reconstruction. First, the field correlation coefficients were calculated between the observed and reconstructed temperature fields during the common period of 1961–1998. As shown in Figure S1, the correlation coefficients varied from 0.97 to 0.99 with little fluctuation and high significance (0.0001). This indicated that the reconstructed temperature spatial pattern was similar to that observed for a given year and the reconstruction was reasonable for the annual spatial distribution. The correlation coefficient between the average observed and reconstructed values was 0.99 (⁠$p=0.0001$⁠). PCA and RPCA were performed to verify the consistency of the spatial patterns between the observed and reconstructed data. The observed data were from 1961 to 2004, and the reconstructed data were for the common period 1730–1998. PCA and RPCA are powerful tools for detecting spatial patterns in regional climate [34, 37]. The PCA results are presented in Figure S2 and Table 2. The total contribution rates of the first three principal components were 86.42% and 76.06% for the observed and reconstructed data, respectively (Table 2). According to the observed temperature field PCA (Figure S2), the eigenvector values of the first principal component in the region were all positive, indicating that the trends of the temperature evolution in the area were similar. The eigenvector of the second (third) principal component which exhibited a northwest-southeast (northeast-southwest) anti-phase, indicating that when the temperature in the northwest (northeast) is higher, the temperature in the southwest (southeast) is lower, and vice versa. The variance contributions of the second and third principal components were 8.97% and 4.92%, respectively. For the reconstructed temperature field, the variance contributions of the first three principal components were 42.99%, 20.46%, and 12.61%, respectively, and the eigenvector field exhibited spatial patterns similar to those observed. For RPCA (Figure S3, Table 2), the results of the observed temperature field showed that the highest loading (>0.45) of the first rotation factor (VF1) was mainly located northwest of the SETP. This may be referred to as the “NW pattern.” Similarly, VF2 refers to the southeast factor, called the “SE pattern,” with the highest loading in the southeast SETP. The highest loading of the VF3 was observed in the northeastern SETP and was named the “NE pattern.” The first three reconstructed RPCA field patterns exhibited the “NW pattern,” “SE pattern,” and “NE pattern,” respectively, which corresponded to the instrumental data. In general, the reconstructed temperature field reflects spatial distribution characteristics.

The ten-year-average average standardized temperatures for each meteorological site are shown in Figure 3(a). The figure exhibits some important climate events, such as recent warming and some low-temperature periods during the little ice age between the 15^th and early 19^th centuries. The main low-temperature periods were during the 1510s–1560s, 1610s–1640s, 1700s–1720s, and 1750s–1800s. The period indicates the entire decade of the year indicated (for example, 1510s refers to the average from A.D. 1510 to A.D. 1519). Figures 3(b)–3(d) showed three typical spatial distributions of the temperature departure in the SETP. For Figure 3(b), the mean temperature departure during the period 1735 to 1751 suggested the northwest-southeast anti-phase, and the temperature was positive in the northwest SETP region, while it was negative in the southeast SETP region. This characterized the second principal component (Figure S2). For the typical periods 1752–1787(Figure 3(c)) and 1986–1998 (Figure 3(d)), all characterized the first principal component (Figure S2), with similar temperature departures (positive or negative) in the whole area. In addition, according to the PCA results (Figure S2), the temperature changes in the SETP were similar. We calculated the regional average temperature during the common period of 1730–1998 (Figure 4), and the correlation coefficients between the reconstructed regional average temperature and each observed series varied from 0.44 and 0.83 (most of which were above 0.6 and all significance tests at level 0.01), while the correlation coefficient between the reconstructed and observed regional temperature was 0.78, with the 0.001 confidence degree reaching a very significant level. The average SETP temperature series suggested that the relatively warm periods occurred during AD 1730–1754, 1814–1829, 1853–1886, 1910–1932, 1955–1966, and 1981–1998, while the relatively cool periods occurred during 1755–1788, 1799–1813, 1830–1852, 1887–1909, 1933–1954, and 1967–1980.

To evaluate the spatiotemporal representativeness of the regional average temperature, spatial correlations were calculated between the reconstructed and observed temperatures and those derived from the CRU gridded dataset (TS 4.03) for the Northern Hemisphere during the period 1961–1998. The spatial correlations exhibited similar patterns for the observed (Figure 5(a)) and reconstructed (Figure 5(b)) temperatures. Both the reconstructed and observed temperatures were significantly correlated with CRU-gridded surface temperatures over a large geographic area. This suggests that the temperature reconstruction accurately captured the regional temperature and climatic variations over a large SETP area. The results of the spatial correlation analysis also confirmed that the SETP temperature was mainly synchronized with the temperature changes in the middle and low latitudes of the Northern Hemisphere (0–30°N). In contrast, the synchronization with the high latitudes was comparatively insignificant, which was similar to the previous research [36]. However, it was still highly synchronous with the temperature changes in Alaska. The oceanic and atmospheric oscillations, such as the El Nĩno-Southern oscillation (ENSO) and Pacific decadal oscillation (PDO), were the most important diver forcing for the large-scale co-varying climate changes [51, 52]. Not only instrumental records [53, 54], but also research based on proxy data [55–57] showed that there were hydroclimate or temperature teleconnections between Asia and North America on interannual to multidecadal timescales. Our results also suggest synchronous temperature changes between the SETP and the Alaska. The Asian Pacific Oscillation (APO) and Pacific decadal oscillation (PDO) may have a key role in large-scale teleconnection. This may occur via atmospheric circulations such as the monsoon, western Pacific subtropical high pressure, and westerlies, to link these co-varying climate changes [54–57].

Multiscale climate signals from the original temperature reconstruction were extracted using the ensemble empirical mode decomposition (EEMD) method (Figure 6). EEMD decomposed the original data into six main intrinsic mode functions (IMF1–IMF6) and one residual term to present the periodic climate signals at different time scales in the reconstructed temperature series. The IMF1–IMF6 were highly correlated with the original temperature data (⁠$P<0.05$⁠), and the variance contribution rates were 45.51%, 13.50%, 8.69%, 12.55%, 3.46%, and 5.33%, respectively (Table 3). The high correlations and variance contributions indicate that these IMFs are well-suited for extracting high- and low-frequency signals from the original data. The correlation coefficients of IMF1 and IMF2 with the original dataset were 0.646 and 0.460, respectively, and their variance contribution rates were the largest, indicating that IMF1 and IMF2 may better represent the cycle characteristics of the original reconstruction than the other components.

The six intrinsic mode functions showed periodic variation characteristics in the reconstructed series at different timescales, and every IMF may be related to certain physical processes. The major IMF1 and IMF2 cycles were 2.6–3.8 and 5.5–8.3 y, respectively, and their contributions account for 59.01% (⁠$45.51+13.50%$⁠) of the total variance. These cycles principally reflect the interannual variations in the summer temperatures. The 2–8-y cycle was consistent with the ENSO cycle, which is an important influencing factor of interannual climate change [58]. Additionally, interannual variations in summer temperatures in China are largely affected by the Asian summer monsoon. Some studies have found that the Asian summer monsoon is influenced by the interannual variability of the Indian Ocean Basin sea surface temperature (SST) anomalies. These anomalies are also referred to as the Indian Ocean Basin mode (IOBM) [59, 60]. The 2–8-y interannual variations exhibited in IMF1 and IMF2 also likely correspond with the IOBM, which has been confirmed to affect the SETP [61]. However, according to the spatial correlation results between the reconstructed temperature and global SST (Figure 7), areas with significant positive correlations (⁠$P<0.05$⁠) were mainly concentrated in the tropical Indian Ocean, tropical Atlantic Ocean, regions of the Antarctic Circumpolar Current, the area near the western Pacific warm pool, and the Kuroshio region. The relationship between the reconstructed temperature and SST in the El Niño area was nearly insignificant, and the correlation coefficient between the reconstructed temperature with the IOBM and ENSO [58] was 0.236 and 0.057 (1854–1998), respectively. This may indicate that IOBM has a more important role in the change in summer temperatures in the SETP than ENSO. Recently, many studies have reported that the interannual and interdecadal variation of the IOBM may exert an important influence on Asian and even Northern Hemisphere climates [62–68]. During the positive phase of the IOBM, a new atmospheric heating source can be induced to enhance the Indian summer monsoon [67], and the local Hadley circulation is enhanced. In addition, the vertical descending motion over the SETP is stronger, causing above-normal temperatures in this area [61]. In contrast, the negative phase of IOBM resulted in below-normal temperatures.

The decadal variabilities shown in IMF3–IMF4 may be related to the 11-y Schwabe cycle and 22-y Hale solar activity cycles, respectively [69]. Many studies have found that solar activity affects climate change [70–73]. The 67.3-y cycle of IMF5 may be connected to the Atlantic multi-decadal oscillation (AMO), which has been identified as an ongoing series of multi-decadal changes occurring in the SST of the North Atlantic Ocean with an estimated period of 65–80 y [74]. The AMO can cause a change in the North Atlantic circulation, which can affect the westerlies and wave trains that modulate the surface temperatures and air pressures in the North Atlantic [75]. Therefore, the summer temperature variability in the SETP may be affected by the AMO because of its downstream wave-train guidance mechanism [76]. The IMF6 component exhibited periodic oscillations on multi-centennial scales, with a relatively small contribution of 5.33%, which may indicate the influence of the Gleissberg cycle [69]. The RES component exhibited a decreasing temperature trend after 1730 and an increasing trend after 1850. Moreover, as indicated by the IMF6 and RES components on a multi-century scale, the current temperatures are in the warming phase. This warming trend has also been observed in some SETP temperature reconstructions [77–79].

Based on a tree-ring network of 53 tree-ring chronologies (one δ¹³C chronology, one MXD chronology, and fifty-one width chronologies), the May–June temperature field of 28 meteorological stations located over the SETP was reconstructed following the modified PPR method. The common period was 1730–1998, and the longest single period was 1480–2002. The correlation coefficients, explained variance, reduction in error, principal component analysis, and rotated principal component analysis indicated that the reconstruction could capture the spatiotemporal characteristics of regional temperature variability. The temperature changes in the SETP were mostly similar, suggesting that relatively warm periods occurred during AD 1730–1754, 1814–1829, 1853–1886, 1910–1932, 1955–1966, and 1981–1998, while relatively cool periods occurred during 1755–1788, 1799–1813, 1830–1852, 1887–1909, 1933–1954, and 1967–1980. Using the EEMD method, six main intrinsic mode functions and one residual term were extracted to show periodic climate signals at different time scales in the reconstructed temperature series. These oscillations may reflect the influence of the Indian Ocean Basin mode, solar activity, and Atlantic multidecadal oscillation on SETP temperature changes.

Dataset is available upon request to the corresponding author or the first author.

The authors declare no conflict of interest.

This research was jointly funded by the Sichuan Central Government guides local science and technology development project (No. 2021ZYD0019), the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (No. 2019QZKK0103), and the National Natural Science Foundation of China (NSFC) (Nos. 41772173 and 41775070).