Multi-Proxy Temperature Reconstruction from the West Qinling Mountains, China, for the Past 500 Years

A total of 290 tree-ring samples, collected from six sites in the West Qinling Mountains of China, were used to develop six new standard tree-ring chronologies. In addition, 73 proxy records were assembled in collaboration with Chinese and international scholars, from 27 publically available proxy records and 40 tree-ring chronologies that are not available in public datasets. These records were used to reconstruct annual mean temperature variability in the West Qinling Mountains over the past 500 years (AD 1500–1995), using a modified point-by-point regression (hybrid PPR) method. The results demonstrate that the hybrid PPR method successfully integrates the temperature signals from different types of proxies, and that the method preserves a high degree of low-frequency variability. The reconstruction shows greater temperature variability in the West Qinling Mountains than has been found in previous studies. Our temperature reconstruction for this region shows: 1) five distinct cold periods, at approximately AD 1520–1535, AD 1560–1575, AD 1610–1620, AD 1850–1875 and AD 1965–1985, and four warm periods, at approximately AD 1645–1660, AD 1705–1725, AD 1785–1795 and AD 1920–1945; 2) that in this region, the 20th century was not the warmest period of the past 500 years; and 3) that a dominant and persistent oscillation of ca. 64 years is significantly identified in the 1640–1790 period.


Introduction
Temperature-sensitive proxy data can be used as a primary basis for understanding temperature variations through time. Tree-ring datasets are particularly important natural proxy records used in climatological research, as they provide accurate dating, annual-scale resolution, and are available from locations distributed worldwide [1]. A major focus in dendroclimatology has been the reconstruction of temperatures at a single location using local tree-ring samples [2][3][4][5][6][7][8][9][10][11][12][13][14][15].
Despite the advantages of tree-ring records as temperature proxies, tree-ring chronologies often include less low-frequency information than do other proxies because of the so-called segment-length curse problem [16], which refers to the difficulty of preserving cyclic signals that are longer than the duration of the age of individual trees. Samples from individual trees rarely span time frames longer than a millennium, and are usually much shorter. The standardization process which is used to eliminate the effects of individual tree growth patterns, tends to weaken lowfrequency climate signals [17]. Although dating errors of other proxy types are usually greater than those obtained from tree ringdata, other proxies are often more accurate in preserving lowfrequency temperature signals [18]. Therefore, in this study, we employed a hybrid reconstruction method that considers the full spectrum of temperature variability, by separating the variability into high-and low-frequency bandwidths. Every proxy record contains both climate-related information and non-climate-related ''noise''. In the proxy records of individual trees, the noise can be large and can even dominate the climate signal; however, when averaging a number of proxies, this noise is reduced (assuming that the noise is random). By assembling a large number of proxy records from the same region, we are more likely to extract a common climate-related signal, and the signal is likely to be less strongly affected by the noise inherent in individual proxies.
The correlations that exist between proxy data and instrumental data in Asia are complex [19]. For example, tree-ring width chronologies from the northeastern Tibetan Plateau (Dulan County), where elevation gradients are large, have been used to reconstruct precipitation [6,20], temperature [3,21], and the Palmer drought severity index (PDSI) [22]. Trees growing near the upper treeline typically respond to temperature, whereas those near the lower treeline respond to precipitation [23]. Furthermore, the complex monsoon system and the thermal forcing of the Tibetan Plateau may strongly influence Asian climate [24]. By using a large number of proxies to rebuild the response function through an ensemble reconstruction approach that averages a large number of ensemble members, the accuracy of the reconstructions may be improved [19].
We developed six new tree-ring chronologies from the West Qinling Mountains, and assembled 67 other proxy records in China. Using a modified point-by-point regression (hybrid-PPR) method, we reconstructed the annual mean temperatures over the past 500 years. Special attention was paid to the reconstruction method, as it differs slightly from those used in previous studies [19,25].

Materials and Methods
he West Qinling Mountains (33u-35uN, 102u-105uE) is situated near the boundary of the East Asia monsoon region, in the transition zone between the Tibetan Plateau and the Loess Plateau of China [26]. A multi-proxy approach is required to accurately reconstruct past climate variability in this region because of the complexity of the climate regime and the variability associated with different climate signals from this area. No specific permits were required for the described field studies. Permission to obtain tree-ring samples was granted by the Bailong River Forestry Bureau to the Lanzhou University research team headed by Dr. Fengmei Yang (first author). The study area belong to the Bailong River Forestry Bureaus This location is not privately owned or protected in any way, and the field studies did not involve endangered or protected species.
For the CRUTEM3v dataset we used a grid cell centred at 32.5uN and 102.5uE, and for the CRUTS3.1 dataset we used the four grid cells nearest to the West Qinling Mountains. We also used a grid cell centred at 33.75uN 103.75uE, as a correlation analysis between the CRUTS3.1 dataset and the averaged 12meteorological-station dataset showed this grid cell to have the highest correlation coefficient for the period 1951-1995 (r = 0.85). In Figure 2A, the three instrumental temperature datasets seem similar, but the amplitude of the low-frequency variability differs greatly, as shown in Figure 2B. The amplitude of the temperature variability in the CRUTEM3v data after AD 1966 is clearly greater than is evident in the averaged data of the 12 meteorological stations. Similarly, the amplitude of the tempera-ture variability in the CRUTS3.1 data after AD 1957 is slightly greater than that observed in the averaged data of the 12 meteorological stations. The variability between the different instrumental datasets is of importance, as the differences have a strong impact on the reconstruction results.

Temperature proxy data
The West Qinling Mountains are relatively free from human disturbance and the trees are mostly healthy, thus providing ideal conditions for obtaining for dendroclimatological records from tree-ring samples. During 2009 and 2010, we collected 290 treering samples obtained as increment cores at six sites in the West Qinling Mountains, and from these data, we developed six new multi-centennial tree-ring width chronologies (see Figure 1). The ring widths in the cross-dated samples were measured and the initial dating (and associated errors) of the samples were crosschecked using the quality-control software COFECHA to check the cross-dating and overall quality of the tree-ring chronologies [28]. The ring-width measurement series were detrended by applying a negative exponential curve; some tree-ring width series which did not exhibit negative exponential growth trends were excluded when the chronologies were developed. The chronologies were calculated from the bi-weighted robust means of the raw measurement data using the dendrochronology program library in R (dplR) [29]. The standard version of the six chronologies in which the expressed population signal (EPS) exceeds a threshold value of 0.85 [30] was assessed for coherence using correlation analysis.
All records were required to have decadal (or higher) temporal resolution. With the exception of some of the tree-ring chronologies, all of the proxies in this study have been previously used to reconstruct temperature variability records and all have been validated as temperature proxies [41][42][43]. We used only tree-ring records showing a strong correlation with the annual mean temperature of the 12 stations during the period AD 1951-1995, at a confidence level of 90% or above. All tree-ring chronologies were required to span the interval AD 1800-1980. Proxies with non-annual temporal resolutions were linearly interpolated to an annual resolution; this interpolation introduced a strong bias by adding red noise to the high-frequency components of the spectrum, thus reducing the verification skill because of the higher autocorrelation [49,50]. The non-annually resolved proxy data were therefore used only to reconstruct low-frequency temperature variability.

Reconstruction methods
The temperature reconstruction method used in this study was a modified PPR method. The PPR is a straightforward principal component regression that has been successfully used to produce high-quality reconstructions of droughts in North America [51], and for generating the North American Drought Atlas [52] and the Monsoon Asia Drought Atlas [19]. The PPR is based on the principle that only proxy records close to the location being reconstructed are good predictors for that location; this premise is particularly important for places with a spatially heterogeneous climate, such as China. The traditional PPR approach uses two variables: search radius and screening probability. The search radius defines the maximum tolerable distance between potential proxy records and the location of the reconstruction. The 450 km search radius was found to be the optimal search radius [53]. The four search radii for the Palmer Drought Severity Index (PDSI) in Asia, 500 km, 1000 km, 2000 km and 3000 km were chosen in that study [19]. The reason for this is the irregular distribution and relatively high noise level of many tree-ring chronologies in China. In this work, search radii of 100, 300, 500 and 1000 km were used for the temperature reconstructions, as smaller search radii are likely to yield a more accurate representation of the local climate information.
The screening probability defines the correlation threshold for including proxy series in the temperature reconstruction. First, all proxies were screened against instrumental data and only those with correlations at or above a 90% confidence level were retained. Because of the high level of noise in much of the proxy data from China, the same weighted model was applied to the records in this study, as was implemented in the creation of the Monsoon Asia Drought Atlas [19].
Our approach differs from the traditional PPR approach by incorporating a ''hybrid'' procedure that separately calibrates information in ''inter-decadal'' and ''inter-annual'' bandwidths; inter-decadal bandwidths are those with periods longer than 10 years (including proxies with annual or decadal resolution) and inter-annual bandwidths are those with periods shorter than or equal to 10 years (including those with only annually resolved proxies). This process is derived from the hybrid frequencydomain calibration [50]. The final annual temperature reconstructions were assembled using the two frequency results. Our approach further differs from the traditional PPR in that we used the regularized errors-in-variables (EIV) method to allow for the assimilation of non-tree-ring data, instead of using the traditional principal component regression method. The steps are described in more detail below. 3.1 Infilling the proxy data. We used a cubic spline to interpolate proxy records with non-annual-to annual-scale resolutions. Proxies that did not extend to AD 1996 were extrapolated to AD 1996 using the regularized expectation maximization (RegEM) algorithm, based on their mutual covariance with the other available proxies over the 1800-1996 period [54]. In this study, 18 tree-ring chronologies, two ice-core records (No. 52 and 53 in Table 1) and a Chinese documentary record (No. 51 in Table 1) were extrapolated. Finally, to avoid biases due to different temporal resolutions and interpolation procedures, all records were filtered to retain frequencies of f,0.1 cycles per year in the reconstruction.
3.2 Screening of proxy data. The screening probability in original PPR frame is a method for selecting proxy records for use in the temperature reconstructions at a particular grid point, based on a threshold defined by a certain significance level for correlation [53]. Firstly, every proxy record was required to show a strong correlation (.90% confidence level; n = 45) with the annual temperature of the 12 meteorological stations for the period AD 1951-1995. Because of the high level of noise in the tree-ring chronologies, the same weighted model with the same coefficients was applied in this study as was used in the creation of the Monsoon Asia Drought Atlas [19]. This is expressed as where P w is the weighted proxy data, P u is the original unweighted proxy data in standard normal deviate form, r is the correlation between the tree-ring chronologies and the instrumental temperature data over the calibration period 1951-1995, and p is some power, 0, 0.5, 0.67, 1, 1.5, 2, used to assign different weights to each proxy dataset, following the ensemble PPR method [19].
3.3 Splitting of proxy and instrumental data. Following the procedure described in the paper [55], the screened proxy and instrumental data were separated into high-and low-frequency bandwidths using a Butterworth IIR filter with a cut-off frequency of 0.1 cycle per year. The screened proxy records were used to reconstruct temperature in separate high-and low-frequency bandwidths. This frequency split makes it possible to assimilate proxy data series with different temporal resolutions.
3.4 Regressing of proxy data. The regularized EIV temperature reconstruction method [55] was used to reconstruct past temperature variability in the two frequency bandwidths, instead of the more traditional principal component regression method. These techniques are described in detail elsewhere [49,54,[56][57][58]. We acknowledge the considerable debate in the literature about the validity of the RegEM/EIV method; however, the method has been utilized in numerous previous studies for reconstructing temperatures on local to global scales [49,55,[58][59][60][61][62], and it is also used to reconstruct continental-scale temperature variability by the PAGES 2K network. We considere it to be a feasible method for rebuilding the transfer function. The code for the RegEM/EIV method is availible at http://www.meteo.psu.edu/holocene/ public_html/supplements/Multiproxy Means07/code/coderecon/. When applying small values of the search radius, the approach used here shares the reconstruction principles that form the basis of the PPR method [53]. Based on the EIV method, the regularized expectation maximization (RegEM) algorithm [49,58] was used to determine a matrix of regression coefficients B k at grid point k and model errorse k , as  T recon,k~Precon,k z(T cal,k zP cal,k ) : where T recon,k is a row vector giving the predicted temperature at grid pointk, T cal,k is a row vector giving the instrumental temperature during the calibration period, P recon,k is a matrix of the proxy data (as the predictors), P cal,k is a matrix of proxy data during the calibration period [53]; this equation is revised from this paper [49]. As a further safeguard against potentially nonrobust results, a minimum of seven predictors is required in implementing the EIV procedure, according to the process described in this paper [55]. This criterion regarding the number of minimum predictors is an empirical constant, which is a compromise between the number of available predictors and the requirement for the model stability. For this reconstruction, the minimum number of predictors is greater than 20.
3.5 Validation of reconstruction. The accuracy and skill of the grid-point reconstructions were assessed using a split-period approach [19,55]. Reconstruction models were validated for the period AD 1921-1950, and the fit of the reconstructions was measured by the square of the Pearson product-moment correlation coefficient (r 2 ), the reduction of error (RE), and the coefficient of efficiency (CE) [19]. The uncertainty of the reconstruction is expressed as in the paper [55], the code for this formula is available at Prof. Mann's website (www.meteo.psu.edu/ mann/supplements/MultiproxyMeans07/code/codeveri/ calc_error.m), and is expressed by where U is the uncertainty, S is the standard deviation of the instrumental data during the calibration period, and r 2 is the squared correlation coefficient, calculated from the data for the verification period.
A wavelet analysis was used to detect cycles in our reconstruction results. The source code for the wavelet analysis is available at this website (http://paos.colorado.edu/research/wavelets/ software.html) [63]. The mother function was set as the Morlet wavelet, which gives a high resolution for the periodicity and includes a complex exponential function modulated by a Gaussian wavelet [63]. A simple red-noise model representing an autoregressive linear Markov process [64] was chosen as an appropriate background spectrum for determining the 95% confidence levels of the global wavelet power spectrum. The autocorrelation parameter of red noise background is 0.72.

Results and Discussion
Reconstruction accuracy as measured by r 2 , RE and CE are shown in Table 2. We used four search radii (100, 300, 500 and 1000 km), and six weighted cases (powers of 0, 0.5, 0.67, 1, 1.5, and 2), giving 24 ensemble members. Our ensemble includes some members for which RE and CE exceed 0, as these values indicate that these parameters are of value [19]. Estimates of the uncertainty are included in Table 2, as based on equation 3. The CE values are mostly negative, indicating that the results are not well verified. The main reasons for this are the very short duration of the instrumental temperature data in the study area and the restricted ability of proxy data to extract the climate signals; thus, temperature is only one factor affecting tree-ring growth patterns in this region [65]. The last warm period (AD 1920(AD -1945) is likely to be the regional expression of 20 th first century global climate warming, the similar example that the Arctic air temperature warming during the AD 1910-1940 period [66].
Our reconstructions show that the 20 th century was not warmer than other periods in the past 500 years in this region of China. The phases of the three reconstructions are mostly consistent, but the reconstruction that was calibrated against the CRUTEM3v dataset shows larger amplitudes than do the other two reconstructions. It is important to note, however, that the choice of instrumental data has a significant effect on the amplitude of the reconstructed temperature variations during the warm periods in the 18 th century. For example, for the CRUTEM3v dataset, the instrumental temperature (filtered using a Butterworth IIR filter with a cut-off frequency of 0.1 cycle per year) exhibits a distinctly and larger amplitude than other datasets ( Figure 2B). When compared with the long-term means, the 1961-1990 reference period was relatively cold in the West Qinling Mountains, in contrast to most other regions in the Northern Hemisphere. Consequently, the reconstructed temperatures are above the 1961-1990 means for most of the past 500 years, which is the reverse of the situation observed in most other regions, including China as a whole. Figure 4 compares of our reconstruction with the other results [67]; both studies used the same instrumental data and CRUTEM3v for the past 500 years. Grey lines show the uncertainties. Two reconstructions are presented, using the different proxy datasets. There are only several proxy records in the study area used in the other reconstruction [67]; however, we used 73 proxy records in this study. The reconstructions show two cold periods at ca. AD 1610-1620 and ca. AD 1965-1985, and two strong warming periods at ca. AD 1645-1660 and ca. AD 1920-1945. Overall, the timing of some of the transition periods from warm to cold conditions in our results coincides with the other results, but the timings of some transition periods do not match. The amplitudes of our reconstruction are larger than those   [67] for the field temperature reconstruction of corresponding geographical regions (see Figure 4). The differences might be attributed to different methodologies and different proxy data. The PPR method has the advantage, at least in theory, of preserving local signals more accurately than do methods based on empirical orthogonal function regression [67]. In addition, in this study, we have used more local proxy data, which can preserve more distinct local climate variability, than the other result [67]. The use of a great number of local and high-quality proxy records will more realistically reflect the local climate variability in a given area.
A wavelet analysis was used to detect cycles in the reconstructions over the full reconstruction period AD 1500-1995. The results ( Figure 5) indicate that a dominant and persistent oscillation of ca. 64 years is present across almost the entire period, and is especially noticeable during ca. AD 1640-1790. This dominant timescale of climatic variability, which has been found in instrumental records [68], other proxy data [69] and simulation results [70], may originate from an internal variability of the ocean-atmosphere system, e.g. the Atlantic multi-decadal oscillation [71]. Figure 6 compares between our temperature reconstruction with other reconstructions of temperature [5,36,72], precipitation [33], the El Niñ o-Southern Oscillation (ENSO) [73], and the Pacific Decadal Oscillation (PDO) [74]. The temperatures reconstructed for Sunan County are representative of the temperature variability for the central Qilian Mountains [5,72]. The Qamdo temperature reconstruction is representative of the temperature variability for the southeast area of the Tibetan Plateau [36]. The precipitation reconstructed for Zhamashike correlates strongly with precipitation in other regions of the Qilian Mountains; this correlation is, significant at the 95% confidence level, indicating that this record is likely to represent a response to a common precipitation signal [33], and therefore used to represent precipitation variability in the Qilian Mountains. The ENSO variance reconstruction, which denotes the Niñ o3 region, is based on the North American Drought Atlas [73] and is derived from the tree-ring datasets in Northern America. The PDO index was developed from 17 Asian tree-ring chronologies [74], filtered to retain frequencies of f ,0.1 cycle per year.
The dominant 1705-1725 warm period ( Figure 6) coincides with a cold and dry period in the central Qilian Mountains and a warm period in the southeastern Tibetan Plateau and is related to lowered ENSO variance and raised PDO values. This period coincides with a weakened Asian summer monsoon [40], which is the main source of precipitation in the study area. Moreover, the temperatures during this period are greater than or equal to those in the 20 th century. These warm periods were evident in the southeastern Tibetan Plateau [36]. The other three warm periods could be associated with dry periods in the Qilian Mountains and a strengthened Asian summer monsoon. This indicates that the Asian monsoonal dynamics may be a major factor affecting climate variability in the study area, but it is not the only factor. The last warm period occurred during the period AD 1920-1945. A widespread drought event occurred in the 1920s and 1930s in northern China [75,76]; this can be observed in our reconstruction, although the duration of the event may have been longer in the study area, possibly representing divergences between local and regional environmental trends.
According to our reconstructions, the first three cold periods occurred at approximately AD 1520-1535, AD 1560-1575, and AD 1610-1620, all of which all belong to period of the Little Ice Age. The last cold period AD 1610-1620 can also be identified in the central Qilian Mountains and on the Tibetan Plateau, but the other cold periods can not be found (see the Figure 6). These patterns indicate that temperatures exhibited different phases and  Reconstructed temperature (blue line) in the Qamdo region of the Tibetan Plateau (blue line), calculated using the tree-ring regional curve standardization (RCS) chronology [36], and the reconstructed temperature (black line), using the tree ring chronology from the middle Qilian Mountains [5,72]. C) Reconstructed precipitation and temperature (green line), calculated using the tree-ring chronology from Zhamashike in the Qilian Mountains [33]. D) Reconstructed ENSO (orange line), calculated using tree-ring chronology [73], and the reconstructed PDO (black line), calculated from tree-ring chronologies [74], filtered to retain frequencies above a cut-off frequency of 0.1 cycle per year. doi:10.1371/journal.pone.0057638.g006 amplitudes in different areas during the Little Ice Age. The precipitation in the central Qilian Mountains was high during AD 1520-1535, but dry conditions prevailed during the next two cold periods. This indicates that warm periods and wet periods are not homogeneous in this area.
Temperatures during the second and third cold periods correspond to the relatively low variances of ENSO during the same periods; however, during the first cold period the temperatures correspond to a high ENSO variance. Moreover, the PDO index compared two lower values and a high value during the same periods. This shows that the temperature in the West Qinling Mountains was affected by multiple factors. The fourth-cold period, which is observed in the reconstruction for the middle of the Qilian Mountains, and the southeastern area of the Tibetan Plateau ( Figure 6), may represent a larger scale cold event; precipitation during this time was low. The ENSO variance and the PDO index show two cold phases, with cold periods evident in the PDO and ENSO records during the 1950s and 1960s, as observed in Figure 6; a cold and wet period also occurs during this interval. However, in contrast to the explanation given for the above mentioned cold periods, this cold period may be explained by a stronger Indian summer monsoon [77]. Thus, local climate anomalies such as these serve to illustrate the importance of increasing our understanding of past climate variability, not only at large spatial scales, but also at local and regional levels.

Conclusions
Tree-ring samples from six sites in the West Qinling Mountains of China were used to develop six new tree-ring width chronologies. In combination with 67 other proxy records of various types, the new chronologies were used to reconstruct temperatures in the West Qinling Mountains from AD 1500 to the present. Our results show that a hybrid-PPR method can effectively assimilate different types of proxy data with different temporal resolutions, and preserve more information about local climate variability than other similar techniques. The reconstruction shows 1) five distinct cold periods, at approximately AD 1520-1535, AD 1560-1575, AD 1610-1620, AD 1850-1875, and AD 1965-1985, and four warm periods, at approximately AD 1645-1660, AD 1705-1725, AD 1785-1795 and AD 1920-1945 during the past 500 years in the West Qinling Mountains; and 2) the 20 th century was not the warmest period during the past 500 years in the West Qinling Mountains; and 3) there is a dominant and persistent oscillation of ca.64 years during the period AD 1640-1790.