Long-Range Memory in Millennium-Long ESM and AOGCM Experiments

Northern Hemisphere (NH) temperature records from a reconstruction and a number of millennium-long climate model experiments are investigated for long-range memory (LRM). The models are two Earth system models and two atmospheric-ocean general circulation models. The periodogram, detrended fluctuation analysis and wavelet variance analysis are applied to examine scaling properties and to estimate a scaling exponent of the temperature records. A simple linear model for the climate response to external forcing is also applied to the reconstruction and the forced climate model runs, and then compared to unforced control runs to extract the LRM generated by internal dynamics of the climate system. With one exception the climate models show strong persistent scaling with power spectral densities of the form S(f)∼ f−β with 0.8< β < 1 on time scales from years to several centuries. This is somewhat stronger persistence than found in the reconstruction (β ≈ 0.7). The exception is the HadCM3 model, which exhibits β ≈ 0.6. We find no indication that LRM found in these model runs are induced by external forcing, which suggests that LRM on sub-decadal to century time scales in NH mean temperatures is a property of the internal dynamics of the climate system. Temperature records for a local site, Reykjanes Ridge, are also studied, showing that strong persistence is found also for local ocean temperature.


Introduction
The presence of long-range memory (LRM) in climatic records is well documented in the geophysics literature.LRM is characterized by an algebraically decaying autocorrelation function lim t→∞ C(t) ∝ t −γ such that ∞ 0 C(t)dt = ∞, i.e., 0 < γ ≤ 1. Equivalently, the power spectral density (PSD) of LRM time series follows a power law lim f →0 S(f ) ∝ f −β , where β = 1 − γ and 0 < β < 1.A typical model for an LRM stochastic process is the persistent fractional Gaussian noise (fGn).This is a stationary process with 0 < β < 1.The cumulative integral (or sum) of such a process has the PSD of the form S(f ) ∼ f −β , but with β → β +2.Such a process with 1 < β < 3 is a non-stationary LRM process called a fractional Brownian motion (fBm).
Because of the noisy nature of PSD estimators like the periodogram, other methods for estimating β are preferred (Beran et al., 2013).In this paper we use the periodogram as the first crude characterization of the data and for detection of spectral peaks indicating lack of power-law scaling, but employ Detrended Fluctuation Analysis (DFA) and Wavelet Variance Analysis (WVA) for parameter estimation.
Most of the LRM studies of climatic time series investigate local time records (e.g., Pelletier, 1997;Weber and Talkner, 2001;Eichner et al., 2003), but LRM has also been found in global observed temperature records (Lennartz and Bunde, 2009) and reconstructed temperature records for the Northern Hemisphere (Rybski et al., 2006;Mills, 2007).Some surface temperature records from AOGCM climate models have been analyzed with the main result that LRM is not reproduced in agreement with that of observational temperature (Syroka and Toumi, 2001;Bunde et al., 2001;Govindan et al., 2001Govindan et al., , 2002;;Bunde and Havlin, 2002).Some of the model experiments produce temperature with multiple scaling regimes, and some of them yield smaller scaling exponents than the observational temperature.However, in (Syroka and Toumi, 2001;Bunde et al., 2001;Govindan et al., 2001) the model experiments all had greenhouse gas forcing as the only dynamic forcing, while remaining external forcings, such as total solar irradiance and volcanic effects, were kept constant.Govindan et al. (2002) and Bunde and Havlin (2002) used experiments where (i) all forcings were fixed, (ii) with fixed forcings except greenhouse gas forcing, and (iii) with fixed forcings except greenhouse gas plus aerosol forcing.Their main conclusion was that the temperature from the model experiments fail to reproduce the scaling behaviour found in observational data, and that the models display large differences in scaling from different sites.Of these scenarios, the one with dynamic greenhouse gas plus aerosol forcing performed better with respect to producing the scaling observed in instrumental temperature records.Global fields of observed and simulated surface temperatures from an AOGCM climate model experiment were studied in Fraedrich and Blender (2003).The experiment was run with fixed forcings.The result from observational data was mostly in agreement with previous studies of temperature in oceanic and coastal regions, but the authors found white noise scaling (β ≈ 0) at continental interiors.Analysis of a 1000-year temperature simulation from the model experiment produced similar scaling exponents to what was found for the observational data in this study.Blender and Fraedrich (2003) made a similar analysis of temperature from two different model experiments with dynamic greenhouse gas forcing, giving results in agreement with Fraedrich and Blender (2003).
Temperature from model experiments with constant forcings, and time-varying greenhouse gas, sulfate aerosol, ozone, solar, volcanic forcing and various combinations was studied in Vyushin et al. (2004).Scaling exponents for temperature at 16 land sites and 16 sites in the Atlantic ocean were estimated.They found that inclusion of volcanic forcing considerably improved the scaling behavior.Rybski et al. (2008) used model experiments with all constant forcing and with dynamic solar, volcanic and greenhouse gas forcing.They analyzed data from grid cells all over the globe, but did not investigate global or hemispheric means.They found that for the forced run experiment the temperature showed a scaling exponent in agreement with observational temperature, while the temperature from the control run showed generally lower persistence.
Studies of LRM in temperature records from climate model experiments mostly use temperature from local sites, and some also use temperature spatially averaged over larger regions.Global mean temperature was studied by Syroka and Toumi (2001), but hemispheric means have not been studied with regards to LRM.For observational and reconstructed temperature, global and hemispheric means are also far less studied than local data.
In the present study we analyze scaling properties of surface temperature for the Northern Hemisphere from paleoclimate simulations and compare to those of temperature reconstruction by Moberg et al. (2005) which spans the last two millennia.Hemispheric temperature records from four different Earth system climate models are analyzed, and both forced runs and control runs are investigated.In order to avoid effects of anthropogenic forcing only data up to the year 1750 is used.This will give an idea of what role other natural external forcing like solar, CO 2 , volcanic and aerosol forcing play in producing LRM, and indicate if LRM can arise from internal dynamics alone.
Separation of the LRM arising from internal dynamics from the LRM induced by external forcing can also be achieved from reconstructed and simulated temperature data if the forcing data are known.The method makes use of a simple linear model for the global temperature response (Rypdal and Rypdal, 2013).The response to the external forcing can then be computed and subtracted from the observed or modeled temperature record to yield a residual which represents the internal variability of the climate system.Analysis of this residual and temperature from forced runs and control runs are compared for those models where temperatures from both forced runs and control runs are available.
This paper is organized as follows: Section 2 describes the DFA and WVA methods and the response model.Information about the models and the data used can be found in Section 3, and the results from the analysis are presented in section 4. Discussion and conclusion follow in section 5.

Detrended Fluctuation Analysis
The Detrended Fluctuation Analysis (DFA) (Peng et al., 1994;Kantelhardt et al., 2001) was explicitly designed to remove polynomial trends in the data.The method can be summarized in four steps.First, we construct the cumulative sum (the "profile") of the temperature time series x(t); Y (i) = i t=1 x(t) − x , where x denotes the mean.In the second step the profile is divided into N τ = N/τ nonoverlapping segments of equal length τ .This is done both starting at the beginning and at the end of the profile, so 2N τ segments are obtained altogether.In the third step, an n'th order polynomial is fitted to, and then subtracted from, each segment.Thus, at this stage we have formed the detrended profile Y τ (i) = Y (i) − p ν (i), where p ν (i) is the polynomial fitted to the ν'th segment.In the final step, the variance of each segment, is computed.The fluctuation function is given by the square root of the average over all the segments, .
The scaling parameter β is found through the relation What we have described is the n'th order detrended fluctuation analysis, denoted DFAn.It has the property of eliminating the effect of an n − 1'th order polynomial trend.In this paper we employ second order DFA, denoted DFA2, which eliminates linear trends.

Wavelet Variance Analysis
The continuous wavelet transform is the convolution between a time series x(t) and the rescaled wavelet Ψ(t/τ ); The mother wavelet Ψ(t) and all rescaled versions of it must fulfill the criteria ∞ −∞ Ψ(t ) dt = 0.For LRM time series, the variance F (τ ) = (1/N ) N t=1 W 2 (t, τ ) scales as a power-law (Flandrin, 1992;Malamud and Turcotte, 1999), The method is therefore known as the wavelet variance analysis (WVA).In this study we have used the Mexican hat wavelet, which is capable of eliminating linear trends, and denote the method WVA2.The properties of the WVA2 analysis are similar to the DFA2 in that it usually yields similar values of β.It is, however, much more sensitive to the presence of additional oscillations in the data, which show up as wavy structures in the function F (τ ).We use it in this paper mainly as a tool (in addition to the periodogram) to detect such oscillations.

The response model residual analysis
For the preindustrial period the most important contributions to the external radiative forcing F (t) are orbital, solar variability, and aerosols from volcanic eruptions.Orbital forcing can be computed with high accuracy, and total solar irradiation has been reconstructed for the last ten millennia.Existing reliable reconstructions of volcanic forcing cover the last millenium.The forcing data used here are further described in Section 3. The evolution of the global mean surface temperature anomaly T on decadal to millennial time scales can tentatively be modeled as a linear response to F (t) in addition to a response to stochastic forcing from unresolved spatiotemporal "turbulence" (e.g., forcing of the seasurface temperature from atmospheric weather systems).A simple stochastic-dynamic model (SDM) with an LRM response function is (Rypdal and Rypdal, 2013): Here B(s) is the Wiener stochastic process whose increments dB(s) is a Gaussian white noise process and σdB(s) represents the stochastic component of the forcing.T (t) is the temperature relative to the temperature T 0 at time t = 0 (the beginning of the record) and F = F + F 0 is perturbed forcing F relative to that of a radiative equilibrium at surface temperature T 0 plus the actual radiative imbalance F 0 at t = 0.By definition F (0) = 0. F 0 is a model parameter which is estimated from the data along with the other model parameters β, µ, and σ.The stochastic part of this solution (the term to the right) has a power spectral density of the form S(f ) ∼ f −β , and is fractional Gaussian noise (a stationary process) if β < 1 and a fractional Brownian motion (nonstationary) if 1 < β < 3. Time-series information about global climate forcing and its various components exists for the instrumental period as well as for the last millennium.This information can be used in conjunction with the observed temperature records to perform maximum-likelihood estimates (MLE) of the parameters of the model.The details of the MLE method applied to this response model are explained in Rypdal and Rypdal (2013).In a short-range memory response model, the powerlaw kernel (t − s) β/2−1 in the response model is replaced with an exponential e −(t−s)/τ , where τ is the time constant.In this case the parameter µ −1 can be interpreted as the effective heat capacity of the climate system.In the LRM response model µ −1 does not have a simple physical interpretation, although it is (in combination with β) a measure of the thermal inertia of the system.The memory parameter β estimated from this model should be interpreted as the LRM parameter for the internal temperature response, and hence the problem of separating the LRM contribution from the forcing and the internal LRM has been eliminated.The β estimated in Rypdal and Rypdal ( 2013) is β ≈ 0.75, which is not much lower than the value estimated for the full temperature record from detrending techniques like DFA and WVA.This shows that the detrending techniques effectively eliminate the contribution to β from the anthropogenic trend.
In the present paper the Crowley forcing (Crowley, 2000) is used for Moberg reconstructed temperature (Moberg et al., 2005) and for the temperature from the ECHO-G forced run experiment.The COSMOS experiment was run with a different forcing, and this forcing is used as input to our response model.For the temperature from the LOVECLIM experiment, solar and volcanic forcings were used together with forcings from CO 2 and tropospheric aerosols corresponding to the Crowley forcing.The full forcing data in these three cases are shown in Figure 1.
3 Data 3.1 The reconstruction of Moberg et al. (2005) The reconstructed temperature presented in Moberg et al. (2005) is a Northern Hemisphere reconstruction covering the time period 1-1979 AD.The reconstruction is created from 11 low-resolution proxy time series (e.g.ice cores and sediments, 1-180 year resolution) and 7 tree-ring records (annual resolution).The 18 local reconstructed temperature time series were first divided into an Eastern and a Western part.Linear interpolation was then applied to all time series in order to create annual mean values.The beginning and end of the time series were padded with surrogate data so that they all covered the time period 300 BC -2300 AD to minimize edge effects of the wavelet transform.The wavelet transform (WT) with the Mexican hat wavelet basis function was then applied using the set of 22 scales to generate 22 time series.For each scale 1-9 (Fourier timescales <80 years), the WT from the tree-ring proxy series were averaged.For the scales 10-22 (Fourier timescales >80 years), the WT from the lowresolution proxy series were averaged.Scale 1-22 were then merged, creating two full WT time series, one for the Eastern and one for the Western Northern Hemisphere.The two subsets were then averaged, and the inverse WT was calculated, creating a dimensionless NH temperature reconstruction.Finally, the mean and variance of the reconstructed temperature time series were calibrated to correspond to the instrumental data available for the time period 1856-1978.

Marine sediment SST reconstruction; Reykjanes Ridge
The local sea surface temperature (SST) reconstruction applied in the following study is presented in detail in Miettinen et al. (2012).Past August SST has been reconstructed by analyzing marine planktonic diatoms from a composite marine sediment core, recovered at the Reykjanes Ridge in the western subpolar North Atlantic, (57 • 27.09'N, 27 • 54.53'W, at 2630 m water depth).The composite core consist of a 54.3 cm long box core, and a 3.725 m long gravity core.The general assumption is that the down-core diatomic microfossil assemblages are related with past environmental conditions at the core site.Marine diatoms are unicellular, photosynthetic algae with siliceous frustules.For this particular analysis, the down-core diatomic assemblages were converted to August SST estimates by the weighted-average partial least squares technique (ter Braak and Juggins, 1993).The SST reconstruction has an average temporal resolution of 2 years for year 1770-2000 (box core), and 8-10 years for year 1000-1770 (gravity core).

SST reconstruction from observations; Reykjanes Ridge
A reconstruction based on instrumental observations was developed in Smith and Reynolds (2005).For the ocean, sea surface temperature (SST) was used, while surface marine air temperatures where left out due to biases in the daytime temperatures.The SST analysis and a separate land surface air temperature analysis were merged to form a monthly merged analysis from 1880 to 1997.The International Comprehensive Ocean-Atmosphere Data Set (ICOADS) SST observations release 2 was the primary SST data, but the combined satellite and in situ SST analysis of Reynolds et al. ( 2002) was also included.The reconstruction was separated into low-and high-frequency components, which were added for the total reconstruction.The low frequency was reconstructed using spatial and temporal filtering, with a time filter of 15 yr.The low-frequency component was subtracted from the data before reconstrucion of the high-frequency component using spatial covariance modes.The method for reconstructing the data is described in detail in Smith and Reynolds (2004).This reconstruction contains improvements over many earlier studies: It is globally complete, incorporates updates in ICOADS, the analysis variance have less dependence on sampling compared to some earlier analysis, and uncertainty estimates indicate when and where the analysis is most reliable.

LOVECLIM model and experiment
The Earth system model LOVECLIM version 1.2 contains a quasi-geostrophic model for the atmosphere (ECBilt2), coupled to an ocean GCM (CLIO3) (Goosse et al., 2010).The two models have 3 and 20 vertical levels, respectively.A thermodynamic sea-ice model is incorporated into the OGCM, and the vegetation model VECODE is used to simulate the dynamics of trees, grasses and deserts.It includes the evolution of the terrestrial carbon cycle, while a separate model LOCH simulates the ocean carbon cycle.Both the solubility and the biological pumps are included in this model.Incorporated in LOVECLIM is also the ice-sheet model AGISM, which consists of 3 modules; ice sheet flow, visco-elastic bedrock and mass balance at the ice-atmospehere and iceocean interfaces.We apply surface temperature data from one experiment with this model; "LOVECLIM Climate Model Simulation Constrained by Mann et al. 2009 Reconstruction" (Goosse et al., 2012).In this experiment, simulations are constrained by the mean surface temperature reconstruction of Mann et al. (2009).External forcing includes TSI (total solar irradiance), volcanic eruptions, land cover changes, orbital forcing, greenhouse gases and aerosols.When we implement the response model to these data, only time series for the solar, volcanic and greenhouse gas forcing are applied.The solar forcing time series is based on the reconstruction by Muscheler et al. (2007).The volcanic activity time series originate from Crowley et al. (2003) , while the greenhouse gas forcing used is obtained from (Crowley, 2000).

COSMOS ESM model and experiments
The COSMOS ESM model consists of GCMs for the atmosphere and the ocean (Jungclaus et al., 2010).The atmospheric model ECHAM5 (Roeckner et al., 2003) has 19 vertical levels, while the ocean model MPIOM (Marsland et al., 2003) has 40.A thermodynamic sea-ice model is incorporated into the OGCM.Additional modules include the ocean biogeochemistry model HAMOCC5 (Wetzel et al., 2006), and the terrestrial biosphere model JSBACH (Raddatz et al., 2007).
The surface temperature data applied in our analysis are extracted from one experiment in a set of experiments referred to as "Ensemble Simulation of the Last Millenium using the Comprehensive COSMOS Earth System Model" (Jungclaus et al., 2010).External forcing used in the forced simulations include TSI, volcanoes, orbital forcing, greenhouse gases and land use change.An unforced control run is also used here in the comparative LRM study.
For the response model, time series for solar, volcanic and greenhouse gas forcing are applied.The forcing time series used are created specifically for this model and experiment (Jungclaus et al., 2010).The solar forcing time series is based on a combination of reconstructions; from the Maunder Minimum (1647-1715 AD) until today the total solar irradiance (TSI) is based on historical sunspot records (Krivova and Solanki, 2007;Balmaceda et al., 2007), and between 800 AD and the Maunder Minimum the TSI is reconstructed from estimates of the solar open magnetic flux based on 14 C con-centrations in tree rings (Solanki et al., 2004;Krivova and Solanki, 2008;Usoskin et al., 2011).An 11-year solar cycle has been superposed on this part of the reconstruction.
The relative radiative forcing from volcanic eruptions is calculated from aerosol optical depth (AOD) and effective radius R eff .Satellite data from the 1991 Mt. Pinatubo eruption is the basis for these estimates.The greenhouse gas forcing includes CO 2 , where concentrations are computed within the model, based on historical records of fossil fuel emissions by Marland et al. (2003).

ECHO-G model and experiments
The coupled model ECHO-G (Legutke and Voss, 1999) version 4 consist of GCMs for the ocean/sea ice and the atmosphere.The atmospheric model ECHAM4 (Roeckner et al., 1996) includes 19 vertical levels, while the ocean model HOPE-G (Legutke and Maier-Reime, 1999) includes 20 levels.External forcing includes volcanoes, solar irradiance and greenhouse gases, all derived from Crowley (2000) and is used in the response model study.Surface temperature from two experiments is used for analysis; one forced run and one control run with forcing values fixed to year 1990 (Zorita et al., 2003;González-Rouco et al., 2003;von Storch et al., 2004).

HadCM3 model and experiment
The Hadley Centre coupled model 3 is an AOGCM (Gordon et al., 2000), with 19 levels in the atmospheric component HADAM3 (Pope et al., 2000) and 20 levels on the ocean HADOM3 component.External forcing is constant for this experiment (Collins et al., 2000).
The complexity, time period coverered, and temporal and spatial resolution for each model experiment are given in Table 1.

Results
When applying DFA2 and WVA2 directly to the temperature reconstruction records, all data up to the year 1750 are used.Because the Crowley forcing record starts at 1000 AD, only data from this year and forward were used in the re-Table 2. Estimated β from applying DFA2 and WVA2 directly to the full temperature record (all temperatures), from DFA2 and WVA2 applied to the residuals from the deterministic response, and from the response model using MLE (temperature reconstruction and temperature from the forced climate model run experiments).The scaling ranges for DFA2 and WVA2 are also shown in years for Moberg and LOVECLIM, and in months otherwise.The same scaling range has been used for the full record and the residual, except for WVA2 applied to the forced ECHO-G experiment, where the upper scale used for the residual is in parenthesis.sponse model residual analysis of the Moberg reconstructed temperature and the temperature from the LOVECLIM experiment.Therefore the scales shown in the plots for DFA2 and WVA2 may differ somewhat between the full record and the residual from the deterministic response.Table 2 shows the resulting β from applying DFA2 and WVA2 directly to the full temperature and to the residuals from the deterministic response.The β estimated using the response model is also given, where the parameters of the response model are estimated with the MLE method (Rypdal and Rypdal, 2013).The response model residual analysis is applied to the temperature reconstructions and the temperature from forced climate model experiments, while the direct analysis also include temperature from control runs.The scaling ranges used to find β with DFA2 and WVA2 are also indicated in this table.Note that the scaling range is given in years for the Moberg and LOVECLIM temperatures, and in months otherwise.In principle the LRM properties due to internal dynamics can be separated from those induced by the external forcing by applying the response model method of Rypdal and Rypdal (2013) described in section 2. This method allows estimation of the model parameters β, σ, and µ from the Crowley forcing data and the Moberg reconstruction record.Then we can compute the deterministic response and the residual obtained by subtracting this deterministic response from the Moberg record.The residual represents the response to the stochastic forcing and hence the internal variability of the climate system.The scaling properties of this residual can be assessed with the DFA method which also provides a consistency test on the maximum-likelihood estimate of β.

Full record
A caveat here is the low-pass filtered nature of the Moberg record.The MLE method tends to emphasize the shorter scales on which the reconstruction record is smooth, and this will spuriously yield β ≈ 1.A way to avoid this could be to coarse grain both temperature and forcing time series by averaging over successive time windows of a certain length t A , such that the temperature series is no longer smooth.This will give a more reasonable maximum-likelihood estimate of β, but the coarse-grained data cannot capture the causal connection between the almost instantaneous volcanic forcing spikes and the temperature response to them.The resulting blurring of the causal connection on time scales shorter than a decade has the effect that the MLE method interprets the variability on these short scales as stochastic, and hence overestimates the stochastic forcing strength σ, and also yields incorrect estimates of µ and β.The lesson to learn from this is that we cannot expect to obtain a correct separation of deterministic and stochastic forcing and correct parameter estimates from the low-resolution reconstruction data.Another approach to circumvent this problem was suggested in Rypdal and Rypdal (2013), where the response model parameters computed from the annual-resolution instrumental data were applied to the millennium-long annual-resolution forcing record to produce a deterministic-response record with annual resolution.The Moberg record and this deterministic response is shown in Figure 3(a).The residual obtained by subtracting the deterministic response from the reconstructed record is shown in Figure 3(b), and provides a good representation of the internal variability on time-scales longer than a decade.On shorter time-scales the residual is strongly influenced by the forced response due to the smooth character of the temperature reconstruction, but we do not need to use those scales to estimate model parameters if we do not insist on using MLE.The PSD of the residual is shown to exhibit good scaling in Figure 3(c), and the DFA2 fluctuation function of this residual on the longer time-scales should provide good estimates of β for the internal variability, as shown in Figure 3(d).

Results from LOVECLIM experiment
The NH temperature record for the period 1000-1750 AD for the LOVECLIM experiment, its power spectral density (PSD), and the DFA2 and WVA2 fluctuation functions are shown in Figure 4. DFA2 and WVA2 show good scaling with β ≈ 0.97 at least on time scales up to a few hundred yr.The response model gives similar value of β, which suggests that the persistence observed in the modeled record is  due to LRM in the internal dynamics and not a reflection of LRM in the forcing.In this model simulation both forcing input and simulation output are given with annual resolution.This allows us to handle volcanic forcing and the response to volcanic eruptions in a realistic manner.The results from the response model estimates with annual resolution are shown in Figure 5.

Results from COSMOS experiment
The temperature from the COSMOS forced run experiments exhibits some oscillations.In particular a prominent peak corresponding to a multiannual mode is seen in the PSD in Figure 6

Results from ECHO-G experiments
The temperature from the forced experiment "Erik1" shows good scaling with β ≈ 0.9 in Figure 9.The temperature from the control run (Figure 10) also scales well with a similar value for β.The response model yields a slightly smaller β (Figure 11).Here a 1-yr coarse graining has been applied before the parameters have been estimated, since the forcing data have 1-yr resolution.

Results from the HadCM3 experiment
The HadCM3 experiment consists of only a control run, and the scaling properties of the NH temperature series from this experiment differs from the other experiments in some important respects.Figure 12 shows a marked cross-over between two scaling regimes for τ ∼ 100 months in the DFA2 plot and τ ∼ 45 months in the WVA2 plot.The two regimes correspond to fBm-scaling with β ≈ 2 (Brownian motion) for the smaller scales, and fGn-scaling with β ≈ 0.6 for the larger scales.Rather than being dominated by an ENSO-like quasi-oscillatory mode up to scales of a few years, as observed in the COSMOS experiments, we observe here a nonstationary random-walk-like process for those scales.On the longer time scales this model also exhibits persistent scaling, but with somewhat lower persistence than observed in data from control runs and response model estimates in the other climate model experiments we have investigated.

Scaling in local data; Reykjanes Ridge
Analysis of instrumental local station data from continental interiors typically yields very low persistence on time scales up to a few decades.On the other hand, our experience is that coastal and oceanic observations in the temperate regions of the Northern Hemisphere present persistent β-values closer to those found for the hemispheric average.We also believe that this feature extends beyond the decadal scales, i.e., that good scaling with strong persistence prevails on scales up to several centuries in the Northern oceans.As an illustration we present in Figure 13 analysis of the Reykjanes ridge reconstruction from marine sediments described in Section 3.2 and the ECHO-G Erik1 experiment for the period 1000-1750 AD, and of the monthly SST reconstruction for the period 1880-1997 as described in Section 3.3.The figure shows DFA2 plots for the three data sets.ECHO-G shows good scaling in the range 1-100 yr with β ≈ 0.67, as compared to β = 0.91 for the NH-average.The marine sediment reconstruction yields β ≈ 0.45, and the instrumentally-based reconstruction β ≈ 0.56.The greatest uncertainty in the βestimate is in the marine sediment reconstruction, for which a very limited range of scales is available for analysis.The record has uneven time spacing, but the time step is mostly almost the same, slightly below 10 years.The data points inconsistent with this are ignored, and DFA2 applied to the remaining record.A maximum-likelihood estimation method for time series with uneven time spacing yields β close to what was found with DFA2.In spite of these uncertainties the analysis demonstrates consistently persistent scaling over time scales from years to centuries in these local data from model experiment and reconstruction data.

Conclusions
The temperatures from all model experiments except HadCM3 a yield higher values of β than the Moberg reconstruction when scales longer than a decade are considered.The reconstruction is said to represent temperature in the Northern Hemisphere, but most of the proxies used are in land and coastal areas.In this sense they may be considered more like representations of land or coastal temperature.Studies of observational data show higher persistence in sea surface temperature than land air temperature, and the value for β found for the Moberg record is more in agreement with the one found for the Northern Hemisphere land temperature than ocean temperature (Eichner et al., 2003;Lennartz and Bunde, 2009).Our estimate is in agreement with Rybski et al. (2006).The temperature from the model experiments is averaged over grid cells from both land and ocean areas, and the influence of the ocean might be what yields the higher value of β than found for the Moberg reconstruction.
The temperatures from the COSMOS experiments (both forced and control run) clearly show an influence of a quasiperiodic variability with a 2-3 year period, which can be associated with ENSO.The ECHO-G experiments show less influence of this oscillation, in HadCM3 it appears more like a crossover between two scaling regimes, and in LOVE-CLIM it is not noticeable.For the reconstructed temperature the ENSO time scales are not resolved, but in instrumental records the tropical oceans show a spectrum similar to that found from HadCM3.
The temperatures from the forced ECHO-G experiment and the LOVECLIM experiment show a more distinct Little Ice Age anomaly, in agreement with the temperature reconstruction, than the temperature from the COSMOS forced run experiment.This anomaly may also influence the estimation of β.
All the NH-averaged temperatures from forced experiments show clear persistent scaling with 0.8 < β < 1 on most of the available scales, i.e., from a decade to several centuries.The control runs and the response model estimates from the forced runs, which reveal the memory properties of the internal climate dynamics, do not show systematically less persistence than obtained directly from the simulated forced temperature records.This observation does not support the notion that the observed long-range memory to great extent is generated by the forcing.Such a suggestion was made by Rybski et al. (2008), based on a global map for the parameter α = (β + 1)/2 computed from both forced runs and control runs of the ECHO-G model.We believe that this discrepancy is caused by the reduction of spatiotemporal noise implied in performing an NH-average.The differences in estimated α between forced and unforced runs for local data may not be in the persistence of the underlying global signal, but rather in differences related to the amplitudes of spatiotemporal modes for the two types of runs.

LFig. 1 .
Fig. 1.The different forcings used as input to the response model, i.e.Crowley forcing used with the Moberg reconstruction (black), forcing used in the COSMOS experiment (red) and forcing used in the LOVECLIM experiment (green).
Figures 2-12 show the analysis of the temperature records.For the full records, the figures show (a) the temperature data, and (b) PSD, (c) DFA2, and (d) WVA2 applied to the data set.For the response model results, the figures show (a) the temperature data and deterministic response, (b) the residual, and (c) PSD and (d) DFA2 applied to the residual.The residual is the deterministic response interpolated to have the original time resolution subtracted from the full-resolution temperature data.For DFA2 and WVA2 95% confidence areas are shown, computed from Monte Carlo ensembles of synthetic fGns.For the response model analysis with DFA2, they are given the β estimated by MLE, otherwise they are given β corresponding to the one found from the DFA2 and WVA2 analysis respectively.The values of β used are indicated in each figure.

Fig. 2 .
Fig. 2. (a) The Moberg reconstructed temperature.(b) PSD, (c) DFA2 and (d) WVA2 applied to the reocrd.The blue areas are the 95% confidence for synthetic fGn with β estimated with DFA2 and WVA2, indicated in the figure.

4. 1 Fig. 3 .
Figure2shows the Moberg reconstructed temperature record, its power spectral density (PSD) and the results of DFA2 and WVA2 applied to the full record.As discussed above the scales up to ∼ 10 years are not representative for the temperature, and this is seen as a cross-over in the slope of both DFA2 and WVA2 fluctuation functions.The deviation from a straight line at the largest scales (lowest frequencies), which is most prominent in the WVA2 fluctuation function, is caused by a nonlinear trend associated with the

Fig. 4 .
Fig. 4. (a) The temperature from the LOVECLIM experiment.(b) PSD, (c) DFA2 and (d) WVA2 applied to the record.The blue areas are the 95% confidence for synthetic fGn with β estimated with DFA2 and WVA2, indicated in the figure.

Fig. 5 .
Fig. 5. (a) Temperature from the LOVECLIM experiment 1year resolution (black) and deterministic response (red).(b) The residual from the deterministic response.(c) PSD and (d) DFA2 applied to the residual.The red area is the 95% confidence for synthetic fGn with β estimated with MLE using the response model, indicated in the figure.Estimating β from the residual with DFA2 yields β = 1.01.

Fig. 6 .
Fig. 6.(a) The temperature from the COSMOS experiment.(b) PSD, (c) DFA2 and (d) WVA2 applied to the record.The blue areas are the 95% confidence for synthetic fGn with β estimated with DFA2 and WVA2, indicated in the figure.

Fig. 7 .
Fig. 7. (a) The temperature from the COSMOS controlrun.(b) PSD, (c) DFA2 and (d) WVA2 applied to the record.The blue areas are the 95% confidence for synthetic fGn with β estimated with DFA2 and WVA2, indicated in the figure.

Fig. 8 .
Fig. 8. (a) 4-year average of the temperature from the COSMOS experiment (black) and deterministic response (red).(b) The residual from the deterministic response.(c) PSD and (d) DFA2 applied to the residual.The red area is the 95% confidence for synthetic fGn with β estimated with MLE using the response model, indicated in the figure.Estimating β from the residual with DFA2 yields β = 0.772.

Fig. 9 .
Fig. 9. (a) The temperature from the ECHO-G experiment Erik1.(b) PSD, (c) DFA2 and (d) WVA2 applied to the record.The blue areas are the 95% confidence for synthetic fGn with β estimated with DFA2 and WVA2, indicated in the figure.

Fig. 10 .
Fig. 10.(a) The temperature from the ECHO-G controlrun.(b) PSD, (c) and (d) WVA2 to the record.The blue areas are the 95% confidence for synthetic fGn with β estimated with DFA2 and WVA2, indicated in the figure.

Fig. 11 .
Fig. 11.(a) 1-year average of the temperature from the experiment Erik1 (black) and deterministic response (red).(b) The residual from the deterministic response.(c) PSD and (d) DFA2 applied to the residual.The red area is the 95% confidence for synthetic fGn with β estimated with MLE using the response model, indicated in the figure.Estimating β from the residual with DFA2 yields β = 0.720.

Fig. 12 .
Fig. 12.(a) The temperature from the HadCM3 experiment.(b) PSD, (c) DFA2 and (d) WVA2 applied to the record.Two scaling regimes are found with β estimated with DFA2 and WVA2, indicated in the figure.

Fig. 13 .
Fig. 13.DFA2 applied to sea surface or air surface temperature at Reykjanes Ridge.The upper curve is the result for the air surface temperature from Erik1 experiment.The lower left curve is based on the monthly reconstructed data for sea surface temperature, and the lower right curve on the marine sediment reconstruction of ocean temperature.The red lines indicate the scales used to estimate β.

Table 1 .
Information on temperature from model experiments