Supplementary Material for “Characterizing decadal to centennial
variability in the equatorial Pacific during the last millennium”

May 28, 2013

Origins of low-frequency variability
In this document we present several additional analyses to diagnose the origin of low-frequency(dec-cen) variability in long unforced and forced climate model simulations, as well as in the reconstructions developed in Emile-Geay et al. [2013a] and further explored in Emile-Geay et al. [2013b](herein referred to as EG13a and EG13b). We also examine the underlying spectral characteristicsof the paleoclimate archives used in the EG13a,b reconstructions. As similar analyses were performed in EG13b, we will begin by summarizing several findings from the earlier study that are
most relevant to our work:• Other recent reconstructions employing a similar methodology [Wilson et al., 2010; Mannet al., 2009] are generally within the uncertainty estimates of EG13b.
• Reconstructed low-frequency variability is more sensitive to the target dataset (e.g., ERSSTv3,
HadSST2i, or Kaplan) than to most other methodological details.
• An alternative method to RegEM (composite plus scale; CPS) was found to be less successful
at recovering low-frequency variability from the multi-proxy network. This differs from theresult of Wilson et al. [2010], which showed very little difference between CPS and RegEM oncentennial timescales. As explained in EG13b, these discrepancies arise because the proxies
used in Wilson et al. [2010] were detrended (the EG13a,b reconstructions were not), whichconsequently eliminates any information in the low-frequency spectral range and makes the1results insensitive to reconstruction method. In the EG13a,b reconstructions, the CPS-based
reconstructions are less energetic on long timescales than the ones employing RegEM, butthey are also demonstrably less robust.
• “Jackknife” experiments to establish uncertainty estimates suggest that no one proxy is driving the result (EG13b, Figure 8), but that the low-frequency variability is in a general sense“coherent” across the network.
• Low-frequency power originates from different proxy records through time (EG13b, Figures 5& 6), which is a feature of the network and the methods employed by EG13a,b. For example,
reconstructions developed using both RegEM and CPS generally agree with each other on
dec-cen timescales, ostensibly meaning that these two different methods are finding the sharedsignal through time. Hence, the provenance of low-frequency power in the reconstruction is
not tied to any one proxy type or region.

Other reconstructions
As in Emile-Geay et al. [2013b], we compare the EG13a,b reconstructions with two recently published results that employ similar methodologies and use some of the same underlying proxies
records (Figure S1). Amplitudes on timescales longer than ∽ 100 years are far weaker in the Wilson et al. [2010] reconstruction, likely reflecting detrending and removal of low-frequency variationsin that study. The amplitude of variability in the [Mann et al., 2009] study is comparable toEG13a,b on centennial timescales, and necessarily much weaker on interannual timescales because
the reconstruction targeted decadally-smoothed fluctuations of global surface temperature (andused low-pass filtered predictors to do so).

Other CMIP5 models
Basic information about the 24 CMIP5 models whose pre-industrial control (i.e., “piControl”)
simulations were used in Figure 2d is summarized in Table S1. Spectra from each individual model
are compared with the LIM in Figure S2. Likewise, Table S2 lists the models with last millennium(“past1000”) simulations, and Figure S3 compares each of these models with the LIM. As in Figure2

3 of the main text, only the period of overlap (1000-1850 CE) of each simulation was used to ensurethat inter-model comparisons were valid.
Results obtained by comparing each available piControl and past1000 simulation with the LIM
support the main conclusions of the text (Figures S2 and S3). Specifically: (1) the amplitude of
dec-cen variations in the control runs of the CMIP5 archive are generally consistent with the linear,stochastically-forced paradigm of tropical Pacific variability; and (2) introducing external forcingsincreases the variance of last millennium simulations on dec-cen timescales, generally (although notentirely, in some cases) beyond what the LIM alone produces. We emphasize, however, that theforced responses of these models do not agree with the reconstruction, as was shown for CCSM4
in the main text. This finding therefore further argues in favor of the idea that even dec-cen
fluctuations may be internally generated or arise as a weak response to external forcings.
On the other hand, evidence presented in EG13b suggests that there may indeed be a significant
role for solar variability to induce dynamical changes on ∽ 200 year timescales in the equatorial
Pacific. Correlations in table S3 support this idea, showing that two out of the three EG13a,b
reconstructions are significantly negatively correlated with reconstructed solar variability. In contrast, the CCSM4 simulation is positively correlated with the solar forcing term, as are most of theother CMIP5 models. Further investigation of the origin and dynamics of these ultra-low frequency
fluctuations in the CMIP5 and paleoclimate archives should be the subject of future work.

The robustness of centennial-scale variability in CCSM4Spectral estimates calculated using an alternative technique (the Blackman-Tukey method) werevery similar to those shown in Figure 2 (Figure S4). However, NINO3.4 spectra calculated fromCCSM4 over several subsets of the full 850-2005 interval highlight the contributions of 20th forcings to the overall spectrum (Figure S6a). These other segments were selected to maximize thetotal amount of forced pre-industrial variability from the CCSM4 simulation (850-1850 CE), theoverlap with other CMIP5 models (1000-1850 CE), and the pre-indusrial overlap with the EG13a,b
reconstructions (1150-1850 CE). It is quite clear that considering only the period of the CCSM4
simulation from before 1850 reduces a substantial portion of the total dec-cen variability from thepower spectrum, and indeed suggests that variations on long timescales in this run may in fact be3

consistent with the LIM of late 20th century instrumental data (red line; CCSM41000−1850 ). Importantly, however, all subsets of the last millennium simulation are above a LIM representation of
CCSM4’s control run, meaning that variance on these timescales in the forced run is very unlikely
to have been generated by internal processes (Figure S6).

The Emile-Geay et al. proxy network
Power spectra from each individual proxy record used in EG13a,b to reconstruct the HadSST2i
NINO3.4 time series are shown in Figure S7. The underling spectral slopes (i.e., power-laws) ofeach record are shown in Figures S8 and S9. Collectively, these three figures illustrate that no
obvious relationship exists between spectral density and proxy type or geography, except that treering spectra exhibit the weakest dec-cen amplitudes and are characterized by the lowest power-laws(Figures S7e and S8). Further, low-frequency variability in the reconstruction is not tied to any
one site, which can plainly be seen in Figure S10.

“Hasselmann” model of forced responses in CCSM4Here we consider the power spectrum of a simple system that is approximately linearly forced byradiative fluctuations derived from the CCSM4 last millennium simulation. Certain aspects of sucha system can be understood in terms of a “Hasselmann” type of model [e.g., Thompson et al., 2009]:

Cd
TRAD (t) = −TRAD (t) + αF (t),dt(1)

where TRAD (t) is a time series of radiatively forced temperature anomalies, C is a constant representing heat capacity, α relates this forcing to a temperature response, and F (t) is the time-evolving
radiative forcing term.
We employ this linear approximation of forced changes in CCSM4 as a diagnostic tool to characterize the underlying shape of the externally forced power spectrum. To do so, we use the downwelling solar flux anomaly (FSNTOA) time series in Figure 1e for F (t). Parameters C and α wereselected to minimize the difference from 850-1850 CE between the forced temperature anomalies(TRAD ) and CCSM4 NINO3.4 smoothed using an 11-year Gaussian filter to de-emphasize ENSO

4

variability. Although we used values of 0.9 (107 J m−2 K −1 ) for C and 0.06 (K W−1 ) for α, wefound that a range of values for C (0.5 to 0.95) and α (0.01 to 0.16) produced spectra that werequite similar in shape.
The time series of our “Hasselmann” model of temperature anomalies (TRAD ) forced by CCSM4
downwelling solar flux anomalies at the top of the atmosphere (FSNTOA) is shown in Figure S11a.
Also shown are the large (e.g., greater than Pinatubo) volcanic injections used to force the model,
as in Figure 1 of the main text. Large volcanic eruptions induce a 1 to 2 o C temperature anomaly
in both the NINO3.4 region and in the TRAD model – a result that agrees well with the global
temperature response to volcanic eruptions in the CCSM4 last millennium simulation [Landrum
et al., 2012].
On dec-cen timescales, the power spectrum of NINO3.4 from CCSM4 and TRAD are in very good
agreement (Figure S11b). This finding lends credibility to the idea that low-frequency variability inthe model is an approximately linear response to radiative forcing. Moreover, because the highest
magnitude radiative anomalies are clearly tied to volcanic events, it follows that most of the lowfrequency variance in CCSM4’s tropical Pacific region arises from this forcing component.

References
Emile-Geay, J., K. M. Cobb, M. E. Mann, and A. T. Wittenberg, Estimating Central Equatorial
Pacific SST variability over the Past Millennium. Part 1: Methodology and Validation, Journalof Climate, doi:10.1175/JCLI-D-11-00510.1, 2013a.
Emile-Geay, J., K. M. Cobb, M. E. Mann, and A. T. Wittenberg, Estimating Central Equatorial
Pacific SST variability over the Past Millennium. Part 2: Reconstructions and Implications,
Journal of Climate, doi:10.1175/JCLI-D-11-00511.1, 2013b.
Landrum, L., B. L. Otto-Bliesner, E. R. Wahl, A. Conley, P. J. Lawrence, N. Rosenbloom, and
H. Teng, Last Millennium Climate and Its Variability in CCSM4, Journal of Climate, 26 (4),1085–1111, doi:10.1175/JCLI-D-11-00326.1, 2012.
Mann, M. E., Z. Zhang, S. Rutherford, R. S. Bradley, M. K. Hughes, D. Shindell, C. Ammann,

5G. Faluvegi, and F. Ni, Global Signatures and Dynamical Origins of the Little Ice Age andMedieval Climate Anomaly, Science, 326 (5957), 1256–1260, doi:10.1126/science.1177303, 2009.
Steinhilber, F., J. Beer, and C. Froehlich, Total solar irradiance during the Holocene, GeophysicalResearch Letters, 36, L19,704, doi:10.1029/2009GL040142, 2009.
Thompson, D. W. J., J. M. Wallace, P. D. Jones, and J. J. Kennedy, Identifying Signatures of
Natural Climate Variability in Time Series of Global-Mean Surface Temperature: Methodology
and Insights, Journal of Climate, 22 (22), 6120–6141, doi:10.1175/2009JCLI3089.1, 2009.
Vieira, L. E. A., and S. K. Solanki, Evolution of the solar magnetic flux on time scales of years tomillenia, Astronomy & Astrophysics, 509, doi:10.1051/0004-6361/200913276, 2010.
Wilson, R., E. Cook, R. D’Arrigo, N. Riedwyl, M. N. Evans, A. Tudhope, and R. Allan, Reconstructing ENSO: the influence of method, proxy data, climate forcing and teleconnections,
Journal of Quaternary Science, 25 (1), 62–78, doi:10.1002/jqs.1297, 2010.

6ModelACCESS1-0
CCSM4
CCSM4
CCSM4
CCSM4
CNRM-CM5CSIRO-Mk3-6-0
CanESM2
FGOALS-g2
GFDL-CM3
GFDL-ESM2G
GFDL-ESM2M
GISS-E2-H
GISS-E2-H
GISS-E2-H
GISS-E2-R
GISS-E2-R
GISS-E2-R
GISS-E2-R
GISS-E2-R
HadGEM2-CC
HadGEM2-ES
IPSL-CM5A-LR
IPSL-CM5A-MR
MIROC-ESM-CHEM
MIROC-ESM
MIROC4h
MIROC5MPI-ESM-LR
MRI-CGCM3
NorESM1-M
bcc-csm1-1
inmcm4run
r1i1p1
r1i1p1r2i1p1
r3i1p1
r4i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p2
r1i1p3
r1i1p141
r1i1p142
r1i1p1
r1i1p2
r1i1p3
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1
r1i1p1Run Length (years)
250
501
156
120
50
850
500
996
900
500
500
500
780
531
531
1163
100
850
531
531
240
577
1000
300
255
531
80
670
1000
500
501
500
500

Variable
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts
ts

Table S1: Summary of pre-industrial control (piControl) runs whose NINO3.4 spectra were shownin Figure 3 of the main text. The variable “ts” corresponds to “skin” temperature (e.g., SST
over the open ocean). All simulations are available online through the Earth System Grid:
http://pcmdi9.llnl.gov/esgf-web-fe/ .

7ModelCCSM4GISS-E2-R
GISS-E2-R
GISS-E2-R
GISS-E2-R
GISS-E2-R
GISS-E2-R
GISS-E2-R
IPSL-CM5A-LR
MIROC-ESM
MPI-ESM-P
bcc-csm1-1

run
r1i1p1
r1i1p121
r1i1p122
r1i1p123
r1i1p124
r1i1p125
r1i1p126
r1i1p127
r1i1p1
r1i1p1
r1i1p1
r1i1p1Run Length (years)
1001
1001
1001
1001
1001
1001
1001
1001
1001
1000
1000
1151

Variable
tas
tas
tas
tas
tas
tas
tas
tas
tas
tas
tas
tas

Table S2: Same as Table S1, but for the runs with available last millennium (“past1000”) simulations. Note that here “surface air temperature” (tas) was used instead of ts because it was available
from a greater number of models.

ERSSTv3
HadSST2i
KaplanCCSM4GISS-E2-R (r1i1p121)
GISS-E2-R (r1i1p122)
GISS-E2-R (r1i1p123)
GISS-E2-R (r1i1p124)
GISS-E2-R (r1i1p125)
GISS-E2-R (r1i1p126)
GISS-E2-R (r1i1p127)
IPSL-CM5A-LR (r1i1p1)
MIROC-ESM (r1i1p1)
MPI-ESM-P (r1i1p1)
bcc-csm1-1 (r1i1p1)

Steinhilber
-0.46
-0.53
-0.66
0.44
0.30
0.30
0.37
0.43
0.43
-0.12
0.36
0.45
-0.19
0.52
0.60

Viera and Solanski
-0.39
-0.58
-0.54
0.30
0.04
0.03
0.30
0.27
0.16
0.10
0.16
0.43
-0.10
0.36
0.48

Table S3: Correlation coefficients between low-pass (>100 year) filtered NINO3.4 time series andsolar forcings. Bold indicates significance at the 95% confidence limit. Correlations were calculated
over the time period 1150-1900, and significance (bold) was estimated using an effective sample
size of 30. The ERSSTv3, HadSST2i, and Kaplan NINO3.4 products are the reconstruction of
EG13a,b, and the reconstructions of solar activity are described in Steinhilber et al. [2009] andVieira and Solanki [2010]. As the Vieira and Solanki [2010] time series was used to force CCSM4,
we only report the correlations between this reconstruction of solar activity the model.

8

CCSM4 (r1i1p1)
GISS-E2-R (r1i1p121)
GISS-E2-R (r1i1p122)
GISS-E2-R (r1i1p123)
GISS-E2-R (r1i1p124)
GISS-E2-R (r1i1p125)
GISS-E2-R (r1i1p126)
GISS-E2-R (r1i1p127)
IPSL-CM5A-LR (r1i1p1)
MIROC-ESM (r1i1p1)
MPI-ESM-P (r1i1p1)
bcc-csm1-1 (r1i1p1)

ERSSTv3
-0.15
0.03
-0.26
0.23
0.09
-0.21
0.69
-0.02
-0.4
-0.62
-0.15
-0.18

HadSST2i
-0.18
0.11
-0.08
0.17
0.06
-0.14
0.43
0.08
-0.4
-0.56
-0.10
-0.37

Kaplan
-0.33
-0.13
-0.38
0.13
-0.11
-0.36
0.67
-0.16
-0.48
-0.42
-0.31
-0.37

Table S4: Correlation coefficients between low-pass filtered (> 100 year) NINO3.4 time series from
CMIP5 models (rows) and the reconstructions (columns). Correlations were calculated over thetime period 1150-1900, and significance (bold) was estimated using an effective sample size of 30.

9

a)
LIM
JEG13
W10 (COA)
W10 (TEL)

1

Power

10

0

10

-1

10

1000

100 5020
Period (years)

10

5

2

b)
LIM
JEG13
M09 (NINO3)

1

Power

10

0

10

-1

10

1000

100 5020
Period (years)

10

5

2

Figure S1: Power spectra of several recent reconstructions of tropical Pacific variability. As in
Figure 2 of the main text, the EG13a,b reconstructions are shown in dark green and the distribution
of NINO3.4 spectra generated by the LIM is shown in gray. Additional reconstructions include:a) Wilson et al. [2010]; and b) Mann et al. [2009]. In Wilson et al. [2010], two reconstructions
were developed to investigate differences between a network of teleconnected predictors (TEL), andpredictors in the “center of action” (COA). The time series generated by Mann et al. [2009] wasproduced by averaging the gridded global temperature reconstruction in that study for the NINO3
region (and because the data were decadally-filtered, the spectrum is only meaningful on the 10
year and longer timescales shown in the figure).

10

ACCESS1-0 (1)CCSM4 (4)

CNRM-CM5 (1)CSIRO-Mk3-6-0 (1)CanESM2 (1)FGOALS-g2 (1)GFDL-CM3 (1)GFDL-ESM2G (1)GFDL-ESM2M (1)GISS-E2-H (3)GISS-E2-R (5)HadGEM2-CC (1)IPSL-CM5A-LR (1)IPSL-CM5A-MR (1)

MIROC-ESM (1)

MIROC4h (1)MIROC5 (1)MPI-ESM-LR (1)MRI-CGCM3 (1)NorESM1-M (1)bcc-csm1-1 (1)inmcm4 (1)1000 100
10 2
Period (years)

1000 100
10 2
Period (years)

1000 100
10 2
Period (years)

1000 100
10 2
Period (years)

power

0

10

-1

10

power

0

10

-1

10

power

0

10

-1

10

HadGEM2-ES (1)
power

0

10

-1

10

MIROC-ESM-CHEM (1)
power

0

10

-1

10

power

0

10

-1

10

Figure S2: Power spectra for each of the control runs from each of the 24 models shown in Figure
2 of the main text. Gray shading indicates the upper and lower bounds of the LIM. Multiple runs
from the same model are plotted on the same panel in slightly different colors of aqua, and thenumber of runs is indicated in parentheses next to the model name.

11

CCSM4 (1)GISS-E2-R (7)

0

0

power

10

10

-1

-1

10

1000

10

100

10

2

1000

IPSL-CM5A-LR (1)
0

0

power

-1

-1

10

100

10

2

1000

MPI-ESM-P (1)

100

10

2

100
10
Period (years)

2bcc-csm1-1 (1)

0

0

10
power

2

10

10

10

-1

-1

10

1000

10

MIROC-ESM (1)

10

1000

100

10

100
10
Period (years)

2

1000

Figure S3: Same as Figure S2, but for the available last millennium (“past1000”) runs. Again,
different runs are shown with varying shades of aqua, and the corresponding control run for eachmodel (if available) is shown in red.

12

a)

b)

1

LIM
CCSM4
Instrumental

1

10
Power

Power

10

0

10

-1

10
100 50 20 10
Period (years)

5

2

c)

1000

100 50 20 10
Period (years)

5

2

d)

1

LIM
Cane-Zebiak (Unforced)
Reconstructions

1

0

10

-1

0

10

-1

10

1000

LIM
Reconstructions
CMIP5 controls

10
Power

10
Power

0

10

-1

10

1000

LIM
Reconstructions
CCSM4 control
CCSM4 last 1ka

10
100 50 20 10
Period (years)

5

2

1000

100 50 20 10
Period (years)

52

Figure S4: Same as Figure 2 of the main text, but using the Blackman-Tukey method of estimating
the spectrum (with a smoothing window of 60 years).

13

a)

b)
LIM
CCSM4
Instrumental

1

Spectral Density (oC2/∆ f)

Spectral Density (oC2/∆ f)

1

10

0

10

-1

10

1000

100 50 20 10
Period (years)

5

2

c)

LIM
Reconstructions
CCSM4 control
CCSM4 last 1ka

0

10

-1

10

1000

100 50 20 10
Period (years)

5

2

d)
10

LIM
Cane-Zebiak (Unforced)
Reconstructions

1

Spectral Density (oC2/∆ f)

1

Spectral Density (oC2/∆ f)

10

0

10

-1

10

1000

100 50 20 10
Period (years)

5

2

10

LIM
Reconstructions
CMIP5 controls

0

10

-1

10

1000

100 50 20 10
Period (years)

52

Figure S5: Same as Figure 2 of the main text, but using a LIM developed from ERSSTv3.

14a) CCSM4 NINO3.4 spectra compared with instrumental LIM
1

Power

10

LIM
JEG131150-1850
CCSM4control

010

CCSM4850-2005
CCSM4850-1850
CCSM41000-1850-1

10

CCSM41150-1850

1000

100
50
20
Period (years)

10

5

2

b) CCSM4 NINO3.4 spectra compared with CCSM4-based LIM
1

Power

10

CCSM4-LIM
JEG131150-1850
CCSM4control

010

CCSM4850-2005
CCSM4850-1850
CCSM41000-1850

-1

10

1000

CCSM41150-1850

100
50
20
Period (years)

10

5

2

Figure S6: Power spectra of CCSM4’s NINO3.4 time series calculated over various time domains
and compared with: a) the LIM described in the main text of late 20th Century instrumental
data; and b) a LIM constructed from a 500 year segment of the CCSM4 control run. Note, thatthe centennial variability in last millennium run is above that produced by a LIM developed fromCCSM4’s own climate (i.e., a LIM constructed from the control run).

15a) Coral
2

Power

10

0

10

-2

10

-3

10

-2

-1

10

10

b) Ice
2

Power

10

Figure S7: Power spectra calculated fromthe individual proxy records used in EG13a,b.
Spectra are shown from the network used to reconstruct HadSST2i, although results are similar for the ERSSTv3 and Kaplan reconstructions. Gray lines in each panel show raw power
spectra from each type of proxy. Mean spectral
densities of each proxy type, calculated as theaverage of log-transformed spectra, are shownin black. Estimates were made using the LombScargle method because some records were notcontinuous through time. All time series werenormalized to exhibit unit variance from 1800
CE onward prior to calculating their spectra.

0

10

-2

10

-3

10

-2

-1

10

10

c) Sedi
2

Power

10

0

10

-2

10

-3

10

-2

-1

10

10

d) Tree
2

Power

10

0

10

-2

10

-3

10

-2

-1

10

10
-1

Frequency (yr )16βa) Overall power law
2.5
2
1.5
1
0.5
0
-0.5
-1

coral

ice

sedi

tree

coral

ice

sedi

tree

coral

ice

sedi

tree

β

b) Interannual power law
2.5
2
1.5
1
0.5
0
-0.5
-1

β

c) Dec-Cen power law
2.5
2
1.5
1
0.5
0
-0.5
-1

Figure S8: Power law coefficients (β) calculated from each proxy type. Power laws relate frequency(f ) to spectral density (S(f )) as S(f ) ∝ f −β : the higher the value, the “redder” spectrum. Herepower laws are shown for the following frequency intervals: a) the overall power spectrum (i.e.,from 1/2∆t to 1/N years, where N is record length and 1/2∆t is the Nyquist frequency); b) theinterannual frequency window (1/2∆t to 1/10 year); and c) dec-cen (1/10 to 1/N year) timescales.

17a) Overall power law

1.5
0.5

β

1
0
-0.5
-1
tree

b) Interannual power law

1.5
0.5

β

1
sedi

0
-0.5

ice

-1
coral

c) Dec-Cen power law

1.5
0.5

β

1
0
-0.5
-1

Figure S9: Power laws (β) of each record used in the EG13a,b reconstructions of HadSST2i. As
in Figure S8, power laws are shown for: a) the overall spectrum; as well as b) interannual, and c)
dec-cen timescales. The color and size indicate the magnitude of β, with warmer colors indicating
redder spectra. Symbol shapes, indicated by the legend on left side of the figure, are the same as
those used in Fiure S8, while symbol size is proportional to the absolute value of β (with small
markers indicating values close to zero).

18

JEG Recons
Proxy recordσa) Palmyra (coral)
1.5
0
-1.5

σ

b) Lake Challa (sedi)
1.5
0
-1.5

σ

c) Cariaco Basin (sedi)
1.5
0
-1.5

σ

d) NADA PDSI (tree)
1.5
0
-1.5

σ

e) Chilean Cordillera (tree)
1.5
0
-1.5

σ

f) Cerro Baraja,MX (tree)
1.5
0
-1.5

σ

g) Cerro Baraja,MX (tree)
1.5
0
-1.5
1200

1300

1400

1500

1600

1700

1800

1900

2000

Time (CE)

Figure S10: Low-pass filtered time series of the eight longest records in the EG13a,b network (black),
shown with the low-pass filtered reconstruction of HadSST2i (gray). A 3rd order Butterworth filter
designed to remove fluctuations shorter than 100 years was applied to each time series to performthe filtering.

19Temp. Anomaly (oC)a) Last millennium temperature anomalies predicted by "Hasselmann" model
2
1
0
-1
-2
900

1000

1100

1200

1300
1400
Time (CE)

1500

1600

1700

50

20

10

5

1800

b) CCSM4 last millennium spectra
0

10

-1

Power

10

-2

10

-3

10

-4

10

1000

CCSM4 forcing850-1850
CCSM4 NINO3.4850-1850
100

2Period (years)

Figure S11: a) Time series of temperature anomalies (TRAD ) simulated by the “Hasselmann” model
described by equation (1) (yellow curve), shown with CCSM4’s NINO3.4 (black). b) Power spectraof TRAD (yellow) and NINO3.4 from CCSM4 (black).

20