Probing the timescale-dependent stability of local and global surface air temperature from climate simulations and reconstructions of the last millennia
PProbing the timescale-dependent stability of local and global surface air temperaturefrom climate simulations and reconstructions of the last millennia
Beatrice Ellerhoff ∗ and Kira Rehfeld Institute of Environmental Physics, Ruprecht-Karls-Universit¨at Heidelberg, INF 229, 69120 Heidelberg, Germany (Dated: February 2, 2021)Earth’s climate can be understood as a dynamical system that changes due to external forcingand internal couplings. Essential climate variables, such as surface air temperature, describe thisdynamics. Our current interglacial, the Holocene (11,700 years ago to today), has been character-ized by great stability of global mean temperature prior to anthropogenic warming. However, themechanisms and spatiotemporal patterns of fluctuations around this mean, called temperature vari-ability, are poorly understood despite their socio-economic relevance for climate change mitigationand adaptation. Here, we examine discrepancies between temperature variability from model sim-ulations and paleoclimate reconstructions by clarifying the stability of local and global surface airtemperature on the timescale of years to centuries. To this end, we contrast power spectral densities(PSD) and their power-law scaling using simulated and observation-based temperature series of thelast six thousand years. We further introduce the spectral gain to disentangle the externally-forcedand internally-generated variability as a function of timescale. It is based on our estimate of thejoint PSD of radiative forcing, which exhibits a scale break around the period of ten years. We findthat local temperature series from paleoclimate reconstructions tend to show unstable behavior onperiods of 10 to 200 years, while simulated temperatures almost exclusively show stable behavior.Conversely, the PSD and spectral gain of global mean temperature are consistent across data sets.Our results point to the limitation of climate models to fully represent local temperature statisticsover decades to centuries. By highlighting the key characteristics of temperature variability, wepave a way to better constrain possible changes in temperature variability with global warming andassess future climate risks.
I. INTRODUCTION
The variability of surface air temperature is presenton all spatial and temporal scales, from synoptic andseasonal changes to long-term variations on periods ofyears to multi-millennia. On the one hand, it arises frominternal processes, such as the El Ni˜no-Southern Oscil-lation (ENSO) [1]. On the other hand, the temperaturevaries due to external forcing, such as the greenhouse ef-fect [2, 3]. Understanding the internally-generated andexternally-forced variability has been suggested to be atleast as necessary for evaluating climate risks for societyand ecosystems as projecting the global mean tempera-ture [4]. Available instrumental observations are limitedto a small time span, leading to challenges in quantifyingtemperature variability. Paleoclimate reconstructions ex-tend the characterization of temperature variability andcan be compared to global circulation models (GCMs)[5–7]. However, discrepancies between model and pale-oclimate data remain to be resolved, especially on thelocal level and on periods between years and centuries[8–12].Characterizing local temperature variability is crucialfor predicting extremes [6], not only to minimize short-term damage but also to design long-term strategies, in-cluding urban planning and food cultivation [13]. Vari-ability of global temperature on periods above years isrelevant to the understanding of long-term changes [14] ∗ beatrice.ellerhoff@iup.uni-heidelberg.de as well as climate sensitivity [15]. It could affect theconfidence in future projections and attribution studies[16, 17]. Therefore, one of the main topics to be investi-gated here is the characteristics of local and global tem-perature variability on periods of years to centuries frommodel simulations and observation-based data of the lastmillennia.To determine how the variability of a temperature se-ries is distributed with timescales τ , we make use ofthe power spectral density (PSD) S ( τ ), known as “spec-trum”. It can be obtained from the Fourier transform ofthe autocorrelation function (Appendix A) [18, 19] andwas shown to often follow a power law S ( τ ) ∼ τ β , (1)with spectral exponent β and period τ [20–26]. We referto this behavior as temporal scaling since the tempera-ture signal has no preferred timescale and is statisticallysimilar within a power-law scaling regime τ ∈ [ τ , τ ].The spectral exponent β , also called scaling coefficient,quantifies the long-range dependence and temporal per-sistence of the process.Long-range memory stochastic processes are suitableto describe temperature signals with temporal scaling[25, 27]. Among those, fractional Gaussian (“white”)noise exhibits a spectral exponent β ∈ ( − ,
1) and astationary variance. Conversely, fractional Brownian(“red”) noise shows β ∈ (1 , σ ( τ ) ∼ τ β − .We summarize this behavior under the term “stability”:A stable temperature signal shows a stationary varianceacross timescales. In the non-stationary, unstable case, a r X i v : . [ phy s i c s . a o - ph ] J a n FIG. 1. Characteristic timescales relevant to surface air tem-perature (T) variability of climatic drivers (dark blue) andclimate subsystems (yellow) [28–30]. The weather and long-term climate shows a non-stationary variance of local andglobal mean temperature (“T unstable”). On interannual tomillennial timescales stability remains to be determined, es-pecially at the local scale. the temperature fluctuates stronger on increasingly longperiods [24].Unstable behavior ( β ≈
2) [22, 24] is typical for theweather regime (hours to weeks) and can be explainedby atmospheric turbulence [31, 32]. In the long-term cli-mate, regional and global mean temperatures show sim-ilar behavior ( β >
1) [23, 24, 33] due to the presence ofnonlinear processes, such as the temperature-albedo feed-back [34]. On timescales between years to millennia, thetemperature is constantly influenced by the interactionof all climate subsystems and by volcanic, solar, as wellas CO forcing (Figure 1). Estimates of the spatially-dependent scaling behavior of local temperature on thesetimescales differ [21, 24]. On the global scale, many stud-ies find β ≈ et al. has iden-tified a change from the so-called macroweather regime( β ≈ . β ≈ . II. DATA
We investigate the timescale-dependent distribution ofsurface air temperature variability using model simula-tions, observation-based data, and radiative forcing re-constructions. The model simulations include ten tran-sient runs from GCM experiments [42]. The observation-based data consists of reanalysis data, instrumental mea-surements, and the paleoclimate reconstructions from thePast Global Changes 2k (PAGES2k) network [43]. Weuse twelve reconstructions of climatic drivers, includingsolar, volcanic, orbital, and CO forcing. All tempera-ture and radiative forcing signals are specific to the Mid-and Late-Holocene (the last 6000 years), with a focus onthe Common Era (0 to 2000 CE). The supplementary (S)Tables S1-3 [44] summarizes their key specifications. A. Model simulations
Each of the ten GCM runs considered features a tran-sient, albeit different forcing and a comparable spa-tiotemporal resolution. The CESM-LME 1 [45] and MPI-M LM [46] experiments serve as representative runs ofthe last millennium. We analyze them at two temporalresolutions (one month, six hours) to capture both thehigh- and low-frequency variability within our availablecomputing capacities (Figure S7 [44]). CESM 1 past 2k[47] is included as a slightly newer run for the CommonEra. To cover the Mid-Holocene, we use simulations fromthe IPSL [48] (denoted IPSL-p6k) and ECHAM5/MPI-OM [49] (denoted ECH5/MPIOM-p6k) of the last 6000years. From the TraCE-21k [50] simulation, we also con-sider only the last 6000 years to retain comparabilityand to avoid spectral artifacts due to significant shifts inthe mean state of climate. The Mid-Holocene runs weretemporally averaged to a bi-monthly resolution to reducecomputational costs. To test for the influence of human-induced climate change on our results, we include theHadCM3 LM1 simulation [12], covering the period from850 to 1850 CE. Furthermore, we compare our resultsto the pre-industrial (PI) control runs from CESM-LME1 and MPI-M LM, as well as the TraCE-21k-ORB run,which is solely forced by orbital changes.
B. Observation-based data
In addition to the simulations, we analyze the monthlyresolved HadCRUT4 (Hadley Centre/Climatic ResearchUnit Temperature) instrumental records, ranging from1850 to 2019 [51]. However, most of the grid-box timeseries are not available as continuous measurements asrequired for spectral analysis. Therefore, we retain onlythose 104 grid boxes with coverage greater than 150 yearsafter interpolating gaps of up to two months. While theNorthern Hemisphere is comparatively well covered upto 72.5 ◦ N, only nine grid boxes remain for the South-ern Hemisphere. Therefore, this selection comes at theexpense of spatial resolution but offers a higher spectralresolution on longer timescales. To further explore thepotential effect of these spatiotemporal constraints, weinclude the ERA5 (European Centre for Medium-RangeWeather Forecasts Reanalysis 5th generation) tempera-ture reanalysis for the years 1979 to 2019 [52]. Along withCESM-LME 1 and MPI-M LM, we analyze the ERA5data at both six-hourly and monthly resolution for thesame reasons as mentioned earlier.In addition to direct temperature observations and re-analysis, we analyze paleoclimate data. Paleoclimaterecords hold preserved biological, chemical, and physi-cal tracers (“proxies”) of past climate. The number oftemperature records from paleoclimate data with sub-centennial resolution is limited. Recent progress has beenmade by improved calibration and pseudo-proxy methodswithin the PAGES2k network [53]. Therefore, we baseour analysis on their newest global multiproxy databasefor temperature reconstructions of the Common Era [43].It gathers 692 records from trees, ice, sediment, corals,speleothems, and documentary evidence with a resolu-tion between weeks and centuries. The records are spreadover 648 locations, including all continental regions andmajor ocean basins.For investigating the variability of global mean surfacetemperature, we use the seven spatially-weighted statis-tical reconstructions for the last 2000 years provided byPAGES2k [43]. To estimate the mean of local spectra,we choose records from the PAGES2k database accord-ing to their resolution ( ≤ ≥ ≥
20 years), as well astheir maximum hiatuses ( ≤
160 years). To reliably de-duce the scaling of the PSD from individual records, weselect the records according to our scales of interest (Ta-ble I), similar to [25, 54]. Ice core records were excludedfrom our analysis since they require additional considera-tion of spectral artifacts and signal-to-noise ratios at thesub-centennial timescales [55, 56].
TABLE I. Requirements on irregularly sampled time series x ( t ) for analyzing power-law scaling on timescales τ ∈ [ τ , τ ].We apply this scheme for τ = 10 and τ = 200 years insection IV B and IV C.Parameter ValueNumber of data points ( N ) ≥ h t i +1 - t i i ) ≤ τ Coverage ( t N - t ) ≥ τ Length of hiatuses (max( t i +1 - t i )) ≤ τ C. Radiative forcing
External forcing contributes significantly to tempera-ture variability and is an essential part of reliable climateprojections [37, 57, 58]. We study its spectral proper-ties using forcing reconstructions, widely implemented inGCM experiments and coordinated within the Palaeocli-mate Model Intercomparison Project (PMIP3/PMIP4)[59, 60]. This includes five solar [61–65], one CO [60]and two volcanic [57, 66] forcing reconstructions as well asBerger’s numerical solution for orbital forcing [67]. Fur-thermore, we calculate diurnal insolation changes fromthe hour angle of the sun [68]. We also use a more re-cently published volcanic [69] and high-resolution solarforcing [70] reconstruction as well as CO measurements[71]. We neglect land-use forcing [72] which is much lowerin amplitude and variability than the other forcings con-sidered here.All forcing reconstructions are rescaled to radiativeforcing equivalents, which express their respective changein the Earth’s radiation balance in Watts per squaremeter (Wm − ). We apply the widely-used formula5 .
35 ln([CO ] / − to rescale CO concentra-tions [CO ], given in parts per million (ppm) [73]. Thestratospheric aerosol optical depth (AOD) from volcaniceruptions is rescaled by ( − − Wm − /AOD [74], how-ever, the optimal conversion factor is still a matter ofdebate [75]. Additional uncertainties arise from the widespread of reconstructions for volcanic and solar forcing.To account for this and the choice of conversion factor,we simulate the joint PSD of radiative forcing by a MonteCarlo approach described in Appendix D. Here, “joint”indicates that the PSD of radiative forcing is calculatedby linear summation of the mean PSD from differenttypes of climatic drivers, rescaled to their radiative forc-ing equivalents. III. METHODS
Spectral analysis is the primary tool used here forstudying the scaling and stability of temperature se-ries. To minimize uncertainties in the spectral analysisof proxy records, we use state-of-the-art approaches forirregularly sampled time series [76]. Statistical estima-tors further test for the agreement between simulationsand paleoclimate data. We apply linear response the-ory to derive the spectral gain and investigate the forcedtemperature response.
A. Spectral analysis
Power spectral analysis requires the assumption thatthe underlying time series can be described as a weakly-stationary, stochastic process [77]. This is reasonablyfulfilled for our temperature time series, given that welinearly detrend them by default. The latter is also suffi-cient to remove potential biases from trends, such as therecent global warming (see supplementary (S) Figures S9and S10 [44]). We use the multitaper method with threewindows [78, 79] and chi-square distributed uncertaintiesto compute the PSD. The two lowest frequencies wereomitted as it is standard for reducing biases of the mul-titaper method [23]. For visual purposes, we apply alogarithmic Gaussian smoothing filter of constant width(0.005 decibels) [80]. Mean spectra were calculated by in-terpolation to the lowest resolution, binning into equallyspaced log-frequency intervals, and taking the averagewith equal weights [23]. This requires the statistical in-dependence of the averaged values [40]. The spectralexponent β is calculated by linear regression to the log-arithm of (1) on periods between τ = 10 and τ = 200years. In the case of seven proxy records with an insuf-ficient resolution, the scaling is estimated on their corre-sponding spectral resolution, but always at least between20 and 200 years (Figure S4 [44]).The uncertainty of thespectral exponent, ∆ β , is given by the standard error ofthe linear regression model ∆ β lm , except for irregularlytemperature series. B. Uncertainties for irregular temperature series
Spectral analysis of proxy records, which are typicallynot sampled in regular time steps, is more prone to er-rors than that of regular time series. We aim to mini-mize biases by accounting for the number of data points,temporal resolution, total coverage, and hiatuses’ lengthwhen selecting the records (Table I). We find that themean temporal resolution of a proxy record approximateswell the optimal interpolation time step. Nevertheless,the interpolation introduces uncertainties which are notcaptured by ∆ β lm . Similar to Laepple et al. [76], wequantify this additional uncertainty ∆ β int in four steps:(1) For each record with spectral exponent β , we simulate N = 100 surrogate time series with annual resolution anda power-law scaling β n ≈ β and n ∈ [1 , N ]. (2) We formthe surrogate’s block-average over the proxy record’s ir-regular time steps and obtain N surrogate time series atrecord resolution. (3) We interpolate the surrogate timeseries, calculate the multitaper spectrum, and extract thescaling exponent β n,lm from linear regression in the same way as for the proxy record (Figure S8 [44]). (4) We cal-culate the mean deviation ∆ β int = N P Nn =1 | β n,lm − β n | of the ensemble. The uncertainty of the individual fits∆ β n,lm is negligible compared to the mean deviation∆ β int . We obtain the uncertainty of the record’s spec-tral exponent from both, the uncertainty of the initialfit ∆ β lm and due to interpolation ∆ β int via quadraticsummation: ∆ β = p (∆ β lm ) + (∆ β int ) . C. Statistical analysis of spectral exponents
We quantify the agreement of simulated and recon-structed β -values using percent agreement, categoricalagreement, and Kappa statistics. Beforehand, we ex-tract the simulated temperature at the proxy record lo-cation by bilinear interpolation of neighboring grid boxesto achieve the best possible comparability between recordand simulation. Percent agreement p gives the percent-age of locations at which the confidence range β ± ∆ β from simulation and reconstruction overlap. The agree-ment by category, here referred to as categorical agree-ment p c , is calculated with the help of ν = 0 .
32, the meanuncertainty of β from all proxy records considered. Wethen assign the three categories “stable” ( β < − ν ), “un-stable” (1+ ν ≤ β ) and “intermediate” (1 − ν ≤ β < ν )to the spectral exponent β . The “intermediate” regimeprevents incorrect assignment. To verify the reliability ofcategorical agreement, we calculate the Kappa statistics κ = ( p c − p e ) / (1 − p e ) (2)with expected percent agreement p e by category [81].The latter can be obtained from p e = N P c =1 n c,m n c,p where c is the category, N the number of locations and n the number of times that models ( m ) and proxy records( p ) have predicted category c . The κ -coefficient quan-tifies the reliability from no agreement beyond chance( κ = 0) to full agreement ( κ = 1). Negative κ indicatesagreement that is beyond change, for example, due tosystematic biases. D. Spectral gain
We investigate how climatic drivers influence the globalmean temperature at period τ by calculating the spectralgain G ( τ ) = S T ( τ ) S F ( τ ) . (3)Here, S T ( τ ) is the PSD of the global mean temperatureand S F ( τ ) the PSD of radiative forcing (see also Ap-pendix B). The gain requires the assumption that theglobal mean temperature can be well approximated as alinear function of the forcing [26, 82, 83] and that differ-ent types of radiative forcing add linearly [84–87]. To thisend, we focus on timescales between years and centuries β =2 β =1global mean(a) -1 -2 -3 -6 -3 period τ (yr) PS D S ( τ ) ( K y r) β =2 β =1(b)local mean -1 -2 -3 model simulations: HadCM3 LM1CESM 1 past2kCESM-LME 1MPI-M LMECH5/MPIOM-p6kTraCE-21k-ORBTraCE-21kIPSL-p6k observation-based:
ERA5HadCRUT4PAGES2k period τ (yr)period τ (yr) FIG. 2. (a) Power spectral densities (PSD) of the global mean temperature from model simulations and observation-baseddata on periods from hours to one thousand years for the Holocene. (b) Same as (a) but for the mean PSD of local temperature.The dashed lines with slope β and arbitrary y-intercept in the log-log graph indicate the scaling behavior for visual comparison.The ensemble means (black solid lines) were formed using equal weights across the model group M (Table S1 [44]). where additivity is a valid assumption and nonlineari-ties in the global mean temperature are sufficiently small[39, 40]. The main practical problem that confronts usis that the gain might be subject to a sampling bias dueto our data sets choice. Therefore, we perform a MonteCarlo simulation of the PSD of radiative forcing and theglobal mean temperature, as well as the spectral gain asdescribed in Appendix D. IV. RESULTS AND DISCUSSIONA. Global mean and mean of local spectra
In order to study the stability of global mean temper-ature, we present its power spectral density in Figure 2(a). It shows the characteristic background continuum,spectral peaks, and higher harmonic artifacts associatedwith the diurnal and annual cycle. Overall, the PSDstend to agree between the data sets, albeit with some dif-ferences on the interannual scale and when compared tothe Trace21k ORB run. The Trace21k-ORB run is solelyforced by orbital changes and therefore shows less vari-ability than the ensemble mean. The broad spectral peakon interannual periods reveals an artificially amplifiedENSO in the shared MPI-M LM and ECHAM5/MPI-OM ocean component [88]. For a better visibility, PIcontrol runs are separately shown in the supplementaryFigure S6 [44]. Overall, the PSD largely agrees amongdifferent data sets, especially towards shorter timescales.We find a power-law scaling of β ≈ β ≈ et al. [25], we find no evidence for a scale break aroundthe centennial scale. One limitation of previous workthat found scale breaks is that the spectra were esti- mated across non-stationary shifts in climate, such asthe deglaciation [33], and with a change in proxies andarchives [23], which can lead to spectral artifacts.We present the area-weighted mean spectra of the lo-cal (grid box) temperature in Figure 2 (b). Compared tothe global mean in (a), the power increases and the spec-tral slope decreases, in line with [89]. The spectra agreeon periods below ten years, except for the ENSO arti-fact mentioned earlier. Moreover, we find an unphysical,narrow peak at 13 years, associated with an unrealisticvariability in the northern North Atlantic of the TraCE-21k run, similar to [90]. Remarkably, the decadal-to-centennial variability of the reconstructed temperature isincreased by one to two orders of magnitude compared tothe simulations. The spectral exponents indicate temper-ature stability ( β <
1) from models, whereas it remainsunclear ( β ≈
1) from paleoclimate data.This finding verifies that models show less regionaltemperature variability and that the mismatch increasestowards longer timescales. The results are robust to sam-pling from the PAGES2k database and the influence ofanthropogenic climate change (Figure S10 [44]). Oneshortcoming of forming the area-weighted mean PSD isthat the uncertainty quantification requires the assump-tion of independent spatial degrees of freedom of the tem-perature field. Due to the presence of spatial correlations,an estimate of the effective spatial degrees of freedom andtheir dependence on the underlying timescale would beneeded to resolve this limitation [91].
B. Spatial patterns of stability
To further investigate the mismatch on local scalingproperties, we compare the spatial dependence of tem-perature stability from simulations and paleoclimate data
HadCM3 LM1(a) CESM 1 past2k(b) CESM-LME 1(c)CESM-LME 1 (cntl)(d) MPI-M LM(e) MPI-M LM (cntl)(f)ECH5/MPIOM-p6k(g) TraCE-21k-ORB(h) TraCE-21k(i)IPSL-p6k(j) archive documents lake/marine sediment tree β -1 -0.5 0 0.5 1 1.5 2 2.5 3 FIG. 3. Local temperature persistence on timescales from τ = 10 to τ = 200 years across multiple climate simulationsand selected proxy records from the PAGES2k database. Colors from blue to red indicate the scaling behavior ranging from β = − β = 3. The background of each panel shows the β -values fitted to the PSD of the local grid box temperature fromsimulations. Symbols indicate the scaling of proxy records from different natural archives. in Figure 3. The simulations largely resemble a station-ary (“stable”) variance ( − < β < β > et al. [89] thatthere is no latitudinal dependence of β , in contrast toprevious results [23]. Inspecting the simulations’ β -values(background of Figure 3), we find a small land-sea con- trast. Strongest scaling occurs in the Southern Oceansin line with previous findings [89]. Ocean-sea ice interac-tions with characteristic timescales of the order of cen-turies and a generally increased internal variability overthe oceans might explain these results.We find generally lower values for the slope β in theENSO and Indo-Pacific region. This could be attributedto the fact that (quasi-)oscillatory signals, such as activemodes of internal variability, are reflected in the PSD asbroad peaks and hence cannot be described by a scalinglaw. On the other hand, this finding is stronger in PI con-trol runs compared to fully-forced runs (Figure 3 (c-f)).Thus, residual effects of the recent global warming trendmight play an additional role [92]. A systematic bias be-comes clear from the spatially almost uniform β -valuesof Trace21k-ORB (Figure 3 (h)). In line with Figure 2,we explain this by the lack of forcing mechanisms on in-terannual to multi-decadal timescales.Marine and lake sediments, as well as the archived doc-uments, follow the general trend of increased β -valuescompared to simulations. Tree ring records agree wellwith most simulations in North America and Siberia, butnot necessarily at the coast of Australia and NorthernEurope. Discrepancies such as those in Southern SouthAmerica could reflect the proxies’ strength in represent-ing local conditions, for example, topography. However,noise sources in the climate signal recording and preser-vation, such as bioturbation, can influence proxy records.Further separating the signal content from noise sourcesin paleoclimate reconstructions can help refine our find-ings [93, 94]. C. Statistical agreement of temperature stability
We further investigate the question of temperature sta-bility by a statistical analysis of β -values from simula-tions and paleoclimate reconstructions. It is based onthe detailed uncertainty quantification outlined in sec-tion III B. Our results show that reconstructions andsimulations agree in less than 30% of locations withinthe scope of uncertainties (Figure 4). To single outthe stability of temperature signals, we study the agree-ment by category. We find approximately 25% of agree-ment within the categories β < − ν (“stable”) and β > ν (“unstable”). Although widely accepted [95],categorical and percentage agreement suffer from the lim-itation to ignore any agreement by chance. Therefore,we investigate the κ -statistics (orange bar in Figure 4)and verify that there is no agreement beyond chance( κ = 0) for almost all models. Only MPI-M LM andHadCM3 LM1 show any, if poor agreement ( κ ≈ . κ <
0) due to its systematic bias.The disagreement could be attributed to both paleocli-mate data and simulations. A systematic bias could arise,for example, through the recent, non-stationary globalwarming trend. Therefore, we repeat our analysis withall time series cut at 1850. In particular, anthropogenicwarming slightly increases long-term temperature vari-ability and thus scaling behavior, but not in a significantway (Figures S6, S9, and S10 [44]). Considering paleocli-mate data, we do not expect other systematic biases sincewe base our results on multiple archives and proxies, andno systematic spatial pattern is discernible (Figure S11[44]). Our uncertainty quantification carefully accountsfor potential errors due to irregular sampling and inter-polation by simulating their influence using surrogates(Figure S8 [44]).The models’ resolutions are another element of un-certainty that impacts variability over a wide range oftimescales [96–98]. We here facilitate inter-model com- parison by using state-of-the-art GCMs with compara-ble spatial and temporal resolutions, but computationalcosts precluded higher resolutions. The latter mightbe necessary to improve the representation of decadalvariability and response to external forcing. In particu-lar, the unstable temperature behavior from paleoclimatedata could indicate that nonlinear processes from an in-teractive carbon cycle and dynamical ice-sheets mightnot be sufficiently represented in models.
IPSL-p6kTraCE-21kTraCE-21k-ORBECH5/MPIOM-p6kMPI-M LM (cntl)MPI-M LMCESM-LME 1 (cntl)CESM-LME 1CESM 1 past2kHadCM3 LM1-0.25 0.00 0.25 0.50 0.75 1.00value p / 100%p c / 100% κ FIG. 4. Percentage agreement p , categorical agreement p c and inter-rater reliability κ of local temperature stability fromsimulations and paleoclimate data. The measures were calcu-lated from a set of bilinearly interpolated simulation recordsand the proxy record at 23 different locations. Missing orangebars indicate no agreement beyond chance and, therefore, zerointer-rater reliability ( κ = 0). D. The forced temperature response
Climatic drivers are not constant in time and thus af-fect the surface air temperature on multiple timescales.To investigate the forced temperature response, wepresent spectra for the main climatic drivers in Figure 5.The PSD of orbital forcing consists of the diurnal and an-nual cycle as well as a background continuum on longertimescales. Unphysical higher harmonics on monthlytimescales were omitted. We calculate the mean volcanic,solar, and CO spectra using an equally-weighted aver-age of spectra from multiple data sets (Figure S5 [44]).The CO forcing follows the orbital forcing. The PSD ofsolar forcing again contains more power and has a pro-nounced peak around the eleven-year solar cycle. Multi-ple theories and paleoclimate reconstructions suggest theincreased variability on centennial to millennial periodsdue to the long-term behavior of solar activity [99]. PS D S ( τ ) period τ (yr) FIG. 5. Power spectral densities from radiative forcings. De-tails on the reconstructions considered here are summarizedin Table S2 and Figure S5 [44].
Volcanic forcing dominates interannual to centennialscales and undergoes a scale break around the period often years. Above decadal scales, it follows a white noisespectrum with constant variance. However, the intermit-tency of volcanic eruptions might have led to spectralartifacts [39]. We verify our results using an analyticalapproach described in Appendix C. Remarkably, the de-rived PSD of an ideal, intermittent time series with Pois-son distributed return times explains our findings. Wefurther demonstrate the scale break by a Monte Carlosimulation of the joint PSD of radiative forcing in Fig-ure 6 (a). This finding raises the question of how thespectrum with a scale break translates into the continu-ous spectrum in Figure 2 (a).We address this question by calculating the spectralgain (3) on periods between years and centuries in Fig-ure 6 (b). Here, observation-based data include Had-CRUT4, ERA5, and PAGES2k again. To account forthe unphysical spectral artifacts explained above, we cal-culate the gain from the model simulation group M and together with group M + (supplementary Table S1[44]). We find that the spectral gain is similar fromobservation-based data and the model simulation groupM , which is the one without artificially amplified ENSO.This suggests that both follow a similar distribution oftimescale-dependent variability, as already indicated byFigure 2(a). Large parts of the gain show constant be-havior, which is most pronounced in M . In a simplifiedway, the gain might be approximated by an ideal linearamplifier/damper of the forcing with comparable inter-nal variability on all timescales. However, we also find adip around decadal scales, which is strongest in the gainfrom measurements. Inspecting Figure 6 (a), this can beexplained by forming the ratio between a spectrum witha scale break at the decadal scale ( β > → β ≈
0) andone with moderate scaling ( β ≈ (a) S F S T -1 -2 -3 PS D S T ( τ ) ( K y r) PS D S F ( τ ) ( W m - y r) joint forcingmodel simulations (M ,M + )model simulations (M )observation-based (b) -1 -2 -3 ga i n G ( τ ) ( K ( W m - )) FIG. 6. Monte Carlo simulation of PSD (a) and spectral gain(b) using temperature and forcing reconstructions as well asmodel simulations. Shaded confidence intervals lie betweenthe 5% and 95 % quantiles. We consider only models fromthe groups M and M + (Table S1 [44]) to exclude spectralartifacts and to represent the historical temperature responsein the best possible way. Notably, M + contains those simu-lations with amplified ENSO [88]. (b) Dashed lines indicatethe mean variance ratio h S T i / h S F i . V. CONCLUSION
In summary, we have investigated the question of tem-perature stability on the timescale of years to centuries.To this end, we have presented power spectral densitiesfor both local and global surface air temperature fromsimulation and observation-based data of the last millen-nia. On this basis, we concluded that locally there is astronger scaling and increased variance in reconstructionsas compared to simulations. Using statistical analysis,we found that local temperature series extracted fromsimulations tend to be stable, whereas the proxy recordshint at unstable behavior. Furthermore, we have largelyextended the spectral analysis of climatic drivers by es-timating the joint PSD from CO , solar, volcanic, andorbital forcing using Monte Carlo simulation. Hereby, wediscovered a scale break at the period of approximatelyten years. Moreover, we have presented the spectral gain,describing the timescale-dependent forced temperatureresponse. We found that it is mostly consistent acrossdata sets and indicates an increasing internal variabilityon timescales of decades to centuries.Our analysis of the spectral gain was limited to globalaverage values and those timescales where linearity canbe reasonably assumed [39, 40, 100]. Nonlinearities areinherent to the climate system, for example, due to thetemperature-albedo feedback. Thus, it will be necessaryto examine their possible effects on multiple spatiotem-poral scales to further extend this work. Studying non-linearities could also shine new light on the mechanismsof scaling in Earth’s climate, which are not yet fully un-derstood and might be linked to nonlinearities as well[6]. Furthermore, we have focused on the current inter-glacial, the Holocene. This is because climate variabilityhas been demonstrated to depend on the mean climatestate [101] and major shifts in climate can lead to spec-tral artifacts. Thus, the conclusions laid out here cannotbe readily applied to other climate states, such as glacialperiods, which is an issue for future studies. Clearly, un-derstanding the dependence of temperature variability onglobal warming demands additional work.Ideally, our findings should be replicated by employ-ing models with increased internal variability on longertimescales and paleoclimate data that provides improvedspatiotemporal resolution. Optimized analysis of noisesources and spectral analysis of (pseudo-)proxy recordscould help to expand our data basis of proxy records withdecadal resolution [56, 102, 103]. Regarding climate mod-els, an improved representation of processes that increaseEarth’s long-term memory, such as an interactive car-bon cycle and dynamical ice-sheets, might strengthen thelong-range dependence and persistence of surface air tem-perature. Future studies could also continue to explorehow internally-generated and externally-forced variabil-ity compares on different spatial scales. Research on theinterrelation between internal and forced changes, as wellas local, regional, and global variability, might prove im-portant and could be conducted using single-forcing ex- periments from ensembles of model simulations.Managing climate risks requires a detailed understand-ing of temperature variability. Locally and on time scalesbetween years and centuries, there is an urgency to ad-dress discrepancies to make further progress in climatemodeling. In this study, we have singled out the key char-acteristics of temperature variability and showed thatobservation-based data and model simulations assess thestability of local temperatures differently. Our resultshave demonstrated the stability and spectral gain as aneasy-to-use yet effective and promising tool for investi-gating variability in Earth’s dynamic climate. ACKNOWLEDGMENTS
This manuscript is based upon data provided by theWorld Climate Research Programme’s Working Groupon Coupled Modelling, which is responsible for CMIPand PMIP. We thank the research groups listed in sup-plementary Tables S1 and S2 for producing and mak-ing available their data from model outputs, measure-ments, paleoclimate, and forcing reconstructions. Thisstudy benefited from discussions within the CVAS work-ing group, a working group of the Past Global Changes(PAGES) project. We thank T. Gasenzer, T. Kunz, andN. Weitzel for discussions and J. B¨uhler, M. Casado,M. Schillinger, and E. Ziegler for helpful comments onthe manuscript. This research has been funded by theHeidelberg Graduate School for Physics, by the PalModproject ( , subproject no. 01LP1926C)and by the Deutsche Forschungsgemeinschaft (DFG, Ger-man Research Foundation) – project no. 395588486.Code to reproduce all figures in this manuscript will bemade available upon publication.
Appendix A: Relation between power spectraldensity and variance
The power spectral density of a weakly-stationary,stochastic process, is given by the Fourier transform ofthe autocorrelation S ( f ) = F{ R ( h ) } with frequency f and lag h = t − t between two points in time [18, 19].For zero lag and zero mean, the integral of the PSD cor-responds to the variance of the signal [77]. Instead offrequency, we use the period τ = 1 /f to express the PSDand spectral gain. Appendix B: Spectral gain for linear systems
In a time-invariant linear system, the output y ( t ) = Z ∞−∞ h ( u ) x ( t − u )d u (B1)is given by the input time series x ( t ) and the impulseresponse function h ( u ) [77]. The Fourier transform0 H ( f ) = F{ h ( u ) } = G ( f )e i φ ( f ) gives the frequency re-sponse function, also called transfer function. G ( f ) and φ ( f ) are the gain and phase, respectively. The inte-gral (B1) corresponds to a product in frequency space F{ y ( t ) } = H ( f ) F{ x ( t ) } . This relates the PSD of theoutput S y ( f ) to the one of the input S x ( f ) via S y ( f ) = | H ( f ) | S x ( f ) = G ( f ) S x ( f ) . (B2) Appendix C: Analytical solution to the PSD ofintermittent volcanic forcing
We investigate the power spectral density of inter-mittent volcanic forcing by approximating the eruptiontime series in a simplified way as a stochastic signal X ( t ) = δ ( t − t i ). This function is zero at all times except t i , when an event of unique amplitude occurs. We de-note T i = t i − t i − the time intervals between two events.We use the fact that the PSD cannot only be calculatefrom the covariance but also from the Laplace transform S ( X, f ) = 2 lim (cid:15) → h|L ( X ( t ) , (cid:15) − π i f ) | i [104]. Based onthis approach, the power spectral density S ( f ) = µ T − | ρ ( f ) | | − ρ ( f ) | , f > ρ ( f ) = F{ ρ ( T ) } and the in-verse mean interval between two events µ T = h T i − [104, 105]. An exponentially decaying probability dis-tribution ρ ( T ) = µ T exp( − µ T T ) for volcanic forcing issuggested [106], and we have checked this for the datasets considered. The Fourier transform reads ρ ( f ) =( µ T +2 π i f ) − such that 1 −| ρ ( f ) | = | − ρ ( f ) | . As a con-sequence, the PSD (C1) takes a constant value. We canobserve this white noise behavior in Figure 5 and Figure 6(a) on timescales longer than a few years, which is on theorder of characteristic return times for eruptions. Belowthese timescales, the variability considerably drops. Thisanalytical result provides an independent verification ofthe PSD for volcanic forcing and its scale break. Appendix D: Monte Carlo sampling of the spectralgain
We simulate the spectral gain (3), as well as the PSD ofglobal mean temperature and the joint PSD of radiativeforcing using a Monte Carlo approach with N = 1000realizations to account for sampling biases. The PSDof global mean temperature is sampled for three groups:the observation-based data, the model simulations fromgroup M , and from M together with M + (Table S1[44]). Here, only models from the groups M and M + areconsidered to exclude spectral artifacts and to representthe historical temperature response in the best possibleway.We sample the simulation-based PSD from the aver-age PSD of the simulations using uniformly distributedrandom weights. To obtain the observation-based PSD,we use the global mean temperature from HadCRUT4,ERA5, and a 7000-member reconstruction ensemble pro-vided by PAGES2k [43]. This ensemble allows us to sam-ple the PSD by randomly selecting one ensemble memberand form the mean of its spectrum with that of the ERA5and HadCRUT4 temperature. The joint PSD of radia-tive forcing is calculated from all forcing reconstructionsconsidered in this work except the Fr¨ohlich et al. so-lar forcing, which has too low temporal resolution aboveinterannual scales (Table S3 and Figure S5 [44]). We as-sume the PSD of CO and orbital forcing as fixed sinceits spectral power is comparatively low on multi-decadalscales. We sample the PSD of solar forcing by usinguniformly distributed weights when forming the averagePSD of all solar reconstructions. Similarly, the PSD ofvolcanic forcing is obtained. In addition, we randomlyvary the conversion factor between ( − − and ( − − Wm − /AOD [75]. The joint PSD of radiative forcing iscalculated by linear summation of the PSD from CO ,orbital, solar, and volcanic forcing.Using this sampling scheme, our Monte Carlo producestwo outcomes: First, we compute the PSD of global meantemperature and the joint PSD of radiative forcing bysimulating an ensemble of N realizations for both forc-ing and response. Second, we sample the spectral gain di-rectly from the quotient (3) in each of the N realizations.In both cases, the average of the generated N -memberensemble and its 5% and 95% quantiles constitute theresult of our Monte Carlo simulation. [1] J. Bjerknes, A possible response of the atmosphericHadley circulation to equatorial anomalies of ocean tem-perature, Tellus , 820 (1966).[2] S. Arrhenius, XXXI. On the influence of carbonic acid in the air upon the temperature of the ground, TheLondon, Edinburgh, and Dublin Philosophical Maga-zine and Journal of Science , 237 (1896).[3] J. B. J. Fourier, Remarques g´en´erales sur les temp´eratures du globe terrestre et des espacesplan´etaires, Annales de Chimie et de Physique , 136(1824).[4] R. W. Katz and B. G. Brown, Extreme events in achanging climate: Variability is more important thanaverages, Climatic Change , 289 (1992).[5] M. Ghil and V. Lucarini, The physics of climate vari-ability and climate change, Reviews of Modern Physics , 10.1103/RevModPhys.92.035002 (2020).[6] C. L. E. Franzke et al., The Structure of Cli-mate Variability Across Scales, Reviews of Geophysics10.1029/2019RG000657 (2020).[7] J. E. Tierney et al., Past climates inform our future,Science , 10.1126/science.aay3701 (2020).[8] T. Laepple and P. Huybers, Global and regional vari-ability in marine surface temperatures, Geophysical Re-search Letters 10.1002/2014GL059345 (2014).[9] T. Laepple and P. Huybers, Ocean surface temperaturevariability: Large model-data differences at decadal andlonger periods, Proceedings of the National Academy ofSciences of the United States of America , 16682(2014).[10] L. A. Parsons, G. R. Loope, J. T. Overpeck, T. R.Ault, R. Stouffer, and J. E. Cole, Temperature and pre-cipitation variance in CMIP5 simulations and paleocli-mate records of the last millennium, Journal of Climate10.1175/JCLI-D-16-0863.1 (2017).[11] F. C. Ljungqvist, Q. Zhang, G. Brattstr¨om, P. J. Krusic,A. Seim, Q. Li, Q. Zhang, and A. Moberg, Centennial-Scale Temperature Change in Last Millennium Simula-tions and Proxy-Based Reconstructions, Journal of Cli-mate , 2441 (2019).[12] J. C. B¨uhler, C. Roesch, M. Kirschner, L. Sime, M. D.Holloway, and K. Rehfeld, Comparison of the oxygenisotope signatures in speleothem records and iHadCM3model simulations for the last millennium, Climate ofthe Past Discussions , 1 (2020).[13] W. B. Anderson, R. Seager, W. Baethgen, M. Cane, andL. You, Synchronous crop failures and climate-forcedproduction variability, Science Advances , 10.1126/sci-adv.aaw1976 (2019).[14] M. Crucifix, A. D. Vernal, and C. Franzke, CentennialTo Millennial Climate Variability, PAGES (2017).[15] M. Rypdal, H. B. Fredriksen, E. Myrvoll-Nilsen,K. Rypdal, and S. H. Sørbye, Emergent scale invarianceand climate sensitivity, Climate , 1 (2018).[16] T. P. Barnett et al., Detection and attribution ofrecent climate change: A status report, Bulletin ofthe American Meteorological Society , 10.1175/1520-0477(1999)080¡2631:DAAORC¿2.0.CO;2 (1999).[17] N. Bindoff et al., Detection and attribution of climatechange: From global to regional, in Climate Change2013: The Physical Science Basis. Contribution ofWorking Group I to the Fifth Assessment Report of theIntergovernmental Panel on Climate Change , edited byI. P. o. C. Change (Cambridge University Press, 2013).[18] N. Wiener, Generalized harmonic analysis, Acta Math-ematica , 117 (1930).[19] A. Khintchine, Korrelationstheorie der station¨arenstochastischen Prozesse, Mathematische Annalen ,604 (1934).[20] C. Wunsch, The spectral description of climate changeincluding the 100 ky energy, Climate Dynamics , 353(2003). [21] K. Fraedrich, U. Luksch, and R. Blender, A 1/f-modelfor long time memory of the ocean surface temperature,Phys. Rev. E , 10.1103/PhysRevE.70.037301 (2003).[22] C. Franzke, Nonlinear trends, long-range dependence,and climate noise properties of surface temperature,Journal of Climate 10.1175/JCLI-D-11-00293.1 (2012).[23] P. Huybers and W. Curry, Links between annual, Mi-lankovitch and continuum temperature variability, Na-ture , 329 (2006).[24] S. Lovejoy, A voyage through scales, a missingquadrillion and why the climate is not what you expect,Climate Dynamics , 3187 (2015).[25] T. Nilsen, K. Rypdal, and H. B. Fredriksen, Arethere multiple scaling regimes in Holocene temperaturerecords?, Earth System Dynamics , 419 (2016).[26] H. B. Fredriksen and M. Rypdal, Long-range persis-tence in global surface temperatures explained by lin-ear multibox energy balance models, Journal of Climate10.1175/JCLI-D-16-0877.1 (2017).[27] B. B. Mandelbrot and J. W. Van Ness, FractionalBrownian Motions, Fractional Noises and Applications,SIAM Review , 10.1137/1010093 (1968).[28] J. P. Peixto and A. H. Oort, Physics of climate, Reviewsof Modern Physics , 365 (1984).[29] E. J. Rohling et al., Making sense of palaeoclimate sen-sitivity, Nature , 683 (2012).[30] E. J. Rohling, G. Marino, G. L. Foster, P. A. Goodwin,A. S. von der Heydt, and P. K¨ohler, Comparing ClimateSensitivity, Past and Present, Annual Review of MarineScience , 261 (2018).[31] F. Zhang, Y. Qiang Sun, L. Magnusson, R. Buizza, S. J.Lin, J. H. Chen, and K. Emanuel, What is the pre-dictability limit of midlatitude weather?, Journal of theAtmospheric Sciences 10.1175/JAS-D-18-0269.1 (2019).[32] J. D. Pelletier, Natural variability of atmospheric tem-peratures and geomagnetic intensity over a wide rangeof time scales, Proceedings of the National Academyof Sciences of the United States of America , 2546(2002).[33] F. Zhu, J. Emile-Geay, N. P. McKay, G. J. Hakim,D. Khider, T. R. Ault, E. J. Steig, S. Dee, and J. W.Kirchner, Climate models can correctly simulate thecontinuum of global-average temperature variability,Proceedings of the National Academy of Sciences of theUnited States of America , 8728 (2019).[34] J. Lohmann and P. D. Ditlevsen, Random and ex-ternally controlled occurrences of Dansgaard–Oeschgerevents, Climate of the Past , 609 (2018).[35] K. Rypdal, L. Østvand, and M. Rypdal, Long-rangememory in Earth’s surface temperature on time scalesfrom months to centuries, Journal of Geophysical Re-search: Atmospheres , 7046 (2013).[36] K. Marvel, G. A. Schmidt, R. L. Miller, and L. S.Nazarenko, Implications for climate sensitivity from theresponse to individual forcings, Nature Climate Change , 386 (2016).[37] A. P. Schurer, G. C. Hegerl, M. E. Mann, S. F. B. Tett,and S. J. Phipps, Separating Forced from Chaotic Cli-mate Variability over the Past Millennium, Journal ofClimate , 6954 (2013).[38] S. Lovejoy and D. Schertzer, Stochastic and scal-ing climate sensitivities: Solar, volcanic and or-bital forcings, Geophysical Research Letters ,https://doi.org/10.1029/2012GL051871 (2012). [39] S. Lovejoy and C. Varotsos, Scaling regimes and lin-ear/nonlinear responses of last millennium climate tovolcanic and solar forcings, Earth System Dynamics10.5194/esd-7-133-2016 (2016).[40] K. Rypdal and M. Rypdal, Comment on “Scalingregimes and linear/nonlinear responses of last millen-nium climate to volcanic and solar forcing” by S. Love-joy and C. Varotsos (2016), Earth System Dynamics ,597 (2016).[41] J. Mitchell, An overview of climatic variability and itscausal mechanisms, Quaternary Research , 481 (1976).[42] A. Henderson-Sellers and K. McGuffie, A climate mod-elling primer (Wiley, Chichester, UK, 1987).[43] PAGES 2k Consortium., Consistent multidecadal vari-ability in global temperature reconstructions and sim-ulations over the Common Era, Nature Geoscience ,643 (2019).[44] See supplementary material for futher information ondatasets and benchmarks.[45] B. L. Otto-Bliesner, E. C. Brady, J. Fasullo, A. Jahn,L. Landrum, S. Stevenson, N. Rosenbloom, A. Mai, andG. Strand, Climate variability and change since 850 cean ensemble approach with the community earth systemmodel, Bulletin of the American Meteorological Society10.1175/BAMS-D-14-00233.1 (2016).[46] J. H. Jungclaus et al., Climate and carbon-cycle vari-ability over the last millennium, Climate of the Past10.5194/cp-6-723-2010 (2010).[47] Y. Zhong, A. Jahn, G. H. Miller, and A. Geirsdottir,Asymmetric Cooling of the Atlantic and Pacific ArcticDuring the Past Two Millennia: A Dual Observation-Modeling Study, Geophysical Research Letters ,12,412 (2018).[48] P. Braconnot, D. Zhu, O. Marti, and J. Servonnat,Strengths and challenges for transient Mid- to LateHolocene simulations with dynamical vegetation, Cli-mate of the Past , 997 (2019).[49] N. Fischer and J. H. Jungclaus, Evolution of the sea-sonal temperature cycle in a transient Holocene simu-lation: Orbital forcing and sea-ice, Climate of the Past10.5194/cp-7-1139-2011 (2011).[50] Z. Liu, Transient simulation of last deglaciation witha new mechanism for Bølling–Allerød warming, Science , 10.1126/science.1171041 (2009).[51] C. P. Morice, J. J. Kennedy, N. A. Rayner, andP. D. Jones, Quantifying uncertainties in global andregional temperature change using an ensemble ofobservational estimates: The HadCRUT4 data set,Journal of Geophysical Research Atmospheres ,10.1029/2011JD017187 (2012).[52] H. Hersbach, B. Bell, P. Berrisford, and S. Hirahara,The ERA5 global reanalysis, Quarterly Journal of theRoyal Meteorological Society , 1999 (2020).[53] PAGES2k Consortium., A global multiproxy databasefor temperature reconstructions of the Common Era,Scientific Data 10.1038/sdata.2017.88 (2017).[54] J. W. Kantelhardt, Fractal and Multifractal Time Se-ries, in Mathematics of Complexity and Dynamical Sys-tems , edited by R. A. Meyers (Springer New York, NewYork, NY, 2011) pp. 463–487.[55] T. Laepple, T. M¨unch, M. Casado, M. Hoerhold,A. Landais, and S. Kipfstuhl, On the similarity and ap-parent cycles of isotopic variations in East Antarcticsnow pits, Cryosphere 10.5194/tc-12-169-2018 (2018). [56] M. Casado, T. M¨unch, and T. Laepple, Climatic infor-mation archived in ice cores: impact of intermittencyand diffusion on the recorded isotopic signal in Antarc-tica, Climate of the Past , 1581 (2020).[57] T. J. Crowley, Causes of climate change over thepast 1000 years, Science 10.1126/science.289.5477.270(2000).[58] G. C. Hegerl, T. J. Crowley, M. Allen, W. T. Hyde,H. N. Pollack, J. Smerdon, and E. Zorita, Detection ofHuman Influence on a New, Validated 1500-Year Tem-perature Reconstruction, Journal of Climate , 650(2007).[59] P. Braconnot, S. P. Harrison, M. Kageyama, P. J.Bartlein, V. Masson-Delmotte, A. Abe-Ouchi, B. Otto-Bliesner, and Y. Zhao, Evaluation of climate models us-ing palaeoclimatic data, Nature Climate Change , 417(2012).[60] G. A. Schmidt et al., Climate forcing reconstructions foruse in PMIP simulations of the Last Millennium (v1.1),Geoscientific Model Development , 185 (2012).[61] G. Delaygue and E. Bard, An Antarctic view ofBeryllium-10 and solar activity for the past millennium,Climate Dynamics , 2201 (2011).[62] F. Steinhilber, J. Beer, and C. Fr¨ohlich, Total solar irra-diance during the Holocene, Geophysical Research Let-ters , 10.1029/2009GL040142 (2009).[63] Y. Wang, J. L. Lean, and N. R. Sheeley, Jr., Modelingthe Sun’s Magnetic Field and Irradiance since 1713, TheAstrophysical Journal 10.1086/429689 (2005).[64] R. Muscheler, F. Joos, J. Beer, S. A. M¨uller, M. Von-moos, and I. Snowball, Solar activity during the last1000yr inferred from radionuclide records, QuaternaryScience Reviews , 82 (2007).[65] L. E. Vieira and S. K. Solanki, Evolution of the solarmagnetic flux on time scales of years to millenia, Astron-omy and Astrophysics 10.1051/0004-6361/200913276(2010).[66] C. Gao, A. Robock, and C. Ammann, Volcanic forcing ofclimate over the past 1500 years: An improved ice core-based index for climate models, Journal of GeophysicalResearch: Atmospheres , 10.1029/2008JD010239(2008).[67] A. L. Berger, Long-term variations of dailyinsolation and Quaternary climatic changes.,Journal of Atmospheric Sciences 10.1175/1520-0469(1978)035¡2362:ltvodi¿2.0.co;2 (1978).[68] M. Crucifix, Palinsol - Package (R) (2016).[69] M. Toohey and M. Sigl, Reconstructed volcanic strato-spheric sulfur injections and aerosol optical depth, 500BCE to 1900 CE, version 2, World Data Center for Cli-mate (WDCC) at DKRZ 10.1594/WDCC/eVolv2k v2(2017).[70] C. Fr¨ohlich, Solar irradiance variability since 1978: Re-vision of the PMOD composite during solar cycle 21,Space Science Reviews , 53 (2006).[71] C. D. Keeling, R. B. Bacastow, and A. E. Bain-bridge, Atmospheric carbon dioxide variations atMauna Loa Observatory, Hawaii, TELLUS 10.3402/tel-lusa.v28i6.11322 (1976).[72] J. Pongratz, C. Reick, T. Raddatz, and M. Claussen,A reconstruction of global agricultural areas and landcover for the last millennium, Global BiogeochemicalCycles , 10.1029/2007GB003153 (2008).[73] G. Myhre, E. J. Highwood, K. P. Shine, and F. Stordal, New estimates of radiative forcing due to well mixedgreenhouse gases, Geophysical Research Letters ,2715 (1998).[74] G. A. Schmidt et al., Climate forcing reconstructions foruse in PMIP simulations of the last millennium (v1.0),Geoscientific Model Development , 33 (2011).[75] G. Myhre et al., Anthropogenic and Natural RadiativeForcing Supplementary Material., in Climate Change2013: The Physical Science Basis. Contribution ofWorking Group I to the Fifth Assessment Report ofthe Intergovernmental Panel on Climate Change , editedby Intergovernmental Panel on Climate Change (Cam-bridge University Press, Cambridge, 2014) pp. 1–30.[76] T. Laepple and P. Huybers, Reconciling discrepanciesbetween Uk37 and Mg/Ca reconstructions of Holocenemarine temperature variability, Earth and PlanetaryScience Letters , 418 (2013).[77] C. Chatfield and H. Xing,
The Analysis of Time Series:An Introduction with R , seventh ed., Chapman & Hall/ CRC Texts in Statistical Science (CRC Press, 2019).[78] D. B. Percival and A. T. Walden,
Spectral Analysisfor Physical Applications: Multitaper and ConventionalUnivariate Techniques (Cambridge University Press,Cambridge, UK, 1993).[79] P. Yiou, E. Baert, and M. F. Loutre, Spectral analysisof climate data, Surveys in Geophysics , 619 (1996).[80] J. W. Kirchner, Aliasing in 1 /f a noise spectra: Origins,consequences, and remedies, Phys. Rev. E , 66110(2005).[81] J. L. Fleiss and J. Cohen, The Equivalence of WeightedKappa and the Intraclass Correlation Coefficient asMeasures of Reliability, Educational and PsychologicalMeasurement , 613 (1973).[82] O. Geoffroy, D. Saint-martin, D. J. Olivi´e, A. Voldoire,G. Bellon, and S. Tyt´eca, Transient climate response ina two-layer energy-balance model. Part I: Analytical so-lution and parameter calibration using CMIP5 AOGCMexperiments, Journal of Climate , 1841 (2013).[83] D. G. MacMynowski, H.-J. Shin, and K. Caldeira, Thefrequency response of temperature and precipitationin a climate model, Geophysical Research Letters ,https://doi.org/10.1029/2011GL048623 (2011).[84] A. Kirkev˚ag, T. Iversen, J. E. Kristj´ansson, Ø. Seland,and J. B. Debernard, On the additivity of climate re-sponse to anthropogenic aerosols and CO2, and the en-hancement of future global warming by carbonaceousaerosols, Tellus, Series A: Dynamic Meteorology andOceanography
60 A , 10.1111/j.1600-0870.2008.00308.x(2008).[85] H. Shiogama, D. A. Stone, T. Nagashima, T. Nozawa,and S. Emori, On the linear additivity of climateforcing-response relationships at global and continen-tal scales, International Journal of Climatology ,10.1002/joc.3607 (2013).[86] G. A. Meehl, W. M. Washington, C. M. Ammann, J. M.Arblaster, T. M. Wigley, and C. Tebaldi, Combina-tions of natural and anthropogenic forcings in twentieth-century climate, Journal of Climate , 10.1175/1520-0442(2004)017¡3721:CONAAF¿2.0.CO;2 (2004).[87] V. Ramaswamy and C. T. Chen, Linear additivityof climate response for combined albedo and green-house perturbations, Geophysical Research Letters ,10.1029/97GL00248 (1997). [88] J. Jungclaus, (private communications) (2020).[89] H. B. Fredriksen and K. Rypdal, Spectral characteris-tics of instrumental and climate model surface temper-atures, Journal of Climate , 1253 (2016).[90] G. Danabasoglu, On Multidecadal Variability of the At-lantic Meridional Overturning Circulation in the Com-munity Climate System Model Version 3, Journal of Cli-mate , 5524 (2008).[91] T. Kunz and T. Laepple, Time-scale dependent esti-mation of spatial degrees of freedom, in EGU GeneralAssembly 2018 (2018).[92] S.-W. Yeh, J.-S. Kug, B. Dewitte, M.-H. Kwon, B. P.Kirtman, and F.-F. Jin, El Ni˜no in a changing climate,Nature , 511 (2009).[93] M. Reschke, K. Rehfeld, and T. Laepple, Empirical es-timate of the signal content of Holocene temperatureproxy records, Climate of the Past , 521 (2019).[94] L. J. L¨ucke, G. C. Hegerl, A. P. Schurer, and R. Wilson,Effects of memory biases on variability of temperaturereconstructions, Journal of Climate , 10.1175/JCLI-D-19-0184.1 (2019).[95] J. L. Fleiss, B. Levin, and M. C. Paik, Statistical Meth-ods for Rates and Proportions , 3rd ed., Wiley Series inProbability and Statistics (Wiley, 2013).[96] B. P. Kirtman et al., Impact of ocean model resolutionon CCSM climate simulations, Climate Dynamics ,1303 (2012).[97] J. M. Klavans, A. Poppick, S. Sun, and E. J. Moyer, Theinfluence of model resolution on temperature variability,Climate Dynamics , 3035 (2017).[98] D. L. Hodson and R. T. Sutton, The impact of reso-lution on the adjustment and decadal variability of theAtlantic meridional overturning circulation in a coupledclimate model, Climate Dynamics , 10.1007/s00382-012-1309-0 (2012).[99] L. J. Gray et al., Solar influences on climate, Reviewsof Geophysics , 10.1029/2009RG000282 (2010).[100] L. Fern´andez-Donado et al., Large-scale temperature re-sponse to external forcing in simulations and reconstruc-tions of the last millennium (2013).[101] K. Rehfeld, T. M¨unch, S. L. Ho, and T. Laepple, Globalpatterns of declining temperature variability from theLast Glacial Maximum to the Holocene, Nature ,356 (2018).[102] T. Kunz, A. M. Dolman, and T. Laepple, A spectralapproach to estimating the timescale-dependent uncer-tainty of paleoclimate records – Part 1: Theoretical con-cept, Climate of the Past , 1469 (2020).[103] A. Dolman, T. Kunz, J. Groeneveld, and T. Laepple,Estimating the timescale-dependent uncertainty of pa-leoclimate records – a spectral approach. Part II: Appli-cation and interpretation, Climate of the Past Discus-sions 10.5194/cp-2019-153 (2020).[104] R. L. Stratonovich, Topics in the Theory of RandomNoise (Gordon and Breach, New York, New York, USA,1967).[105] B. Lindner, Superposition of many independent spiketrains is generally not a Poisson process, Physical Re-view E - Statistical, Nonlinear, and Soft Matter Physics , 1 (2006).[106] P. Papale, Global time-size distribution of volcanic erup-tions on Earth, Scientific Reports , 6838 (2018). upplementary Material for Probing the timescale-dependent stability of local and globalsurface air temperature from climate simulations andreconstructions of the last millennia
Beatrice Ellerhoff and Kira Rehfeld
LIST OF TABLES
S1 Key specifications of model simulation and observation-based data . . . . . . . . . 2S2 Key specification of climatic drivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3S3 Key specification of paleoclimate data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
LIST OF FIGURES
S1 Global mean temperature over the Common Era . . . . . . . . . . . . . . . . . . . . . . . . 7S2 Climatic drivers over the last millennium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7S3 Temperature signals from proxy records over the last millennia . . . . . . . . . . . . 8S4 PSD of proxy records . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9S5 Mean PSDs for radiative forcing reconstructions . . . . . . . . . . . . . . . . . . . . . . . . . 10S6 PSD from last millennium runs and their PI control runs . . . . . . . . . . . . . . . . . 10S7 Mean PSD from high- and low-frequency spectrum . . . . . . . . . . . . . . . . . . . . . . 11S8 Uncertainty quantification for irregularly sampled temperature proxies . . . . . 11S9 The effect of sampling on the mean of local PSDs from proxy records . . . . . . . 12S10 The effect of global warming on the spectral exponent β . . . . . . . . . . . . . . . . . . 12S11 Scaling of local temperature from simulation and proxies . . . . . . . . . . . . . . . . . 13S12 Goodness of fit for local scaling coefficients from model simulations . . . . . . . . 141 a r X i v : . [ phy s i c s . a o - ph ] J a n ABLE S1. Key specifications of model simulations and observation-based data used to estimatetemperature variability. We give the main references (Doc.) of the simulation runs, observation-based data sets, and their forcing (Forc.). The groups M + and M were assigned to better dis-tinguish spectral properties when performing a Monte Carlo simulation of the spectral gain in themain manuscript. The temporal resolution ∆ t is given in months (m), hours (hr), years (yrs), anddays. The spatial resolutions of atmosphere and ocean are denoted by the subscripts ◦ deg at. and ◦ deg oc. , respectively. Name Doc. Model ∆ t ◦ deg at. ◦ deg oc. Forc. Time (CE)
Model simulations
CESM 1 past2k a [1] CESM1 1 m 2 1 [2] 1-2005CESM-LME 1 a [3] CESM1 1 m, 6 hr 2 1 [4] 850-2006CESM-LME 1 cntl [3] CESM1 1 m 2 1 [4] 850-2006MPI-M LM b [5] MPI-ESM 1 m, 6 hr 3.75 GR30 c [5] 800-2005MPI-M LM cntl [5] MPI-ESM 1 m 3.75 GR30 c [5] 800-2005HadCM3 LM1 a [6] iHadCM3 1 m 2.5 x 3.75 1.25 [7] 800-1950IPSL-p6k a [8] IPSL-CM5A 2 m 2.5 x 1.27 2 [8] 4000 BCE - 2000TraCE-21k [9] CCSM3 2 m ≈ .
75 3.6 x v d [10, 11] 4,000 BCE - 1990TraCE-21k-ORB [9] CCSM3 10 yrs ≈ .
75 3.6 x v d [10, 11] 4,000 BCE - 1990ECH5/MPIOM-p6k b [12] ECHAM5/MPI-OM 1 m 3.75 2 [10] 4000 BCE - 2000 Observation-based
ERA5 [13] 1 m, 6 hr 2 2 1979-2019HadCRUT4 [14] 7 days 5 5 1850-2019PAGES2k [15] 1 yr [15] 0-2000 a assigned to group M assigned to M +c curvilinear grid with nominal resolution of 3.0 ◦ d the latitudinal resolution is variable (v), with finer resolution near the equator ( ≈ . ◦ ) ABLE S2. Key specifications of climatic drivers used to estimate the power spectral density ofradiative forcing as well as the gain function of the forced temperature response.Name ∆ t Time (CE)
Volcanic forcing
Crowley et al. a [17] 10 days 500BCE-1900Gao et al. a [18] 1 yr 850-2000Toohey et al. [19] 850-2000 Total solar irradiance
Delaygue et al. ab [20] 1 yr 850-1850Muscheler et al. ab [21] 1 yr 850-1850Steinhilber et al. a [22] 1 yr 850-1850Vieira et al. a [23, 24] 1 yr 850-1850Wang et al. ab [25] 1 yr 1610-2009Fr¨ohlich et al. [26] 36 hr 1978-2017 CO Schmidt et al. a [16] 1 yr 850-2000Keeling et al. [27] 1 hr 1970-2016 Insolation at 65 ◦ N Berger ac [10] 1 yr, 3 hr 0-2000 a from the PMIP simulations of the Last Millennium [16] b when multiple versions were provided by PMIP, the “with-background”-version was considered here c computed with the R -package ”Palinsol” [28] A B L E S . K e y s p ec i fi c a t i o n o f p r o xy r ec o r d s u s e d t o e s t i m a t e l o c a l t e m p e r a t u r e v a r i a b ili t y . T h e fi r s t s i x c o l u m n s ( “ I D ”” - “ P r o xy ”” ) a r e t a k e n f r o m t h e P A G E S k d a t a b a s e [ ]. T h e I D i s a dd i t i o n a ll y m a r k e d w i t h a s t a r t( ? ) i n c a s e t h e p r o xy w a s u s e d f o r c a l c u l a t i n g t h e ag r ee m e n t i n s c a li n g b e h a v i o r b e t w ee n m o d e l s a ndd a t a . T h ec u t o ff i d e n t i fi e s w h e t h e r t h e p r o xy w a s u s e d f o r a n a l y z i n g t h e p o t e n t i a li m p a c t o f t h e r ece n t g l o b a l w a r m i n g t r e nd ( F i g u r e S nd F i g u r e S ) . I t i nd i c a t e s w h e t h e r t h e n ece ss a r y c r i t e r i a f o r t i m e s e r i e ss e l ec t i o n w e r e f u l fi ll e db e f o r e ( P I ) a nd /o r o v e r t h e f u ll h i s t o r i c a l ( h i s t) p e r i o d . S i m il a r l y , t h e s e l ec t i o n c o l u m n g i v e s c r i t e r i a f o r c a l c u l a t i n g t h e m e a n o f l o c a l s p ec t r aa nd t h e i r c o m p a r i s o n ( F i g u r e S ) . T h e t a b l e s p a n s m u l t i p l e p ag e s . I D N a m e L a t L o n E l e v m a s l A r c h i v e P r o xy c u t o ff s e l ec t i o n A f r i c a - L a k e T a n ga n y i. T i e r n e y . − . l a k e s e d i m e n t T E X P I / h i s t l oo s e A f r i c a - M a l a w i. P o w e r s . − . l a k e s e d i m e n t T E X P I / h i s t l oo s e A f r i c a - P - . T i e r n e y . . m a r i n e s e d i m e n t T E X P I / h i s t s t r o n g/ l oo s e A r c - B r a y a S o . D A nd r e a . − . − l a k e s e d i m e n t a l k e n o n e P I / h i s t l oo s e A r c - C l e gg201061 . − . − l a k e s e d i m e n t m i d g e P I / h i s t l oo s e A r c - G u l f o f A l a s k a . W il e s . − . − t r ee T R W P I / h i s t s t r o n g/ l oo s e A r c - H a ll e t L a k e . M c K a y . . − . − l a k e s e d i m e n t B S i P I / h i s t s t r o n g/ l oo s e A r c - I ce l a nd . B e r g t h o r ss o n . . − . − d o c u m e n t s h i s t o r i c P I / h i s t l oo s e A r c - K o n g r e ss v a t n . D A nd r e a . . l a k e s e d i m e n t a l k e n o n e P I / h i s t l oo s e A r c - L a k e . R o ll a nd . . − . − l a k e s e d i m e n t c h i r o n o m i d P I / h i s t l oo s e A r c - L a k e E . D A nd r e a . − . − l a k e s e d i m e n t a l k e n o n e P I / h i s t l oo s e A r c - L a k e P i e n i - K a . L u o t o . . . l a k e s e d i m e n t c h i r o n o m i d P I / h i s t l oo s e A r c - L u o t o200960 . . l a k e s e d i m e n t m i d g e P I / h i s t s t r o n g/ l oo s e A r c - M a c k e n z i e D e l t a . P o r t e r . . − . − t r ee m i ss i n g h i s t s t r o n g/ l oo s e A r c - M D . J i a n g . . − . − m a r i n e s e d i m e n t d i a t o m P I / h i s t s t r o n g/ l oo s e A r c - S o p e r L a k e B a f . H u g h e n . . − . − l a k e s e d i m e n t v a r v e t h i c k n e ss P I / h i s t s t r o n g/ l oo s e A r c - S t o r e gga S li d . S e j r up . . . m a r i n e s e d i m e n t f o r a m d O P I / h i s t s t r o n g/ l oo s e A r c - T h o m a s . − . − l a k e s e d i m e n t v a r v e t h i c k n e ss P I / h i s t s t r o n g/ l oo s e A r c - T o r n e t r a s k . M e l v i n . . . t r ee T R W P I / h i s t s t r o n g/ l oo s e A r c - Y a m a li a . B r i ff a . . t r ee T R W P I / h i s t s t r o n g/ l oo s e A s i a - C e n t r a l C h i n a . W a n g . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - C hu . . S i h a il o n g w a n l a k e . . l a k e s e d i m e n t a l k e n o n e P I / h i s t s t r o n g/ l oo s e A s i a - E a s t C h i n a . W a n g . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - E a s t C h i n a r e g . W a n g . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - F u j i a n a nd T a i. W a n g . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - G u a n g d o n g . Z h e n g . . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - G u a n g d o n ga nd . Z h a n g . . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - H un a n J i a n g s u . Z h a n g . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - K un a s h i r I s l a . D e m ez h k o . . h y b r i dh y b r i d P I / h i s t s t r o n g/ l oo s e A s i a - L o w e rr e a c h e s . Z h a n g . . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - M i dd l e r e a c h e . Z h a n g . . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - S O K G S O . M un z . . . m a r i n e s e d i m e n t f o r a m i n i f e r a P I / h i s t s t r o n g/ l oo s e A s i a - S o u r t h a nd M i d . D e m ez h k o . . b o r e h o l e b o r e h o l e P I / h i s t l oo s e D N a m e L a t L o n E l e v m a s l A r c h i v e P r o xy c u t o ff s e l ec t i o n A s i a - S o u t h C h i n a . W a n g . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A s i a - Z h e j i a n ga nd F . Z h a n g . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e A u s - M t R e a d . C oo k . − . . t r ee T R W P I / h i s t s t r o n g/ l oo s e A u s - O r o k o . C oo k . − . . t r ee T R W P I / h i s t s t r o n g/ l oo s e E u r - C e n t r a l a nd E a . P l a . . . l a k e s e d i m e n t c h r y s o ph y t e P I / h i s t l oo s e E u r - C e n t r a l E u . D o b r o v o l n y . d o c u m e n t s D o c u m e n t a r y P I / h i s t s t r o n g/ l oo s e E u r - C oa s t o f P o r t u . A b r a n t e s . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t s t r o n g/ l oo s e E u r - F i nn i s h L a k e s . H e l a m a . . t r ee M X D P I / h i s t s t r o n g/ l oo s e E u r - L a k e S il v a p l a . L a r o c q u e - T o b l e r . . . l a k e s e d i m e n t c h i r o n o m i d P I / h i s t s t r o n g/ l oo s e E u r - L a k e S il v a p l a . T r a c h s e l. . . l a k e s e d i m e n t r e fl ec t a n ce P I / h i s t s t r o n g/ l oo s e E u r - N o r t h I ce l a nd . R a n . . − . − m a r i n e s e d i m e n t d i a t o m P I / h i s t s t r o n g/ l oo s e E u r - S ee b e r g s ee . L a r o c q u e - T o b l e r . . . l a k e s e d i m e n t m i d g e P I / h i s t s t r o n g/ l oo s e E u r - S t o c k h o l m . L e i j o nhu f v ud . . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e E u r - T a lli nn . T a r a nd . . . d o c u m e n t s h i s t o r i c P I / h i s t s t r o n g/ l oo s e NA m - B a s i n P . − . − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e NA m - C l e a r P . − − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e NA m - C o n r o y L . − . − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e NA m - D a r k L . − . − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e NA m - H e ll K t . − . − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e NA m - L a k e M i n a45 . − . − l a k e s e d i m e n t p o ll e n P I / h i s t s t r o n g/ l oo s e NA m - L C l o ud s − − l a k e s e d i m e n t p o ll e n P I / h i s t s t r o n g/ l oo s e NA m - L i tt l e P i n e L . − . − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e NA m - L N o i r . − . − l a k e s e d i m e n t p o ll e n P I / h i s t s t r o n g/ l oo s e NA m - R ub y L . − . − l a k e s e d i m e n t p o ll e n P I / h i s t l oo s e O k L R - A l b o r a nS e a - TT R - B . N i e t o - M o r e n o . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - A l b o r a nS e a - TT R - B . N i e t o - M o r e n o . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - A r a b i a nS e a . D oo s e - R o li n s k i. . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - C a p e G h i r . K i m . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - C a p e G h i r . M c G r e go r . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - C a p e H a tt e r a s . C l r o u x . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t s t r o n g/ l oo s e O k L R - C a r i a c o B a s i n . B l a c k . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t s t r o n g/ l oo s e O k L R - C h il e a n M a r g i n . L a m y . − − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - D r y T o r t u ga s . L und . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e D N a m e L a t L o n E l e v m a s l A r c h i v e P r o xy c u t o ff s e l ec t i o n O k L R - D r y T o r t u ga s A . L und . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - E a s t e r n T r o p i c a l N o r t h A t l a n t i c . K uhn e r t . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - E m e r a l d B a s i n . K e i g w i n . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - F e n i D r i f t R i c h t e r . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - F i s k B a s i n . R i c h e y . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t s t r o n g/ l oo s e O k L R - G a rr i s o n B a s i n . R i c h e y . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a h i s t l oo s e O k L R - G r e a t B a h a m a B a n k . L und . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - G r e a t B a h a m a b A t l a n t i c L und M C . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - G r e a t B a rr i e r . H e nd y . − . . c o r a l C o r a l S r / C a P I / h i s t s t r o n g/ l oo s e O k L R - G u l f o f G u i n e a . W e l d e a b . . . m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - J a c a f F j o r d . S e pu l v e d a . − . − − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - K u r o s h i o C u rr e n t . I s o n o . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - M a k a ss a r S t r a i t - M D - . N e w t o n . . . m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t s t r o n g/ l oo s e O k L R - M a k a ss a r S t r a i t . L i n s l e y . − . . m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - M a k a ss a r S t r a i t . O pp o . − . . m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - M D . C a l v o . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - M i n o r c a . M o r e n o . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - N o r t h I ce l a nd . S i c r e . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t s t r o n g/ l oo s e O k L R - N W P a c i fi c . H a r a d a . . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - O k i n a w a T r o u g h . W u . . . m a r i n e s e d i m e n t T E X P I / h i s t l oo s e O k L R - P h ili pp i n e s - M D - . S t o tt . . . m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e O k L R - P i g m y B a s i n . R i c h e y . . − . − m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t s t r o n g/ l oo s e O k L R - S a n t a B a r b a r a . Z h ao . . − − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t s t r o n g/ l oo s e O k L R - S o u t h A t l a n t i c . L e du c . − . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - S o u t h C h i n a S e a . Z h ao . . . m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - S o u t h e r n C h il e . M o h t a d i. − . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t l oo s e O k L R - S o u t h I ce l a nd . S i c r e . . − . − m a r i n e s e d i m e n t a l k e n o n e P I / h i s t s t r o n g/ l oo s e O k L R - Sub T r o p i c a l E a s t e r n A t l a n t i c . d e M e n o c a l. . − . − m a r i n e s e d i m e n t p l a n k t o n i c f o r a m i n i f e r a P I / h i s t l oo s e O k L R - W e s t e r nS v a l b a r d . Sp i e l h ag e n . . . m a r i n e s e d i m e n t p l a n k t o n i c f o r a m i n i f e r a P I / h i s t l oo s e O k L R - W e s t Sp i t z b e r g e n . B o nn e t . . m a r i n e s e d i m e n t d y n o c i s t M A T P I / h i s t l oo s e O k L R M a k a ss a r S t r a i t - M D - . N e w t o n . − . . m a r i n e s e d i m e n t f o r a m M g/ C a P I / h i s t l oo s e S A m - L ag un a A c u l e o . v o n G un t e n . − . − . − l a k e s e d i m e n t r e fl ec t a n ce P I / h i s t s t r o n g/ l oo s e S A m - L ag un a C h e p i c a l. d e J o n g . − . − . − l a k e s e d i m e n t r e fl ec t a n ce P I / h i s t s t r o n g/ l oo s e S A m - L ag un a E s c o nd i d a . E l b e r t . − . − . − l a k e s e d i m e n t B S i P I / h i s t s t r o n g/ l oo s e AGES2k FIG. S1. Evolution of global mean temperature over the Common Era from model simulationsand observation-based data as in Table S1. Anomalies are given with respect to the referenceperiod 1961-1990 (HadCM3 LM1:1800-1850) and as a running average of five years. A O D CrowleyGaoToohey R F ( W / m ) DelaygueMuscheler SteinhilberVieira Wang c on c . ( pp m ) Schmidt R F ( W / m ) Berger time (yr CE) s o l a r v o l c an i c C O i n s o l a t i on ° N R F ( W / m ) Fröhlich(a) (b) c on c . ( pp m ) Keelingtime (yr CE) s o l a r C O FIG. S2. (a) Reconstruction of climate drivers over the Common Era used to estimate the PSDof radiative forcing (RF). Labels indicate the data reference as given in Table S2. (b) Additionalobservational data for solar and CO forcing used to obtain high-frequency spectral estimates.Highly resolved insolation changes due to the diurnal and annual cycle are not shown here.
54 655 656469 514 541 550 554352 355 357 359 466152 273 342 343 35049 82 83 85 88 T e m pe r a t u r e ano m a l y ( K ) FIG. S3. Temperature anomalies from proxy records used to estimate local temperature variability.The labels give the ID from the PAGES2k database as indicated by a star ( ? ) in Table S3. IG. S4. As Figure S3, but showing the PSD of the proxy records. The shaded area indicates theperiods between 10 and 200 years, used to estimate the scaling relationship. Scaling exponents forspectra that do not cover the full period were estimated on their corresponding timescales, but atleast between 20 and 200 years. IG. S5. Power spectral densities of solar (a), CO (b), and volcanic (c) from multiple data setsand their mean spectra (black). Low- and high-frequency tails, likely prone to spectral biases, wereomitted in the average.FIG. S6. PSD for global mean and mean of local spectra of surface air temperature from lastmillennium runs and their PI control (cntl). IG. S7. Joint PSD (dashed orange line) from the mean spectrum of a highly (hr) and lowerresolved time series on the example of surface air temperature from the CESM-LME 1 run. Toimprove computational efficiency, we formed this joint PSD for ERA5, MPI-M LM, and CESM-LME 1. The high-frequency component was estimated from data at hourly resolved temperatureseries from 1981 to 1990 CE since this is a part of the climatic period that is shared among thethree data sets. int signalseriessurrogate PSD
FIG. S8. Step-wise estimation of uncertainty on β -scaling relation for irregularly sampled proxyrecords on the example of the PSD from the “NAm-LakeMina”-record. When interpolated toits mean temporal resolution ∆ t , the raw signal has a scaling relation β = 1 .
31 (orange dashedline). The uncertainties were calculated from 100 surrogates that are random time series withapproximately the same power-law scaling β (black dotted line). The PSD is altered after formingthe block-average of the surrogate time series and interpolation to ∆ t (blue solid line). Linearregression on this spectrum (black dashed line) gives the scaling coefficient β int of the surrogateseries. uto ff : FIG. S9. Global mean of local PSD from proxy records subject to sampling from the PAGES2kdata base and global warming. Dashed and solid lines refer to the mean spectra of local proxyrecords cut at 1850 CE (PI) and 2000 CE (hist) respectively. The plot compares strong (blue) andloose (orange) criteria. Loose criteria are described in the data section of the main manuscript andwere used to calculate the mean of local spectra from reconstructions. For the stronger criteria werequire the mean temporal resolution ( h t i +1 - t i i ) ≤ t N - t ) ≥ N ≥
30. Hiatus are tolerated up to max( t i +1 - t i ) ≥
40 years.
ProxyCESM 1 past2kCESM-LME 1MPI-M LMECH5/MPIOM-p6kTraCE-21k IPSL-p6k -0.2 0.0 0.2 β hist − β PI FIG. S10. Normalized density plots of the deviations of the spectral exponent β extracted fromtime series up to 2020 CE ( β hist ) and time series up to 1850 CE ( β PI ). The densities were computedfrom the local power-law scaling used to estimate the statistical agreement between models anddata (Table S3).
55 (-70.5:-32.3) 654 (-70.9:-33.9) 342 (145.5:-41.8) 343 (170.3:-43.2) 656 (-71.8:-45.5)469 (-75.1:45.8) 350 (-8.9:41.1) 550 (-120:34.2) 152 (120:34) 273 (65.9:24.8) 514 (-64.8:10.8)352 (28.3:62) 49 (-146.6:61) 554 (-27.9:57.5) 355 (9.8:46.5) 359 (7.5:46.1) 466 (-95.5:45.9)83 (-68.8:69.9) 85 (19.6:68.3) 88 (68:66.8) 541 (-17.4:66.5) 357 (-17.7:66.5) 82 (5.3:63.8)
ProxyCESM 1 past2kCESM-LME 1CESM-LME 1 (cntl)MPI-M LMMPI-M (cntl)TraCE-21k-ORBTraCE-21kHadCM3 LM1ECH5/MPIOM-p6kIPSL-p6k -1 0 1 2 3 -1 0 1 2 3 -1 0 1 2 3 -1 0 1 2 3 -1 0 1 2 3 -1 0 1 2 3 scaling coefficient β modeltreelake/marinesedimentdocumentsProxyCESM 1 past2kCESM-LME 1CESM-LME 1 (cntl)MPI-M LMMPI-M (cntl)TraCE-21k-ORBTraCE-21kHadCM3 LM1ECH5/MPIOM-p6kIPSL-p6kProxyCESM 1 past2kCESM-LME 1CESM-LME 1 (cntl)MPI-M LMMPI-M (cntl)TraCE-21k-ORBTraCE-21kHadCM3 LM1ECH5/MPIOM-p6kIPSL-p6kProxyCESM 1 past2kCESM-LME 1CESM-LME 1 (cntl)MPI-M LMMPI-M (cntl)TraCE-21k-ORBTraCE-21kHadCM3 LM1ECH5/MPIOM-p6kIPSL-p6k FIG. S11. Local scaling coefficients and its confidence intervals from simulations and paleoclimatedata. Each panel denotes the ID within the PAGES2k database and its coordinates (longitude ◦ E,latitude ◦ N). Symbols indicate the archive of the proxy record or whether the scaling exponentbelongs to a model simulation or not. Red shading refers to the unstable regime ( β > β < adCM3 LM1(a) CESM 1 past2k(b) CESM-LME 1(c)CESM-LME 1 (cntl)(d) MPI-M LM(e) MPI-M LM (cntl)(f)ECH5/MPIOM-p6k(g) TraCE-21k-ORB(h) TraCE-21k(i) IPSL-p6k(j) std. Δβ FIG. S12. Standard deviation of the least-square regression of local scaling coefficient β for allmodel simulations. The goodness of fit is generally lower for models with a comparatively lowtemporal resolution (Trace21k-ORB) or coverage (HadCM3 LM1). Furthermore, there is a slightlyincreased deviation in areas of active modes of variability, such as the ENSO region ((b), (c), (e),(j)). This could be due to the fact, that spectral peaks describe these (quasi-)periodic signals betterthan power-law scaling.
1] Y. Zhong, A. Jahn, G. H. Miller, and A. Geirsdottir, Asymmetric Cooling of the Atlanticand Pacific Arctic During the Past Two Millennia: A Dual Observation-Modeling Study,Geophysical Research Letters , 12,412 (2018).[2] J. H. Jungclaus et al., The PMIP4 contribution to CMIP6 – Part 3: The last millennium,scientific objective, and experimental design for the PMIP4 simulations, Geoscientific ModelDevelopment , 4005 (2017).[3] B. L. Otto-Bliesner, E. C. Brady, J. Fasullo, A. Jahn, L. Landrum, S. Stevenson, N. Rosen-bloom, A. Mai, and G. Strand, Climate variability and change since 850 ce an ensembleapproach with the community earth system model, Bulletin of the American MeteorologicalSociety 10.1175/BAMS-D-14-00233.1 (2016).[4] G. A. Schmidt et al., Climate forcing reconstructions for use in PMIP simulations of the lastmillennium (v1.0), Geoscientific Model Development , 33 (2011).[5] J. H. Jungclaus et al., Climate and carbon-cycle variability over the last millennium, Climateof the Past 10.5194/cp-6-723-2010 (2010).[6] J. C. B¨uhler, C. Roesch, M. Kirschner, L. Sime, M. D. Holloway, and K. Rehfeld, Comparisonof the oxygen isotope signatures in speleothem records and iHadCM3 model simulations forthe last millennium, Climate of the Past Discussions , 1 (2020).[7] A. P. Schurer, S. F. B. Tett, and G. C. Hegerl, Small influence of solar variability on climateover the past millennium, Nature Geoscience , 104 (2014).[8] P. Braconnot, D. Zhu, O. Marti, and J. Servonnat, Strengths and challenges for transient Mid-to Late Holocene simulations with dynamical vegetation, Climate of the Past , 997 (2019).[9] Z. Liu, Transient simulation of last deglaciation with a new mechanism for Bølling–Allerødwarming, Science , 10.1126/science.1171041 (2009).[10] A. L. Berger, Long-term variations of daily insolation and Quaternary climatic changes., Jour-nal of Atmospheric Sciences 10.1175/1520-0469(1978)035¡2362:ltvodi¿2.0.co;2 (1978).[11] F. Joos and R. Spahni, Rates of change in natural and anthropogenic radiative forcing overthe past 20,000 years, Proceedings of the National Academy of Sciences , 1425 LP (2008).[12] N. Fischer and J. H. Jungclaus, Evolution of the seasonal temperature cycle in a transientHolocene simulation: Orbital forcing and sea-ice, Climate of the Past 10.5194/cp-7-1139-2011(2011).[13] H. Hersbach, B. Bell, P. Berrisford, and S. Hirahara, The ERA5 global reanalysis, QuarterlyJournal of the Royal Meteorological Society , 1999 (2020).[14] C. P. Morice, J. J. Kennedy, N. A. Rayner, and P. D. Jones, Quantifying uncertainties in globaland regional temperature change using an ensemble of observational estimates: The Had-CRUT4 data set, Journal of Geophysical Research Atmospheres , 10.1029/2011JD017187(2012).[15] PAGES 2k Consortium., Consistent multidecadal variability in global temperature reconstruc-tions and simulations over the Common Era, Nature Geoscience , 643 (2019).[16] G. A. Schmidt et al., Climate forcing reconstructions for use in PMIP simulations of the LastMillennium (v1.1), Geoscientific Model Development , 185 (2012).[17] T. J. Crowley and M. B. Unterman, Technical details concerning development of a 1200 yr roxy index for global volcanism, Earth System Science Data , 187 (2013).[18] C. Gao, A. Robock, and C. Ammann, Volcanic forcing of climate over the past 1500 years:An improved ice core-based index for climate models, Journal of Geophysical Research: At-mospheres , 10.1029/2008JD010239 (2008).[19] M. Toohey and M. Sigl, Reconstructed volcanic stratospheric sulfur injections and aerosoloptical depth, 500 BCE to 1900 CE, version 2, World Data Center for Climate (WDCC) atDKRZ 10.1594/WDCC/eVolv2k v2 (2017).[20] G. Delaygue and E. Bard, An Antarctic view of Beryllium-10 and solar activity for the pastmillennium, Climate Dynamics , 2201 (2011).[21] R. Muscheler, F. Joos, J. Beer, S. A. M¨uller, M. Vonmoos, and I. Snowball, Solar activityduring the last 1000yr inferred from radionuclide records, Quaternary Science Reviews , 82(2007).[22] F. Steinhilber, J. Beer, and C. Fr¨ohlich, Total solar irradiance during the Holocene, Geophys-ical Research Letters , 10.1029/2009GL040142 (2009).[23] L. E. Vieira and S. K. Solanki, Evolution of the solar magnetic flux on time scales of years tomillenia, Astronomy and Astrophysics 10.1051/0004-6361/200913276 (2010).[24] N. A. Krivova, L. Balmaceda, and S. K. Solanki, Reconstruction of solar total irradi-ance since 1700 from the surface magnetic flux, Astronomy and Astrophysics 10.1051/0004-6361:20066725 (2007).[25] Y. Wang, J. L. Lean, and N. R. Sheeley, Jr., Modeling the Sun’s Magnetic Field and Irradiancesince 1713, The Astrophysical Journal 10.1086/429689 (2005).[26] C. Fr¨ohlich, Solar irradiance variability since 1978: Revision of the PMOD composite duringsolar cycle 21, Space Science Reviews , 53 (2006).[27] C. D. Keeling, R. B. Bacastow, and A. E. Bainbridge, Atmospheric carbon dioxide variationsat Mauna Loa Observatory, Hawaii, TELLUS 10.3402/tellusa.v28i6.11322 (1976).[28] M. Crucifix, Palinsol - Package (R) (2016)., 53 (2006).[27] C. D. Keeling, R. B. Bacastow, and A. E. Bainbridge, Atmospheric carbon dioxide variationsat Mauna Loa Observatory, Hawaii, TELLUS 10.3402/tellusa.v28i6.11322 (1976).[28] M. Crucifix, Palinsol - Package (R) (2016).