Abstract Large-scale climate history of the past millennium reconstructed solely from tree-ring data is prone to underestimate the amplitude of low-frequency variability. In this paper, we aimed at solving this problem by utilizing a novel method termed “MDVM”, which was a combination of the ensemble empirical mode decomposition (EEMD) and variance matching techniques. We compiled a set of 211 tree-ring records from the extratropical Northern Hemisphere (30–90°N) in an effort to develop a new reconstruction of the annual mean temperature by the MDVM method. Among these dataset, a number of 126 records were screened out to reconstruct temperature variability longer than decadal scale for the period 850–2000 AD. The MDVM reconstruction depicted significant low-frequency variability in the past millennium with evident Medieval Warm Period (MWP) over the interval 950–1150 AD and pronounced Little Ice Age (LIA) cumulating in 1450–1850 AD. In the context of 1150-year reconstruction, the accelerating warming in 20th century was likely unprecedented, and the coldest decades appeared in the 1640s, 1600s and 1580s, whereas the warmest decades occurred in the 1990s, 1940s and 1930s. Additionally, the MDVM reconstruction covaried broadly with changes in natural radiative forcing, and especially showed distinct footprints of multiple volcanic eruptions in the last millennium. Comparisons of our results with previous reconstructions and model simulations showed the efficiency of the MDVM method on capturing low-frequency variability, particularly much colder signals of the LIA relative to the reference period. Our results demonstrated that the MDVM method has advantages in studying large-scale and low-frequency climate signals using pure tree-ring data.

Citation: Xing P, Chen X, Luo Y, Nie S, Zhao Z, Huang J, et al. (2016) The Extratropical Northern Hemisphere Temperature Reconstruction during the Last Millennium Based on a Novel Method. PLoS ONE 11(1): e0146776. https://doi.org/10.1371/journal.pone.0146776 Editor: Juan A. Añel, Universidade de Vigo, SPAIN Received: April 30, 2015; Accepted: December 22, 2015; Published: January 11, 2016 Copyright: © 2016 Xing et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All relevant data are within the paper and its Supporting Information files. All tree-ring raw measurement data used in this paper are available from the International Tree-Ring Data Bank (http://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring). Funding: This study is supported by National Natural Science Foundation of China (41175066), China Postdoctoral Science Foundation (2014M550711), National Natural Science Foundation of China (41275076), and China Meteorological Administration Special Public Welfare Research Fund (GYHY201306019). The funders played no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Knowledge of climate change in the past is essential for understanding the response of the climate system to natural and anthropogenic forcings, and also useful for evaluating the observed 20th century warming objectively in a long-term perspective. As the worldwide instrumental records rarely extend back further than the mid-19th century, pre-instrumental data of temperature is often deduced from proxy records derived from natural archives. In recent years, a number of proxy-based high-resolution temperature reconstructions on hemispheric scale have been carried out (e.g. [1–12]). Nevertheless, a great range has been observed in the reconstructed amplitudes among the currently existing temperature reconstructions for the Northern Hemisphere (NH). In particular, reconstructions based solely on tree-ring data are prone to underestimate the temperature amplitude during the Medieval Warm Period (MWP) and the Little Ice Age (LIA) (e.g. [3]). In order to better capture the low-frequency variability, lower temporal resolution proxies (e.g. ice cores, sporopollen, speleothems and lake sediments, etc.) are usually used to combine with tree-ring data for large-scale temperature reconstructions (e.g. [7, 12]). However, the quality of these non-annually resolved records is hard to assess since they are difficult to calibrate against the annual meteorological data. Thus, the combination of tree ring and other proxies in low resolution may bring about significant uncertainties in the final reconstructions [13, 14]. Fortunately, some recent studies have demonstrated the feasibility of reflecting longer-term temperature trends using pure tree-ring data, provided that data are appropriately processed in respects of detrending in chronology developing, data screening, and statistical techniques [4, 10, 15, 16]. Firstly, it is known that different detrending methods usually result in obvious differences of the final chronologies in signal extracting. Compared with traditional data-adaptive curve-fitting methods (e.g. negative exponential curve or smoothing spline), the superiorities of Regional Curve Standardization [17] and signal-free method [18] have been fully demonstrated. Both methods can avoid the “segment length curse” [19], and capture more low-frequency variance in excess of the mean length of individual samples used in chronology development. Secondly, with regard to data screening, some criteria (e.g. length of the candidate chronologies, and threshold of the correlation with local/large-scale instrumental series) are usually employed in the screening process to ensure that the temperature signal is as strong as possible in the final large-scale temperature reconstructions [4, 5, 7, 9]. Previous research definitely emphasized the importance of the quality of available proxy data rather than quantity for further robustness and improvements in large-scale proxy-based reconstructions [20]. Nevertheless, as the criteria of screening procedure are more thorough and strict, fewer records could pass the screening process, easily leading to uneven spatial distribution of the tree-ring data and relatively poor spatial representativeness of the final reconstructions [1]. Thirdly, in terms of the statistical techniques used for calibration in large-scale reconstructions, the traditional method is based on the linear regression model (ordinary least squares). There are other two main methods based on variance matching (e.g. [8, 21]) and error-in-variables (EIV) regression model (total least squares) (e.g. [9, 22, 23]). The comparisons of considerable reconstructions have suggested that the application of inverse-regression-based methods, like those used in some researches [2, 3], is apt to result in large underestimation of low-frequency NH temperature variability [15, 24]. As indicated above, it is thus important and feasible to explore an effective method for large-scale temperature reconstruction on the basis of improvements in terms of chronology developing, data screening and calibrating procedures. In this paper, a novel method applicable to pure tree-ring data was developed in order to fully make use of the tree-ring records and better capture low-frequency variability. The method we proposed was termed “MDVM”, as its core was based on ensemble empirical mode decomposition (EEMD) and variance matching techniques. Its detailed steps would be expatiated in the following sections. The objectives of this study were (1) to introduce the MDVM method, (2) to evaluate its performance by establishing a reconstruction of the extratropical NH mean temperature for the past millennium, and (3) to analyze the characteristics of reconstructed temperatures in terms of trend, amplitude and response to external forcings.

Data Instrumental data We made use of the monthly land surface air temperature dataset (CRUTEM4) [25] for the period 1850–2000, with a horizontal resolution of 5°×5° longitude-latitude grid. The extratropical (30–90°N) land-only mean annual (January-December) temperature anomaly series relative to 1961–1990 was extracted as the target series for subsequent reconstruction. Tree-ring data Considering that different detrending methods would bring about discrepancy of final chronologies especially on the low-frequency variation, tree-ring raw measurements data rather than the ready-made chronologies was compiled from the International Tree-Ring Data Bank (ITRDB, http://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring). We also collected several chronologies (a few were reconstructed series) from some published literatures related on regional temperature signal. The length of these tree-ring records (width data with limited density data) was mostly longer than 800 years. Additionally, we usually utilized the most up-to-date data rather than the old records. Model simulation data For further analyzing, we extracted extratropical land-only mean annual simulated surface temperature from the outputs of 7 climate models for Last Millennium experiment in CMIP5 (Phase 5 of the Coupled Model Intercomparison Project) [26]: CCSM4 [27], HadCM3 [28], MPI-ESM-P [29], FGOALS-s2 [30], BCC-CSM1.1 [31], IPSL-CM5A-LR [32] and CSIRO-MK3L-1.2 [33]. These results are available via the Earth System Grid Federation (http://pcmdi9.llnl.gov/esgf-web-fe/). As described in CMIP5, the Last Millennium transient simulation is conducted by coupled climate system model (or earth system model) with external forcings and uniformly spans from 850 to 1850 AD. The simulated series were simply adjusted by removing their respective mean of the whole period (850–1850 AD) for subsequent comparison with the final reconstructed NH temperature.

Methods Chronology development The signal-free Regional Curve Standardization (SF-RCS) method [18] was used for the collected tree-ring raw measurements data. In each sampling site, the individual series were aligned by cambial age to estimate the site-specific mean biological growth curve. It has been borne out that SF-RCS is adept at preserving the low-frequency trends and can better mitigate the trend distortion problem by several recent studies on regional climate reconstructions [34–36]. The reliable period of each chronology we established was determined by expressed population signal (EPS), and the chronologies were truncated when their EPS values were less than 0.80 in consecutive 40 years [37]. MDVM method The two keys of MDVM method are ensemble empirical mode decomposition (EEMD) and variance matching. EEMD is a data self-adapted signal processing approach, which can relief the mode-mixing problem and bring about more stable and physically interpretable results [38]. EEMD is a powerful tool for decomposing complicated dataset into a finite number of monocomponent intrinsic mode functions (IMFs). Through large numbers of trials to add white noise, only the persistent surviving IMFs could be extracted, which should contain more physical meanings [38]. EEMD has been widely used in climate research [39–41] and also applied in a few dendrochronological studies [42, 43]. Variance matching method was employed in the calibration to adjust the mean value and variance agreeing with the target series over the common period [8, 21]. It has been proved that variance matching is able to avoid the problem with underestimation of low-frequency variability associated with direct-regression-based calibration methods [15, 24]. Some researches pointed out that the indirect-regression method can avoid this problem in more effective way, and this variance matching method should fall between the direct- and indirect- regression methods [6, 24]. As shown in Fig 1, the itemized steps of MDVM method were as follows: PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Flow chart illustrating the MDVM method developed in this study for the extratropical Northern Hemisphere temperature reconstruction. https://doi.org/10.1371/journal.pone.0146776.g001 Stage 1. Decompose each series into four components on annual, decadal, multi-decadal and centennial timescales, respectively. The instrumental extratropical NH mean temperature and each tree-ring chronology (or reconstructed series) was decomposed into IMFs using the EEMD method, respectively [Eq (1)]. Then, the power spectral analysis was applied for each IMF and the residual series to diagnose the dominant period, upon which to classify the IMFs (or residual series) into certain components on four typical timescales, i.e. annual scale component (less than 10 year), decadal scale component (10–30 year), multi-decadal scale (30–90 year) and centennial scale component (longer than 90 year) [Eq (2)]. (1) (2) (3) where N i was the total number of IMFs for i-th proxy (or instrumental CRU) time series [Ser i (t)]. Subscript j denoted the timescale of each component [comp i (t)]: j = 1 for annual scale, j = 2 for decadal scale, j = 3 for multi-decadal scale, and j = 4 for centennial scale. Stage 2. Screen individual components with corresponding benchmark series on decadal, multi-decadal and centennial scales, respectively. Considering that the regional tree-ring record coincided with the hemispheric mean temperature mainly on decadal and longer timescales, the decomposed variability of tree ring and temperature series on annual scale was ignored in this study, which meant that the smoothed series would consist of 3 components [Eq (3)]. Correlation analysis was employed between tree-ring and the benchmark series on decadal, multi-decadal and centennial scale, respectively. It is worth noting how to accurately estimate the threshold of significant correlation. Because the filtered time series have a reduction on the degree of freedom (DOF) comparing with the unfiltered data, more vigorous statistical tests should be utilized to demonstrate the statistical significance. In this study, we applied the Monte Carlo simulation method [44], which compared the actual correlations with correlations between instrumental temperatures and the simulated random time series generated from the original tree-ring series. The white noise was superimposed to the original series so as to generate the samples and these series would also be decomposed into components on the four timescales (details of Monte Carlo simulation please see the Supporting Information). The threshold for the value of the significant correlation coefficient was estimated using 5000 member Monte Carlo simulations with bootstrap resampling (400 trials) at a 95% level of confidence. The tree-ring component would be selected to establish the final reconstruction, if its correlation with the benchmark series on the corresponding timescale exceeded the threshold value. It was noteworthy that, owing to its relatively short length, the CRU temperature series was not able to exhibit complete signal on the centennial scale. To better reconstruct the long-term variability, an extra low-frequency series of NH temperature reconstructed in paper [8] (hereafter Moberg05-LF for short) spanning from 133 AD to 1925 AD, was used instead of CRU variability on the centennial scale, which was largely different from the traditional approach. Moberg05-LF rather than other low-frequency series was selected owing to the following two reasons. Firstly, Moberg05-LF, which was inferred from low resolution proxy records (pollen, borehole, stalagmite, etc.), was totally independent from tree-ring data. Secondly, the timescale of Moberg05-LF variability was longer than 80 year, which was similar to the centennial scale defined in this study. Stage 3. Composite the screened out components, and after scaling on corresponding timescales, combine them into the final reconstruction. Weight each screened tree-ring component on the same timescale in proportion to their explained variance for the benchmark series (i.e. CRU and Moberg05-LF) [43]. Then, composite these screened out components by weight on decadal, multi-decadal and centennial scale, respectively. The function was as follows [Eq (4)]: (4) where the subscript i, j denoted i-th standardized time series and j-th timescale; N j was the total number of screened components at j-th timescale, and CSj(t) was the corresponding composite series. The weight was given by the square of each correlation coefficient (r i,j ) relative to the benchmark series at j-th timescale. After that three composite series were obtained, and then scale them with the benchmark series (CRU and Moberg05-LF) on the corresponding timescales, respectively, by variance matching method [Eq (5)]. Finally, combine the scaled three composite series on decadal, multi-decadal and centennial scale, and the final reconstruction of the extratropical NH mean temperature was established [Eq (6)]. (5) (6) where reconj(t) represented reconstructions at different timescales. M Pcal and M Ical were mean values of proxy composite series and instrumental series during the calibration period at j-th timescale. and were the standard deviation of proxy composite series and instrumental series at j-th timescale. Validation of MDVM reconstruction Pseudo-proxy experiment (PPE). Pseudo-proxy experiment (PPE), as a more objective approach, has been widely applied to test reconstruction methods based on the climate model simulated results. The main purpose of PPE is to put the reconstruction methods to a common framework, and the model simulated temperatures can provide a longer timescale background for testing the low-frequency signal. In this research, PPE has been performed basing on the ensemble results of the CMIP5 last millennium simulations and their corresponding extended simulations. In total, seven sets of results from five models were included: BCC-CSM1.1, CSIRO-Mk3L-1.2, HadCM3, MPI-ESM-P, and GISS-E2-R (p121, p124 and p125, http://data.giss.nasa.gov/modelE/ar5/). The corresponding annual mean surface temperature of the ensemble was extracted to be the known model target. The experiment was designed as follows. Firstly, using the bilinear interpolation scheme, the model simulation results were interpolated onto the 126 tree-ring sites (i.e. as same as the tree-ring records used for the real-world MDVM reconstruction, see S3 Table). Then Gaussian white noise with signal-to-noise ratio (SNR) of 0.1, 0.25, 0.5 and 1.0 was respectively added to the interpolated series to mimic the pseudo-proxy series over the entire 1150 years (851–2000 AD). Following the same steps of MDVM method (Stage 1–3), the variability on annual scale was discarded, and the reconstructed composite series at other three typical timescales (decadal, multi-decadal and centennial scale) were calculated basing on the pseudo-proxy series from the same proxy networks. The only difference when performing PPE was that the target series for calibration/validation was the ideal series of pure simulated NH mean surface temperature (851–2000 AD) instead. Two trials of verification tests with exchanged calibration and validation periods were conducted: Exp-A with 851–1500 AD as calibration period and 1501–2000 AD as validation period; Exp-B with 1501–2000 AD as calibration period and 851–1500 AD as validation period. Furthermore, in order to evaluate the performance of MDVM on centennial timescale, two extra metrics were calculated as well, i.e. the ratio of standard deviation of reconstructed centennial variability to simulated counterpart (low frequency stand deviation ratio, hereafter LF-Std-Ratio for short) [Eq (7)]. (7) Leave-one-out cross-validation. The “leave-one-out” cross-validation [45] was performed to verify the accuracy of MDVM reconstruction against the instrumental temperature data (CRU) over the period 1850–2000. The following statistics were examined in the validation: reduction of error (RE), coefficient of efficiency (CE), and root-mean-square error (RMSE). The uncertainty of the reconstruction was estimated using the standard deviation (std ver ) of the instrumental extratropical NH temperature and the correlation coefficient (r) between the reconstructed and instrumental series in validation period [46] [Eq (8)]: (8) Monte Carlo simulation process was also applied to test the significance of RE, CE scores and the uncertainty, and samples for the testing were generated by first-order autoregressive model (AR1) plus red noise [9].

Conclusion We introduced a novel method termed “MDVM”, which was a combination of the EEMD and variance matching techniques. This method allowed proxy data to take decomposition, screening, and calibration on separating timescales. Utilizing the MDVM method, a new reconstruction of the extratropical NH mean temperature for the period 850–2000 AD was established. 126 records with relatively fine spatial coverage were finally screened out and utilized to reconstruct temperature variability longer than decadal scale. The MDVM reconstruction depicted significant low-frequency variability in the past millennium with evident LIA and MWP. Temperature in 20th century was likely unprecedented in an 1150-year context. Additionally, the MDVM reconstruction covaried broadly with changes in natural radiative forcing, and particularly showed distinct footprints of multiple volcanic eruptions in the last millennium. The comparisons with other versions of NH temperature reconstructions (using same dataset by different approaches, and previous versions developed by other people), as well as some model simulation data, further demonstrated the efficiency of MDVM method on studying large-scale and low-frequency climate signals using pure tree-ring data.

Acknowledgments We gratefully acknowledge the contributors to the International Tree-Ring Data Bank as well as all the scholars that willingly shared their published but un-archived proxy data with us. We also thank all the reviewers for their helpful suggestions on the quality improvement of this paper.

Author Contributions Conceived and designed the experiments: YL PX XC. Performed the experiments: PX XC. Analyzed the data: PX XC SN ZZ JH. Contributed reagents/materials/analysis tools: PX XC JH SW. Wrote the paper: PX XC YL SN ZZ.