Modified climate with long term memory in tree ring proxies

Long term memory (LTM) scaling behavior in worldwide tree-ring proxies and subsequent climate reconstructions is analyzed for and compared with the memory structure inherent to instrumental temperature and precipitation data. Detrended fluctuation analysis is employed to detect LTM, and its scaling exponent α is used to evaluate LTM. The results show that temperature and precipitation reconstructions based on ring width measurements (mean &agr; = 0.8 ?> ) contain more memory than records based on maximum latewood density (mean &agr; = 0.7 ?> ). Both exceed the memory inherent to regional instrumental data ( &agr; = 0.6 ?> for temperature, &agr; = 0.5 ?> for precipitation) in the time scales ranging from 1 year up to 50 years. We compare memory-free ( &agr; = 0.5 ?> ) pseudo-instrumental precipitation data with pseudo-reconstructed precipitation data with LTM ( &agr; > 0.5 ?> ), and demonstrate the biasing influences of LTM on climate reconstructions. We call for attention to statistical analysis with regard to the variability of proxy-based chronologies or reconstructions, particularly with respect to the contained (i) trends, (ii) past warm/cold period and wet/dry periods; and (iii) extreme events.


Abstract
Long term memory (LTM) scaling behavior in worldwide tree-ring proxies and subsequent climate reconstructions is analyzed for and compared with the memory structure inherent to instrumental temperature and precipitation data. Detrended fluctuation analysis is employed to detect LTM, and its scaling exponent α is used to evaluate LTM. The results show that temperature and precipitation reconstructions based on ring width measurements (mean 0.8 α = ) contain more memory than records based on maximum latewood density (mean 0.7 α = ). Both exceed the memory inherent to regional instrumental data ( 0.6 α = for temperature, 0.5 α = for precipitation) in the time scales ranging from 1 year up to 50 years. We compare memory-free ( 0.5 α = ) pseudo-instrumental precipitation data with pseudo-reconstructed precipitation data with LTM ( 0.5 α > ), and demonstrate the biasing influences of LTM on climate reconstructions. We call for attention to statistical analysis with regard to the variability of proxy-based chronologies or reconstructions, particularly with respect to the contained (i) trends, (ii) past warm/cold period and wet/dry periods; and (iii) extreme events.

Introduction
Natural time series sometimes have a certain selfdependent structure, i.e. when an observation in a time series influences the following observations over a long period of time, this behavior is known as long term memory (LTM). Hurst (1951) first introduced the concept of LTM using the rescaled range (R/S) method in analyzing Nile river discharge. Since then, LTM has then been detected in precipitation, temperature and other climate and more general environmental data (e.g., Koscielny-Bunde et al 1998, Marani 2003, Bunde et al 2013, Yuan et al 2014. In mathematical terms, a time series with LTM has an autocorrelation function C(n) following a power law C n n ( ) ∼ γ , with n as time scales, γ as the autocorrelation exponent, and its mean correlation time diverges for infinitely long series. This scaling behavior can also be expressed as a fluctuation function F n n ( ) ∼ α with α as the DFA exponent, obtained from detrended fluctuation analysis (DFA), or the power law decay of the power spectrum function S f f ( ) ∼ β with f as frequencies and β as the power spectrum exponent. β, α, and γ can be used to measure LTM and their relations are: Values of 0.5 1 α < < indicate LTM; 0.5 α = an uncorrelated process; 0 0.5 α < < anti-correlation; 1 α = is characteristic of f 1 noise; and 1 α > marks a non-stationary and unbounded series.
LTM in global instrumental temperature data yields 0.5 α ≈ over inner continental areas in North America and central Asia, 0.65 α ≈ in coastal areas and 0.8 α > over ocean for a time range from months to a few decades (Koscielny-Bunde et al 1998, Bunde and Havlin 2002, Monetti et al 2003. Other work revealed distinct scaling regimes for precipitation over land, and generally smaller DFA exponents close to 0.5 α = for the time range higher than days (Fraedrich and Larnder 1993, Marani 2003, Kantelhardt et al 2006. Due to the restricted length of instrumental data, the range of LTM in precipitation and temperature time series has been evaluated for periods up to a few years. However, the power spectra are reproduced in millennium long simulation of coupled atmosphere-ocean general circulation model, and the simulated LTM extends up to decadal and centennial time scales . Furthermore, by using Greenland ice cores data and simulations, the LTM for surface temperatures can even be extended to millennium scale . Therefore, we assume no breaking down of the persistence law from inter-annual to centennial scales. As LTM is considered as a robust feature of internal atmosphere-ocean dynamics and independent on time-dependent external forcing (Hunt 1998 and the variability is considered as a integral of a number of nonlinear processes evolving on a wide range of time scales (Ashkenazy et al 2003), our analysis only focuses on the whole data period rather than individual periods.
Memory in proxy data has been identified almost half a century ago (Matalas 1962) by comparing firstorder autocorrelation in the tree-ring series with firstorder autocorrelation in theoretical precipitation/ drought reconstructions. It was then further addressed in Stockton and Boggess (1980) and Meko (1981) introducing methods to correct the spectral biases. For example, lagging schemes were used by Stockton and Boggess (1980) to adjust for differences in autocorrelation between tree ring and river discharge data. Meko (1981) considered Box-Jenkins methods based on moving-mean autoregressive processes ARMA, to correct the behavior of autocorrelations at certain time-lags. The method was subsequently used in other tree-ring based climate reconstruction including Cook et al (1999) in their drought reconstruction. However, the memory was addressed more like short-term, as the autocorrelation was only considered at certain time-lags. Recently several papers have pointed out LTM in proxy data (e.g., Fraedrich et al 2009, Bunde et al 2013, Franke et al 2013, Barboza et al 2014, Kim et al 2014. Franke et al (2013) report that proxy data, including tree-ring width (TRW) and maximum latewood density (MXD) chronologies, show different scaling characteristics including 0.75 α = for precipitation sensitive TRW chronologies, and 0.85 α = for temperature sensitive TRW chronologies. Bunde et al (2013) demonstrated the existence of LTM in TR based precipitation reconstructions. Their work revealed TR chronologies to contain larger LTM compared to observed temperature and precipitation. Since the power spectrum exponents β in Franke et al (2013) may contain relatively high uncertainties (Kim et al 2014), and Bunde et al (2013) and other studies only analyzed few selected proxies or temperature/ precipitation reconstructions, the character and range of LTM in TR proxies (TRW and MXD) still need to be assessed.
Here, we aim at (1) confirming the existence of climate-relevant or -irrelevant processes introducing LTM in TR data; (2) analyzing the scaling behavior in TR based climate reconstructions; (3) assessing the potential influences of modified LTM on climate reconstructions based on pseudo-proxy data. We therefore assess the LTM in a total of 231 world-wide distributed climate sensitive TR chronologies, as well as regional mean climate reconstructions from China by employing the DFA. China is an important region to reconstruct past climate, because it has a wide-range of climatic conditions and offers a variety of proxy data, primarily from tree-rings (e.g. Cook et al 2010).
To demonstrate the impact of LTM on reconstructed climate variability, we compare pseudoinstrumental data (without memory) and pseudoreconstructions (with artificial LTM) derived from simulated (1155 years) precipitation in southeastern Asia (see section data for more details). For detailed temporal comparisons, we give an example using the data at grid-point [105E, 23.3N] for the full period . This grid point is located at the border of in Yangtze River region, one of the most important East-Asia Monsoon regions. For the spatial comparison of trends and anomalies, we show the evaluation from the pseudo-data in the period 1901-2005 years over southeastern Asia.

Data
Proxy and reconstructed data We use 231 TR chronologies from across the globe, which include 133 TRW and 98 MXD site records starting before 1700 AD and ending in the late 20th century (mostly in 1998 AD). 15 (118) TRW datasets are interpreted as proxies for annual temperature (precipitation) variation; see table S1 for detailed information and figure S1 for their locations. All TRW chronologies are directly extracted from international tree ring data bank (ITRDB, Grissino-Mayer and Fritts 1997) and MXD mostly from Cook et al (2010). They are significantly correlated with gridded instrumental data (details see 'table S1 and 'table S2 in the appendix of Briffa et al 2002). Besides, we also use one TRW chronology representing precipitation in northeastern Tibet over the past 3,500 years (Yang et al 2014, here referred to as 'Yang2014), and reconstructions of regional mean temperature/PDSI from various recent publications from China and East Asia (Cook et al 2010, Shi et al 2012, Ge et al 2013, and Yang et al 2002, here referred to as 'Yang2012', 'Cook2010', 'Shi2012-EIV', 'Ge2013', and 'Yang2002', respectively. 'Cook2010' is the regional mean of gridded annual PDSI reconstructions over East Asia. 'SHI2012-EIV' is the regional mean of gridded annual temperature reconstructions over China. 'SHI2012-EIV' is a decadally resolved temperature reconstruction over China derived from merging multiple proxies from nine different climate regions using an area weighting scheme. 'Ge2013' is also a decadally resolved regional mean temperature reconstruction over China derived from multi-proxies using a principal component regression method. 'Yang2012', 'Cook2010' and 'Shi2012-EIV' use only or primarily TR data; 'Yang2002' and 'Ge2013' use multi-proxy data.

Pseudo-instrumental and pseudo-reconstructed rainfall
The pseudo-instrumental precipitation data is directly extracted from model simulations at 1.875 1.875°×°r esolution in the area mostly covering Asia (10-50 N, 65-135E). The simulations are carried out using the Max Planck Institute Earth System Model for paleo application (MPI-ESM-P, Jungclaus et al 2014). The model consists of the spectral atmospheric model ECHAM6 (Stevens et al 2013) in T63L47 resolution, the ocean model MPIOM (Marsland et al 2003, Jungclaus 2013 in GR1.5L40 resolution and prescribed vegetation maps (Pongratz et al 2008). The model is run under PMIP3/CMIP5 protocols for the last millennium (850-1849, 'past1000') and the historical period (1850-2005, 'historical'). The runs 'past1000-r1' and 'historical-r1' over these periods were combined and are similar in length as most regional reconstructions from China.
The pseudo-reconstructed precipitation data are obtained by adding LTM to the pseudo-instrumental precipitation data. We use first Fourier-transform the time series, then change the scaling in the frequency domain S f f ( ) ∼ β , from the original value 0 β = to 1 β = , and then inverse-transform back into the time domain. Note that this procedure of adding long-term persistence only slightly decreases the correlation between time series. We also tested a fractional integrated statistical model and generated time series with similar properties (details not shown here).
Even though, proxy time series normally contain additional noise giving them a lower signal-to-noise ratio compared to instrumental data, we avoid this additional complication here, but focus on the effects of the modified LTM on TR-based reconstructions. Both pseudo data sets are normalized by subtracting the mean and dividing the standard deviation over the whole period . In this letter, extremes are defined here as annual values exceeding the upper and lower 5th percentiles of the time series.

Methodology
Considering the ability in dealing with non-stationary time series, and more accurate statistical outputs (Kim et al 2014), we apply the DFA method to diagnose LTM in all the data we use in this letter (Peng et al 1994, Kantelhardt et al 2001. Supposing we have N records x i , we first subtract the mean and calculate the accumulated sum (profile)  is least squares fitting of variable order polynomials. Normally, quadratic polynomial fitting is enough for the accurate estimation of LTM in temperature/ rainfall records, therefore in this study we choose to use the quadratic polynomial fitting. We can determine the 'detrended walk' as the difference between the original profile Y j v , and the local trend Y j v ,  and calculate the variance as while N s is the number of the windows, and v N 1, 2, 3, , s = ⋯ . If F DFA (s) follows by a power law, F s s ( ) DFA ∼ α . As stated in the introduction, an exponent 0.5 α > indicates a record x i is long-term correlated. While 0.5 α < indicates x i is long-term anti-correlated. 0.5 α = reveals white noise with no auto-correlation. Here, we use this exponent α to measure LTM, which is a generalization of the Hurst exponent and identifies long-time correlation, as well as the stationary/non-stationary nature of the data. DFA was applied to assess LTM in 15 temperature sensitive and 118 precipitation sensitive TRW records, as well as 98 MXD records (figure 1).

Results and Discussion
3.1. Long term persistence in single tree-ring chronologies To illustrate how the DFA exponent is measured, we provide an example of precipitation sensitive TRW records, 'Yang2014', in northeastern Tibet ( figure 1(a)). The fluctuation function shows a power law F n n ( ) ∼ α with 0.8 α = , ranging from 1 to more than 100 years. The DFA exponent of each TRW chronology is measured in this way indicating a mean of 0.86 α = and 0.82 α = from 1-50 years for the temperature and precipitation sensitive TRW time series, respectively. The scaling of some chronologies break at the time scale less than 50 years; this possibly depends on the length of tree core measurement series and the applied detrending methods (Briffa et al 1992, Cook et al 1995. The main climate drivers, temperature and precipitation, appear to have only limited influence on the TR proxy LTM as both land temperature and precipitation are characterized by weak-memory ( 0.6 α ≈ ) and no-memory processes ( 0.5 α ≈ ), respectively. The same analysis using with 98 MXD times series reveals a mean 0.7 α = for 1-50 years, consistent with the power spectrum exponent reported in Franke et al (2013). This indicates the existence of LTM in both TRW and MXD based reconstructions with TRW chronologies generally having higher DFA exponents compared to the MXD. One reason for lower DFA exponent values in the density data is that cell wall growth, as reflected by MXD, occurs over a shorter period in high and late summer, compared to TRW. TRW represent and integral of the cambial activity throughout the growing season, though weighted towards early summer, and is additionally influenced by storage effects of carbohydrates during previous years. These biological carry-over effects increase the autocorrelation in TRW timeseries (Frank et al 2007).
3.2. Long term persistence in regional mean climate reconstructions over China Except for climate reconstructions at single locations, intensified LTM is also detected in regional mean climate reconstructions integrating annually resolved TR data and decadally resolved multi-proxy data. The TRW based East-Asia mean annual PDSI reconstruction ('Cook2010' in figure 2) shows a nearly constant LTM 1.3 α = , on time scales between 1 year and 100 years. The reconstruction of annual mean temperatures over China ('SHI2012-EIV' in figure 2) shows two different scaling properties with 2.5 α = from 1-20 years, and 1.3 α = from 21-100 years. An explanation for the memory change is that high frequencies have been filtered out (using a 10-year low pass filter) in the initial data processing procedure (Shi et al 2012), which reddens the spectrum in the range 1-10 years. For the time range 20-100 years, the scaling exponent is consistent with the memory in 'Cook2010'.
For Yang2002, we found a DFA exponent 1.6 α = , in line with the power spectrum exponent 2 β = reported in Zhang et al (2011). Finally, 'Ge2013-PLS' has a scaling with 1 α = (figure 2). These findings reveal that the regional PDSI and temperature reconstructions over China contain similar scaling exponents a ranging from 1.0 to 1.6 over time scales from 1 year to 100-200years. Averaging single reconstructions over a large-scale area may increase a in the regional mean multi-proxy based temperature and PDSI reconstructions. Similarly, LTM of instrumental temperature over land is increased from 0.65 α = at local and regional scales to 0.75 α = at the global scale (Rypdal et al 2013). A possible explanation is that for these larger scale climate products the spatial correlation among time series might increase the correlation with external forcing and decrease the high-frequency signal from the individual/ local variability (Mann 2011).

Possible influences of intensified LTM on climate reconstruction
To demonstrate how the LTM can influence climate reconstructions, we perform an experiment to compare the no-memory pseudo-instrumental precipitation data with the with-memory pseudo-reconstructed precipitation in both time and spatial domain over southeastern Asia. The pseudo-proxy reconstructions are generated by artificially increasing low frequency variances of the pseudo-instrumental data at each grid point through the power spectrum transform as described in section 2.1.2. Figure 3 shows an example of comparisons between pseudo-reconstructed precipitation time series and pseudo-instrumental, both extracted from the grid at [105E, 23.3N]. The pseudo-instrumental precipitation (thin blue lines in figure 3(c)) has a power-spectrum exponent 0 β = and DFA exponent 0.5 α = ( figure 3(a)). In the pseudo-reconstructed precipitation both β and α equal 1. We smooth the data by using 20 yr moving average, as it is commonly used in reconstruction analysis to make the long-term fluctuations in the series stand out more clearly. Compared to the pseudoinstrumental precipitation, the pseudo-reconstructed time series contains several differing characteristics: the record at this specific grid point shows stronger low frequency variance and includes a generally wetter period 850-1500 followed by rather drier conditions 1500-2005. Wet/dry deviations appear more persistent with an increased chance that a wet year follows a wet year, and a dry year follows a dry year, a finding in agreement the assessment of precipitation reconstructions in Bunde et al (2013). Also, extreme events (the events passing the horizontal lines denoting the upper and lower 5th percentiles in figure 3(c)) more likely cluster in short periods as one main consequence of LTM (Bunde et al 2013), instead of being distributed randomly over time. Third, the pseudoreconstructed precipitation contains a decreasing trend ( decade 2% − ) according to linear regression analysis, compared to no significant trend in the pseudo-instrumental precipitation time series, as well as substantial trends in certain shorter periods, e.g. decade 8% compared to decade 1% − over the last 50 years. To show the influences of modified LTM on reconstructed climate in space, we use the pseudo data to calculate anomalies and linear trends over the last 105 years in southeastern Asia. Figure 4 shows the values of the anomalies and trends are intensified in most regions. For example, while the anomalies of pseudo-instrumental precipitation from 1901 to 2005, relative to the 850-2005 climatology, range from 60% − to 20% − in southern China, the anomalies of pseudo-reconstructed precipitation reach 100% <− . The pseudo-reconstructed precipitation shows much stronger trends compared to the pseudo-instrumental precipitation over most of the region. For example, the pseudo-reconstructed precipitation has a positive trend ranging from decade 20% 60% − in the Yangtze river region, but the pseudo-reconstructed precipitation shows a positive trend exceeding decade 80% in this region. These findings show that climate reconstructions based on proxy data containing LTM can lead to exaggerated climate anomalies and trends in certain regions and specific time periods.

Conclusion
In this paper, we demonstrate that TR chronologies contain more LTM than climate observations. DFA is applied to measure the strength of LTM in worldwidedistributed TR proxies and climate reconstructions (PDSI and temperature) from China. Our results reveal increased DFA exponent including mean value 0.8 α = for TRW proxies and 0.7 α = for maximum density data. For observational temperature and precipitation data, these values are smaller equaling 0.6 α ≈ and 0.5 α ≈ , respectively. We consider that tree physiological processes influencing TRW are biasing the LTM in tree-ring proxies. LTM in single tree-ring (and other) proxies not only propagates into regional mean climate reconstructions, but is even accelerated in these larger scale products.
One main characteristic of data containing LTM is that the successive increments and layers are positively correlated, causing several effects which requiring attention. By comparing the pseudo-instrumental and pseudo-reconstructed precipitation we demonstrated that modified LTM may stimulate (i) additional trend during certain intervals, as well as over the whole reconstruction period, while the trends from external forcing tangling with the trends from intensified climate low-frequency variations; In this case, a new approach is suggested by Lennartz and Bunde (2011) and   Tamazian et al (2015) to re-estimate the significance of trends in time series with LTM.
(ii) increased low-frequency variance, and overestimated mean state of climate anomalies during certain periods; (iii) unrealistic deviations and extreme events in certain regions in spatially resolved reconstructions.
To modify the LTM in TR data, models based on fractionally integral techniques could be used as LTM is in fact a fractal phenomenon, such as the autoregressive fractionally integrated moving-average model (ARFIMA) (Wei 1994) and the fractional integrated statistical model (Yuan et al 2014). Besides, since one TR chronology is typically developed by averaging data from multiple trees of limited lengths, first-order auto-regression or ARMA as a pre-whitening method (Meko 1981, Cook et al 1995 may be used during TR detrending to deduce long-term persistence. We here evaluate the influence of different detrending methods, which are commonly applied in dendrochronology (Cook et al 1995, Briffa et al 1992, and combine these approaches with an AR1 prewhitening procedure to assess LTM in a temperaturesensitive TR chronologies from Siberia covering past 2,000 years. Results show that SLP300 combined with an AR1 pre-whitening procedure could change the LTM to approximate values found in instrumental temperatures (Figures not shown here). However, this analysis needs to be applied to more millennium long TR chronologies to further test whether the combination of AR1 and 300-year cubic smoothing splines (SLP300) can be generally used to adjust LTM in temperature-sensitive TR data.
The processes generating LTM in tree-ring data are still not fully understood. Likely, it is influenced by soil moisture supply, as precipitation-sensitive trees rather response to changes in root-zone soil moisture.
Soil moisture has an integrative behavior, and is possibly long-term persistent in dry areas Fraedrich 2006, Wang et al 2010). Apart from that, it is also likely influenced by biological processes including carbohydrate storage and remobilization during cell wall construction and potential external disturbances and recovery of the forest structure (Frank et al 2007).