ABSTRACT

Title of thesis: NO CHANGE IN ENSO HYDROCLIMATE
VARIABILITY AFTER THE INDUSTRIAL
REVOLUTION AS RECORDED IN 𝛿18O OF
TECTONA GRANDIS L.F. FROM SOUTHEAST
SULAWESI, INDONESIA

Sandy Hardian Susanto Herho, Master of Science,
2023

Thesis directed by: Dr. Michael N. Evans
Department of Geology, University of Maryland
College Park

El Niño-Southern Oscillation (ENSO) is a quasi-periodic interannual oscillation of the

ocean-atmosphere system in the tropical Pacific which greatly influences global climate

variability. However, the long-term response to greenhouse gas forcing is still controversial.

In this study, we measured the oxygen isotopic composition of 𝛼-cellulose samples at intra-

annual resolution from independently crossdated teak cores (Tectona grandis L. f.) collected

at Muna, Indonesia (5.3ºS, 123ºE, elev. 10m). The site and observation has been previously

shown to provide an indirect measure of ENSO activity via local precipitation amount

variations associated with ENSO. We created an ensembled composite of the interannual

variability for the period 1680-2005 (316 years) using empirical high pass filtering and

random sampling of intra-annual resolution measurements. In processing this time series

composite, we used Singular Spectrum Analysis (SSA) to high pass filter the data for

the interannual variability associated with ENSO. The annually-resolved composite time

series of 𝛿18O that we constructed has a higher resolution than other studies that have been



conducted to reconstruct ENSO-hydroclimate activities in the western tropical Pacific region

over this period. Using this 𝛿18O composite, we compared the distribution of events in the

period before and after the industrial revolution using the two-sample Kolmogorov-Smirnov

(KS) test. We found no statistically significant change in the distribution of 𝛿18O anomalies.

The same statistical test was applied to the Niño 3.4 reconstruction from the Last Millennium

Reanalysis (LMR). The results of this study suggest that if there is indeed a forced response

of ENSO, it is as yet indetectable. This may be because the forcing is not yet large enough

or the forced response is small relative to the unforced variability. Additional factors that

might explain this result in the 𝛿18O composite include its observational and interpretational

uncertainty, and in the LMR reconstruction, the scarcity of tropical observational constraints

and systematic error in the representation of ENSO in climate simulations.



NO CHANGE IN ENSO HYDROCLIMATE VARIABILITY AFTER THE
INDUSTRIAL REVOLUTION AS RECORDED IN 𝛿18O OF TECTONA

GRANDIS L.F. FROM SOUTHEAST SULAWESI, INDONESIA

by

Sandy Hardian Susanto Herho

Thesis submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Master of Science

2023

Advisory Committee:
Dr. Michael N. Evans, Chair/Advisor
Dr. Karen L. Prestegaard
Dr. Raden D. Susanto



© Copyright by
Sandy Hardian Susanto Herho

2023



Acknowledgments

I am grateful to Austin Carter, Austin Mickels, Cristy Ho, Jenna Wollney, and Mary

Gbenro for laboratory support. I also benefitted from helpful discussions with Faiz Fajary,

Karen Prestegaard, Michael N. Evans, and R. Dwi Susanto. I also thank Ann Wylie (through

the Green Fellowship in Global Climate Change) and Neera Gupta, who partially financed

my study at UMD. This study was funded by the NSF grant AGS1903626.

ii



Table of Contents

Acknowledgements ii

1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Is ENSO changing? Evidence from climate models . . . . . . . . . . . . . 1
1.3 Is ENSO changing? Evidence from paleoclimatology . . . . . . . . . . . . 3
1.4 This study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.4.1 Study location . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4.2 Data modeling: 𝛿18O of 𝛼−cellulose . . . . . . . . . . . . . . . . . 6
1.4.3 Goals of this study . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Materials and methods 9
2.1 Stable isotope analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1.1 Imaging, intra-annual sampling, and 𝛼−cellulose extraction . . . . . 9
2.1.2 𝛿18𝑂𝑐 data acquisition . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Time-series compositing . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2 Nonparametric statistical test . . . . . . . . . . . . . . . . . . . . . 13

3 Results 15

4 Discussion 22
4.1 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.2 Comparison with prior studies . . . . . . . . . . . . . . . . . . . . . . . . 22

5 Conclusion 24

A Statistical analyses for other embedding dimensions of the high-pass filtered time
series 26

Bibliography 33

iii



Chapter 1: Introduction

1.1 Motivation

The El Niño/Southern Oscillation (ENSO) phenomenon produces the most influential

impacts on the climate of the tropical Pacific, where it occurs, as well as having far-

reaching global consequences [11, 22, 60]. One of the areas that most affected by ENSO

is the western equatorial Pacific, since ENSO is a major contributor to the interannual

and interseasonal atmospheric convergence and precipitation in this region [50]. While

the interannual variability of ENSO under normal conditions (without external radiative

forcings) has been extensively studied, it is currently uncertain whether the patterns of its

occurrence are significantly affected by external radiative forcing that arise from human

activities, such as greenhouses gas (GHG) emissions [16, 17].

1.2 Is ENSO changing? Evidence from climate models

Several studies using climate models indicate an increase in ENSO cold and warm

phase event amplitude, suspected to be related to anthropogenic forcing since the industrial

revolution [39, 44, 48]. This occurs because of a balance between the release of latent

heat in precipitation and radiative cooling in the troposphere [2]. This fact results in an

increase in ENSO-driven precipitation in the western equatorial Pacific as a consequence of

an increase in GHG emissions. However, the forcing of cold phase events by GHG increases

1



is small relative to what has been estimated from natural forcings [44]. This ENSO-driven

precipitation intensification process in the western equatorial Pacific can be explained by the

ocean-dynamical thermostat theory [5, 10, 14, 44]. Increased warming at the surface will

have more impact on the western equatorial Pacific than on the eastern equatorial Pacific, this

is because warming in the eastern region is offset by the upwelling process. This increase

in the SST gradient then exerts an additional force on the pressure gradient which causes

the Walker circulation to become stronger, and ultimately increases the SST gradient even

more. This mechanism is known as Bjerknes feedback [44]. This finding is still debatable

because the consensus results from the simulations of the Coupled Model Intercomparison

Project-Phase 3 and 5 (CMIP3 and CMIP5) that the Walker circulation weakens, which

will reduce trade winds and make the equatorial thermocline lower, resulting in more ”El

Niño-like” phenomena [62, 63, 69]. Another study conducted using the Community Earth

System Model-Last Millennium Ensemble (CESM-LME) project by Stevenson et al. [59]

revealed that the dominant ENSO variability in the last millennium was governed by internal

variability. Although there appears to be an increase in ENSO variability in the 20th century,

this signal is not robust. However, using the Last Millennium experiments of CMIP5 and

Paleoclimate Modelling Intercomparison Project-Phase 3 (PMIP3), Lewis and LeGrande

[42] and Hope et al. [37] instead found an increase in ENSO variability in the historical

period (1906 - 2005) as a result of an increase in GHG concentrations and land-use changes

after the industrial revolution. It can be concluded that there is still debate about how ENSO

changes in response to anthropogenic forcings as projected by various climate models after

the industrial revolution [17, 32, 41, 72].

2



1.3 Is ENSO changing? Evidence from paleoclimatology

Additional studies have analyzed paleoclimatic observations and reconstructions of

ENSO indicators for changes in properties before and after the anthropogenic increase in

radiative forcing of climate, and to establish the range of natural, or unforced variability[e.

g. 15, 43, 45, 46, 53, 56]. For example, many study area selections include broad networks

of observations [e. g. 43], but are based solely on statistical correlation [e. g. 45] rather

than on the forward mechanistic modeling of proxies [e. g. 24]. Others [e. g. 20] are from

centers of action for ENSO, but have relatively low reconstruction skill. Leveraging forward

modeling of proxies with respect to their environmental controls [26], paleoclimatic data

assimilation can be used to update climate model simulations with proxy observations [e.

g. 33]. In the approach of Hakim et al. [33] and related studies [73, 74],

x𝑎 = x𝑝 +K(y−H(x𝑝)), (1.1)

, where x𝑎 is the posterior estimate of Niño 3.4 reconstruction, x𝑝 is the prior estimation from

climate simulations, K is a Kalman gain matrix, y is the paleoclimatic observations, and H

is the forward operator used to map x𝑝 to the proxy domain. Note that in this formulation,

variance change over time in these products may be convolved with changes associated

with observational availability y over time [e. g. 66]. In the work of Zhu et al. [73, 74],

CCSM4-LME simulations [40] were used as the prior, and only coral-based observations

and tree-ring records previously shown to record ENSO activity [19, 21, 43] were used

for observations. In this study, Zhu et al. [73] found no evidence of an ENSO-mediated

response to radiative forcing associated with explosive volcanic events. Historical analysis

of gridded SST data, based on dense observations relative to paleoclimatic datasets, together

with ENSO modeling suggest no forced response for the period 1856-2003 [12]. Overall,

3



there does not appear to be consensus from these studies on the anthropogenically forced

response to ENSO since the Industrial Revolution. Thus, for the paleoclimatic strategy, it

would be beneficial to acquire paleoclimatic data, whose interpretation is underpinned by

mechanistic understanding, and which is from key locations in which there is an ENSO-

related environmental response, to improve assessments of the sensitivity of ENSO to

post-industrial anthropogenic forcing.

1.4 This study

1.4.1 Study location

We addressed this problem by using the annually resolved oxygen isotopic composi-

tion of 𝛼−cellulose (𝛿18𝑂𝑐) from teak (Tectona grandis L. f.) increment tree cores from

Muna, Southeast Sulawesi, Indonesia. These samples were previously and independently

dendrochronologically crossdated by D’Arrigo et al. [20] spanning the period 1656-2005

with a sample depth for each year of at least six increment cores per year [49]. The oxygen

isotopic composition of cellulose in tree rings in Southeast Asia has been successfully used

for reconstructing precipitation and ENSO-driven hydroclimatic responses [51, 58, 68, 75].

Muna is located in in the western equatorial Pacific, known as the Maritime Continent

(MC) [52](20◦S - 20◦N, 90◦E - 160◦E, Figure 1.1). This area is directly affected by ENSO

from a hydroclimatological perspective. Many studies have been conducted that specifically

discuss the close relationship between the variability of precipitation over the MC and the

ENSO [29, 34, 47, 54, 70]. The negative precipitation anomaly over the MC is considered

to have a strong relationship with the ENSO warm phase. The precipitation deficit over the

MC is generally associated with a large-scale air mass subsidence anomaly, a weakening

of the Walker circulation, and a negative anomaly of sea surface temperature (SST) in the

oceans surrounding the MC [54]. Muna is located on the border of region A and region C

4



in the study of monthly precipitation zone delineation using the double correlation method

(DCM) by Aldrian and Susanto [1]. Region A has one peak, as a result of the wet northwest

monsoon from November to April, and one trough, which occurs as a consequence of the

dry southeast monsoon in the period March to May, in its average annual precipitation

pattern. Meanwhile, region C has peak rainfall from June to July and one trough from

November to February. Based on the precipitation spectral analysis conducted by Aldrian

and Susanto [1], it is known that ENSO has a strong interannual effect on the two regions

which are considered to represent the Muna region.

Figure 1.1: Map showing the study location (red dots) in a tropical lowland rain forest (10 m.a.s.l)
in the southeastern part of Sulawesi (5.3◦S, 123◦E).

5



1.4.2 Data modeling: 𝛿18O of 𝛼−cellulose

Nurhati et al. [49] used a forward model of 𝛿18𝑂𝑐 and comparison of simulations with

triplicate observations for the period 1971-2005 to show that ENSO-driven precipitation

acts as the main control of 𝛿18𝑂𝑐 in this location. They used the 𝛿18𝑂𝑐 model of Barbour

et al. [4], as modified by Evans [24] for tropical regions,

𝛿18𝑂𝑐 = 𝑓
(
𝑇, 𝑝𝑟, 𝑅𝐻

)
(1.2)

, where

• 𝑇 : monthly average air temperature (◦C),

• 𝑝𝑟: monthly precipitation (mm),

• 𝑅𝐻: relative humidity (%),

• 𝑓 : forward-deterministic model function as described in Evans [24].

Here, 𝑓 parameterizes the isotopic composition of precipitation (𝛿18𝑂𝑝) as a linear function

of precipitation amount, and of atmospheric water vapor (𝛿18𝑂𝑣) as a constant value of 8‰

lower than 𝛿18𝑂𝑝.

𝛿18𝑂𝑐 can be modeled deterministically into the following equations,

Δ𝑐 = Δ𝑙
(
1− 𝑝𝑥 𝑝𝑒𝑥

)
+ Y𝑐 (1.3)

Δ𝑙 =
{
(1+ Y∗) [1+ Y𝑘 + (Δ𝑤𝑣𝑎 − Y^

𝑒𝑠

𝑒𝑖
) −1]

} (
1− 𝑒𝑃𝑒
𝑃𝑒

)
. (1.4)

Here, Δ is defined as the oxygen isotopic composition relative to the oxygen isotopic

composition of the source water absorbed by plants, which is assumed to be equal to

6



unmodified 𝛿18𝑂𝑝. The
(
1− 𝑝𝑥 𝑝𝑒𝑥

)
fraction of the leaf water isotopic fractionation relative

to source water (Δ𝑙) is an isotope modification caused by the evaporation from leaves

and by a partial re-equilibration process between the leaves and unmodified stem water

(equation 1.3). The leaf evaporative enrichment (Δ𝑙) is calculated as a function of the in-

situ modified oxygen isotopic composition, temperature-dependent equilibrium and kinetic

fractionation factors, and water vapor isotopic fractionation driven by leaf-air differences

in specific humidity, modified by the Péclet effect (equation 1.4). Within this model, 𝛿18𝑂𝑐

is predicted from the oxygen isotopic composition of precipitation, calculated evaporative

enrichment at the leaf, partial re-equilibration of leaf water by Péclet back-diffusion, partial

re-equilibration of photosynthate with unmodified stem water, and fractionation at cellulose

synthesis.

Evans [24] adapted the model from Barbour et al. [4] by parameterizing leaf temperature

from air temperature, 𝛿18𝑂𝑝 as a linear function of precipitation rate, cloud temperature

as a linear function of moist adiabatic lapse rate and condensation level, and leaf stomatal

conductance as an inverse function of leaf-air saturation vapor pressure deficit, calculable

from temperature and relative humidity [24]. Thus, this model can be run with the inputs

in equation (1.2).

The temporal correlation between triplicate-observed Muna 𝛿18𝑂𝑐, and November-

April average 𝛿18𝑂𝑐 simulation, produced from Climate Research Unit (CRU) TS3.01

temperature, precipitation and humidity [35] data at Muna, is statistically significant (1971-

2005: 𝑟 = 0.36, 𝑒𝑑𝑓 = 35, 𝑝 = 0.03) [49]. As predicted by the 𝛿18𝑂𝑝 anomalies produced for

ENSO events in an isotopically enabled reanalysis of atmospheric climate [71], the observed

Muna 𝛿18𝑂𝑐 record is significantly correlated with the Niño 3.4 index (SST anomalies

(averaged over 5◦N - 5◦S, 170◦W - 120◦W; 1971-2005, Sep-Oct annual averages: 𝑟 = 0.52,

𝑑𝑓 = 35, 𝑝 < 0.01), and the statistically significant pattern correlation of reanalysis 𝛿18𝑂𝑝

[71] is positive over the MC and negative over the central equatorial Pacific [49], consistent

7



with climatological understanding of ENSO as described earlier.

1.4.3 Goals of this study

Here we assume that the same mechanistically developed interpretation applies to inter-

pretation of the interannual component of Muna composite 𝛿18𝑂𝑐 developed for the period

extending into the pre-industrial past. By comparative analysis of data for the periods before

and after the exponential increase in CO2 emissions starting in the mid-19th century, we can

assess whether there was an associated change in ENSO activity. Specifically, we can test

the prevailing hypothesis [7–9] that the probability of ENSO-hydroclimate extremes due to

global warming has increased in the industrial period [7–9]. To test the sensitivity of the

results and interpretation to uncertainty in the observation and interpretation of the Muna

𝛿18𝑂𝑐 record, we can perform the same analysis but using as input the independently devel-

oped Niño 3.4 index from the Last Millennium Reanalysis ”turbo” (LMRt) paleoclimatic

data assimilation [73, 74], for which 1941-2000 Niño 3.4 validation 𝑟2 = 0.62).

8



Chapter 2: Materials and methods

2.1 Stable isotope analysis

In this study, we sampled teak (Tectona grandis L. f.) increment cores in Muna, Southeast

Sulawesi, Indonesia (5.3◦S, 123◦E) (Figure 1.1) which was dendrochronologically cross-

dated by D’Arrigo et al. [20]. We analyzed three crossdated samples, tg01gc (1712 -

1995), tg11a (1680 - 1848), and mun6.3 (1798 - 1987). Procedures for preparing samples

for isotopic analysis included sample imaging, microtoming, extracting 𝛼−cellulose, and

weighing and wrapping samples in silver capsules. To acquire 𝛿18𝑂𝑐 data, we thermally

converted samples to analyte gas in an elemental analyzer (EA), introduced the analyte into

a stable isotope ratio mass spectrometer (IRMS) via a continuous flow interface, collected

carbon and oxygen data, and calibrated the data using a 2-point correction for mean and

variance to international reference scales [25]. All laboratory work was performed at the

Paleoclimate CoLaboratory of the University of Maryland, College Park.

2.1.1 Imaging, intra-annual sampling, and 𝛼−cellulose extraction

Sample cores were scanned at 1200 dpi to match the core names, orientations, and scales

before being separated from their wooden bases. We used a combination of rotary microtome

and manual sampling using a razor blade to subsample each crossdated growing season into

at least six sequential samples per year. Samples were placed into 1.5 ml centrifuge tubes

for 𝛼−cellulose extraction following a modified Brendel method; a detailed description is

9



in Anchukaitis et al. [3]. This procedure was done to isolate 𝛼−cellulose, for which we

have a well-developed and tested mechanistic interpretative model (see subsection 1.4.2)

[24, 27] from other wood components. This extraction process was carried out by adding

nitric acid and acetic acid in a ratio of 1:10 to the sample tubes, then capped sample tubes

were heated to 119 - 124◦C for 30 minutes. Next, samples were sequentially rinsed with

ethanol, then water, then ethanol, then acetone, before being dried at 45◦C for 30 minutes.

Extracted 𝛼−cellulose was then stored in a vacuum dessicator for 24 hours or until isotopic

analysis. Then each of these dried samples was wrapped and weighted to 200±20 `g or

100±10 `g (depending on the sample size) in a silver capsule to be analyzed for isotopic

composition.

2.1.2 𝛿18𝑂𝑐 data acquisition

For the analysis of 𝛿18𝑂𝑐, both samples and standards were run so data could be

normalized, allowing drift in the instrument to be accounted for and comparisons to be

made with other datasets. Samples and standards were loaded into the autosampler in order,

and names and weights were filled into the batch processing file, with care taken to ensure

the information matched the order of the samples and standards in the autosampler. The

autosampler was then purged with helium to remove any remaining traces of air. Samples

were left under helium overnight to make sure they were dry before running the analysis.

Once each run was started, C and O from the samples were converted into CO via thermal

conversion at 1080◦C, and CO2 and H2O were removed with an acid/water trap before

sample CO was admitted to the IRMS by continuous flow interface. Currents measured for

masses 28 (12C16O), 29 (13C16O), and 30 (12C18O) were used to determine isotope ratios.

Two independent working standards, SAC (𝛿13𝐶 = −25.13± 0.19 ‰; 𝛿18𝑂 = 31.0± 0.24

‰) and AKC (𝛿13𝐶 = −13.80±0.10 ‰; 𝛿18𝑂 = 23.61±0.24 ‰) were used for correction

10



of systematic mean and variance biases introduced by drift and amplitude through each

sample batch using scripted functions in MATLAB® [25]. Diagnostic plots were used to

check the correction algorithm produced valid corrections. Uncertainty in the isotopic data

was estimated as the standard deviation of replicate corrected working standard data across

all working standards, and is estimated as 0.4 ‰ averaged across 36 batches of samples and

working standards.

2.2 Statistical analysis

2.2.1 Time-series compositing

A compositing step is required to produce an annual unidimensional time series of

𝛿18𝑂𝑐 based on the three cores of differing time intervals and intra-annual resolutions. The

first step was random sampling with replacement of a sub-annual 𝛿18𝑂𝑐 value from each

available annual increment. This step produced 1000 realizations of each of the three time

series. We then subtracted each series mean, producing anomaly series, and to high-pass

filter for analysis of ENSO-band variations, we used singular spectrum analysis (SSA).

SSA is a nonparametric time-series analysis method which can be implemented on

series with missing data [28, 30, 36, 57, 61]. SSA can be divided into two main steps: the

decomposition, which consists of embedding and singular value decomposition (SVD), and

the reconstruction, which consists of grouping and diagonal averaging [30, 36].

The embedding was implemented for each anomaly realization, a unidimensional series

Y𝑁 = {𝑦1, · · · , 𝑦𝑁 } to be mapped into a multi-dimensional series X1, · · · ,X𝐾 , with vector

x𝑖 = (𝑦𝑖, · · · , 𝑦𝑖+𝐿−1)𝑇 ∈ R𝐿 , where 𝐾 = 𝑁 − 𝐿 + 1. The window length 𝐿 has to satisfy the

requirement: 2 ≤ 𝐿 ≤ 𝑁/2. The result of this embedding calculation is a trajectory matrix

X with size 𝐿×𝐾 which is a Hankel matrix, where all anti-diagonal elements have the same

11



value, expressed in the following form,

X =
(
x𝑖 𝑗

)𝐿,𝐾 =

©­­­­­­­­­«

𝑦1 𝑦2 . . . 𝑦𝐾

𝑦2 𝑦3 . . . 𝑦𝐾+1
...

...
. . .

...

𝑦𝐿 𝑦𝐿+1 . . . 𝑦𝑁

ª®®®®®®®®®¬
. (2.1)

After the embedding calculation, the decomposition step was continued by applying

SVD to the trajectory matrix X. We began SVD analysis by defining the covariance matrix

S = XX𝑇 . Then using the following equation,

det(S−_I) (2.2)

, we determined the eigenvectors U𝑖 and their corresponding eigenvalues _𝑖. This procedure

was followed by constructing a singular value matrix 𝚺,

𝚺 =

©­­­­­­«

√
_1 . . . 0
...

. . .
...

0 . . .
√
_𝐿

ª®®®®®®¬
. (2.3)

The principal component matrix is the transpose of the V𝑖 matrix, which was calculated

through the following procedure,

V𝑇
𝑖 =

(
X𝑇U𝑖√
_𝑖

)𝑇
. (2.4)

The final result of the SVD is the decomposition of the trajectory matrix X into eigentriple

12



as follows,

X = U𝑖𝚺V𝑇
𝑖 . (2.5)

The matrix X𝑖 = X1 +X2 + · · · +X𝑑 from the eigentriple was then used to reconstruct the

time series through a grouping process by grouping the set of indices 𝑖 = {1,2, · · · , 𝑑} into

a subset of disjoint 𝐼 = {𝐼1, 𝐼2, · · · 𝐼𝑚}, with 𝑚 = 𝑑, so that it could be expanded into this

following form,

X𝐼 = X𝐼1 +X𝐼2 + · · ·X𝐼𝑚 . (2.6)

The final step to reconstruct the time series based on the power of the spectrum is to perform

diagonal averaging. This stage was performed by transforming each matrix in the equation

2.6 into a Hankel matrix to produce the initial time series form as follows,

y𝑡 =
𝑚∑︁
𝑘=1

�̃�𝑚𝑡 , , where 𝑡 = 1, · · ·𝑚. (2.7)

Because we aimed to do high-pass filtering on each anomaly realization to get interannual

variability, we substituted𝑚 = 7 years and removed �̃�1. The ensemble of 3×1000 high-pass

filtered 𝛿18𝑂𝑐 anomaly timeseries, spanning the entire observational period 1712-1995,

were then used to estimate the composite anomaly median and 95th percentile confidence

intervals. We implemented all of the time-series compositing processes in the MATLAB®

computing environment.

2.2.2 Nonparametric statistical test

The Kolmogorov-Smirnov (KS) test is a nonparametric test of continuous, one-dimensional

probability distributions that can be used to compare a sample with a reference probability

distribution (one-sample KS test), or to compare two samples (two-sample KS test) [65, 67].

In this study, we used the two-sample KS test. A nonparametric statistical test was used

13



because it does not assume normality in the time series distribution. The KS test was used to

test the null hypothesis that two independent samples come from populations with identical

distributions.

• 𝐻0: 𝐹 (𝑥) = 𝐺 (𝑥), for all 𝑥.

• 𝐻1: 𝐹 (𝑥) ≠ 𝐺 (𝑥), for at least one value of 𝑥.

, where 𝐹 (𝑥) and 𝐺 (𝑥) are the cumulative distribution functions (CDFs) of highpassed

𝛿18O𝑐 anomalies (hereafter annual 𝛿18𝑂𝑐 variability) for the population before and after

the industrial revolution, respectively. Let 𝑆𝑛1(𝑋) and 𝑆𝑛2(𝑋) be the empirical cumulative

distribution functions of the 𝛿18O𝑐 anomalies before and after the baseline periods, respec-

tively, given by 𝑆𝑛1 = 𝑘/𝑛1 and 𝑆𝑛2 = 𝑘/𝑛2, where 𝑘 is the same or less than 𝑋 , and 𝑛1 and

𝑛2 are the number of data points in each observation. The KS statistic is

𝐷 = max |𝑆𝑛1(𝑋) − 𝑆𝑛2(𝑋) |. (2.8)

To test the sensitivity of the test to definition of the pre-industrial and industrial eras, we

used two defensible definitions: (1) 1850 and (2) 1970, defined by the initial point at which

anthropogenic greenhouse forcing is estimated to have begun, and the point at which the

decadally averaged rate of increase in greenhouse radiative forcing suddenly increased,

respectively [72]. We also tested the sensitivity of the results to sample size imbalance

by repeating the latter test for independent 25-year time windows before and after 1970.

We used the scipy library [64] in the Python computing environment for test calculations,

and consider p𝑐𝑟𝑖𝑡 = 0.05 for inference. We applied the same test and definitions to the

independent December-February (DJF) Niño-3.4 index from the LMRt paleoclimate data

assimilation product.1

1Data available at: https://zenodo.org/record/5893781.ZAqnh3WYXIE.

14

https://zenodo.org/record/5893781##.ZAqnh3WYXIE


Chapter 3: Results

There are 1,356 sub-annual measurements of 𝛿18𝑂𝑐 on tg01c, 788 on tg11a, and 1,048

on mun6.3. Intra-annual resolution 𝛿18𝑂𝑐 measurements are shown as time series in Figure

3.1. From tg01c, the maximum 𝛿18𝑂𝑐 is 30.7 ‰ in 1928, and the minimum is 19.3 ‰

in 1765 (Figure 3.1a). From tg11a, the maximum value of 𝛿18𝑂𝑐 is 28.7 ‰ in 1686, and

the minimum value is 22.2 ‰ in 1720 (Figure 3.1b). From mun6.3, the maximum 𝛿18𝑂𝑐

is 32.2 ‰ in 1981, and the minimum value is 23.0 ‰ in 1934 (Figure 3.1c). The median

and 2.5th - 97.5th percentile for each core are shown in Table 3.1 (𝛿18𝑂𝑐 values). Series

medians differ by almost 2 permil (Table 3.1), and the 95th percentile ranges of values are

about 4 ‰, 4 ‰, and 5 ‰, respectively (Table 3.1, 𝛿18𝑂𝑐 values).

core name 𝛿18O𝑐 𝛿18O𝑐 variability

median (‰) [2.5th, 97.5th] (‰) median (‰) [2.5th, 97.5th] (‰)

tg01c 25.53 [23.56, 27.52] 0.01 [-1.43, 1.4]
tg11a 25.32 [23.41, 27.27] 0.03 [-1.54, 1.34]

mun6.3 27.33 [24.95, 29.75] -0.01 [-1.94, 2.22]

Table 3.1: Median and 2.5th - 97.5th percentile confindence limits for each core. Columns 2 - 3
are statistics for the intra-annual 𝛿18O𝑐 series. Columns 4 - 5 are statistics for 𝛿18O𝑐 variability by
sample core prior to compositing across cores.

The high pass 𝛿18𝑂𝑐 composite is shown in Figure 3.2. The solid green line indicates

the median values. Meanwhile, the green shaded area shows the 95% percentile interval

obtained from 3×1000 randomized realizations. The highest median of the 𝛿18𝑂𝑐 variability

is 1.40 ‰ in 1712, and the lowest is -1.90 ‰ in 1706. The grand median of 𝛿18𝑂𝑐 variability

15



(a)

(b)

(c)

Figure 3.1: Intra-annual resolution 𝛿18𝑂𝑐 measurements collected from Muna increment core
samples: (a) tg01c, (b) tg11a, and (c) mun6.3. Vertical axis shows oxygen isotopic composition
in permil units relative to the Vienna Standard Mean Ocean Water standard (VSMOW), with each
scale range = 10 ‰ centered on the median series value. The single point in black at 1700, 27.33
shows observational uncertainty for intra-annual measurement analytical uncertainty, as ±2 standard
deviations of the mean of corrected standard measurements. Sequential sample data were assigned
sequential dates for Nov of the prior calendar year through April of the following calendar year by
linear interpolation; horizontal axis shows time (years CE).

16



is 0.00 ‰, and the ensemble 95th percentile confidence interval is [-0.98 ‰, 0.88 ‰]. The

ensemble 𝛿18𝑂𝑐 variability is our best estimate of the site composite 𝛿18𝑂𝑐 record at

interannual timescales, sampling observational uncertainty and filtering for mechanistically

supported timescales of variation.

Figure 3.2: 𝛿18𝑂𝑐 variability composite at annual resolution. Solid green line shows the median,
while shaded region denotes the central 95% quantiles from the 3×1000 random realizations. Error
bar at point (0,1670) shows ±2 SE of the observational uncertainty for the annual 𝛿18𝑂𝑐 variability
composite, estimated as 𝑠/

√
𝑛, 𝑛 = 6 samples/year, and neglecting the reduction of variance associated

with highpass filtering.

The results of empirical CDF calculations of the 𝛿18𝑂𝑐 variability from the pre- and

post-industrial revolution periods are shown in Figure 3.3. Figure 3.3a shows the CDFs

of the 𝛿18𝑂𝑐 variability for the 1680 - 1849 (as the pre-industrial period) and the 1850

- 1995 (as the post-industrial period). Meanwhile, Figure 3.3b shows the CDFs of the

𝛿18𝑂𝑐 variability for the 1695 - 1969 (as the pre-industrial period) and the 1970 - 1995

(as the post-industrial period). The results of the two-sample KS test showed that there

17



is no significant difference in the CDFs between the pre-and post-industrial periods, both

when the industrial revolution began in 1850 (with statistical values: 𝐷 = 0.072, 𝑝−value

= 0.780), or when the industrial revolution began in 1970 (Table 3.2). This result was not

sensitive to the choice of SSA embedding dimension for 𝑚 = 5,7,9,11 (Appendix A).

Table 3.2: Statistical results of the two-sample KS test of 𝛿18𝑂𝑐 anomaly for every 25 years compared
to the period after the industrial revolution (1970 - 1995).

period 𝐷𝑚𝑎𝑥 𝑝−value
1695 - 1719 0.177 0.733
1720 - 1744 0.142 0.910
1745 - 1769 0.217 0.495
1770 - 1794 0.268 0.233
1795 - 1819 0.226 0.424
1820 - 1844 0.188 0.652
1845 - 1869 0.148 0.882
1870 - 1894 0.137 0.929
1895 - 1919 0.134 0.942
1920 - 1944 0.266 0.241
1945 - 1969 0.257 0.297

As a comparison, we also applied the two-sample KS test to the LMRt DJF Niño

3.4 reconstruction [73, 74] during the same period (1680 - 1995) (Figure 3.4a). Figure

3.4b shows the CDFs in 1680 - 1849 (as the pre-industrial period) and 1850-1995 (as the

post-industrial period). Figure 3.4c shows the CDFs in the period 1680 - 1969 (as the

pre-industrial period) and 1970-1995 (as the post-industrial period). Based on the a-priori

𝑝 = 0.05 critical level, we infer that there is no significant difference between the CDFs for

the periods 1680 - 1849 and 1850 - 1995, with statistical values: 𝐷 = 0.140 and 𝑝−value

= 0.082 (Figure 3.4b). We make the same inference if the industrial revolution period is

defined to start in 1970 (Table 3.3, Figure 3.4c).

18



Table 3.3: Statistical results of the two-sample KS test of Niño 3.4 reconstruction for every 25 years
compared to the period after the industrial revolution (1970 - 1995).

period 𝐷𝑚𝑎𝑥 𝑝−value
1695 - 1719 0.345 0.060
1720 - 1744 0.346 0.059
1745 - 1769 0.228 0.413
1770 - 1794 0.268 0.233
1795 - 1819 0.305 0.128
1820 - 1844 0.268 0.233
1845 - 1869 0.305 0.128
1870 - 1894 0.300 0.147
1895 - 1919 0.285 0.216
1920 - 1944 0.331 0.098
1945 - 1969 0.337 0.079

19



(a)

(b)

Figure 3.3: Empirical cumulative probability density function (CDF) curves of the annual 𝛿18𝑂𝑐

variability. The red line shows the period after the industrial revolution, while the blue lines show
the periods before the industrial revolution. We used two different definitions of the anthropogenic
era to show that the results are not sensitive to its definition of the beginning of the exponential
increase in the anthropogenic forcings. In (a) the industrial revolution began in 1850, while in (b)
the industrial revolution began in 1970.

20



(a)

(b)

(c)

Figure 3.4: (a) DJF Niño 3.4 reconstruction (1680 - 1995) using the LMR turbo (LMRt) framework
[74] from the analysis of Zhu et al. [73], solid red line shows median values; shaded area denotes
the 95% quantiles. (b - c) same as the Figure 3.3, but for the LMR data.

21



Chapter 4: Discussion

4.1 Summary of results

We constructed an annually resolved site composite record from three temporally over-

lapping, intra-annually resolved 𝛿18𝑂𝑐 series. Based on prior work demonstrating a link to

ENSO activity on interannual timescales, we infer no change in the inferred distribution of

ENSO events through the 316 year (1680-1995) composite record. This result is sensitive

to neither the parameters in data compositing and highpass filtering, nor the choice of

startpoint of significant greenhouse gas radiative forcing as either 1850 or 1970. Our study

shows no difference in ENSO variability as a result of changes in anthropogenic forcings

regardless of the initial period of the exponential increase.

4.2 Comparison with prior studies

Our results are independently replicated by reaching the same inference from analysis

of LMRt Niño-3.4 index reconstructions for the same intervals. However, we also find that

in these results, some intervals give significances that are only slightly above the canonical

statistically significant 𝑝−value (𝑝 < 0.05) (Table 3.3, Figure 3.4b, c). This result, which

is most apparent for the 17th and 20th centuries (Table 3.3), may be due to artifacts in

the LMR product. Historical observations [12]1 do not support the ≈1𝑜C warming in the

1https://iridl.ldeo.columbia.edu/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/

figviewer.html?plottype=line.

22

https://iridl.ldeo.columbia.edu/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/figviewer.html?plottype=line
https://iridl.ldeo.columbia.edu/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/figviewer.html?plottype=line


reconstructed Niño 3.4 index. One factor that might explain the 17th century change in

distribution could be the tendency of data assimilation products to resolve less variance as

observational coverage decreases. For the spurious 20th century warming and associated

change in distribution of Niño3.4 anomalies, we speculate that the forward operators used in

the data assimilation may have transformed a calibration period trend into a reconstructed

warming. Alternatively, this trend could arise from biased representation of ENSO, or

the response of ENSO to radiative forcing, in the CCSM4-LME simulations. Further

investigation is needed, both statistically (e.g. by using the False Discovery Rate (FDR)

[18]) and dynamically (e.g. by correcting the cold tongue bias representation of the input

model used in the LMRt framework [31]), to address this issue. Our findings do differ

from the those of Liu et al. [45]. In their analysis of 𝛿18𝑂𝑐 observations from Taiwan for

1190-2007 CE, they inferred an increase in ENSO variability at the end of the 20th century.

However, Liu et al. [45] did not test statistical significance of variance changes in their

study.

23



Chapter 5: Conclusion

Our inference is consistent with most studies of paleo ENSO since the Pliocene, which

suggests that the external radiatively forced response of ENSO may be small relative to

the unforced internal variability [23]. This conclusion is provisional because there are

considerable uncertainties in the observations, reconstructions and modeling of ENSO and

its forced response. Nurhati et al. [49] showed that the interpretation of the Muna composite

𝛿18𝑂𝑐 variability record explains less than 1/3 of the interannual variance associated with

ENSO, and the Muna record is from only one location. Further confidence can be developed

from paleoclimatic data assimilation, but products such as LMRt also have uncertainties

associated with data modeling, changes in observational coverage over time, and biases in

the prior climate states used. Finally, there is as yet no model consensus on the forced

response of ENSO. At present, the evolution of ENSO amplitude, frequency, and variability

of anthropogenic forcings is still poorly constrained, and there is no inter-model consensus

on projected changes in the dynamics of ENSO diversity in response to the exponential

increase in GHG since the industrial revolution [6, 13, 17]. A more in-depth physical

investigation is needed on the relationship of anthropogenic forcings to changes in ENSO

variability, for instance concerning ENSO diversity, represented by Central Pacific (CP)

and Eastern Pacific (EP) events [38].

Historical observations suggest a high degree of inter-event variations from the mid-

19th century to the present, consistent with what an intermediate compexity ENSO forecast

model might produce without anthropogenic forcings [12]. Paleoclimatic observations

24



and reconstructions extend this result further into the pre-anthropogenic past, but their

analysis fails to disprove the hypothesis that the forced signal, if present, is difficult to detect

within the large spread of the natural variability. Such efforts will be improved by the

establishment of more high resolution observations in other ENSO-sensitive regions, such

as equatorial South America [55], data modeling, and the further development of spatially

resolved data assimilation products. Until either the climate simulations produce a clearer

consensus, and actual radiative forcing and observed variability exceeds the historical and

paleoclimatogically observed or reconstructed envelope, the effect of anthropogenic forcings

on ENSO variability will remain unclear.

25



Appendix A: Statistical analyses for other embedding dimensions of the

high-pass filtered time series

To test that the computational results of our high-pass filtered time series are not sensitive

to the arbitrary embedding dimension (𝑚 = 7), we performed high-pass SSA filtering with

𝑚 = 5, 9, and 11. The results are shown in Figure A.1. After extracting 𝛿18𝑂𝑐 variabilities

with different embedding dimensions, we applied the two-sample KS test as we did in

subsection 2.2.2. Empirical CDFs for 𝑚 = 5, 9, and 11, respectively shown in Figures

A.2, A.3, and A.4. KS test results comparing CDFs (for 𝑚 = 5, 9, and 11) in 1680 - 1850

with 1850 - 1995 are shown in Table A.1. The results of the KS tests for the post-1970

period for 𝑚 = 5, 9, and 11 are shown in Tables A.2, A.3, and A.4, respectively. The KS

tests with these different embedding dimensions show no statistically significant differences

(𝑝 < 0.05) between 𝛿18𝑂𝑐 variabilities before and after the industrial revolution.

Table A.1: Statistical results of the two-sample KS test of 𝛿18𝑂𝑐 variabilities (𝑚 = 5,9,11) for 1680
- 1849 compared to the period after the industrial revolution (1850 - 1995).

embedding dimension (𝑚) 𝐷𝑚𝑎𝑥 𝑝−value
5 0.061 0.91
9 0.078 0.688

11 0.091 0.487

26



Table A.2: Statistical results of the two-sample KS test of 𝛿18𝑂𝑐 variability (𝑚 = 5) for every 25
years compared to the period after the industrial revolution (1970 - 1995).

period 𝐷𝑚𝑎𝑥 𝑝−value
1695 - 1719 0.138 0.923
1720 - 1744 0.178 0.721
1745 - 1769 0.218 0.483
1770 - 1794 0.191 0.627
1795 - 1819 0.226 0.424
1820 - 1844 0.186 0.664
1845 - 1869 0.111 0.987
1870 - 1894 0.138 0.923
1895 - 1919 0.092 0.999
1920 - 1944 0.226 0.424
1945 - 1969 0.302 0.140

Table A.3: Statistical results of the two-sample KS test of 𝛿18𝑂𝑐 variability (𝑚 = 9) for every 25
years compared to the period after the industrial revolution (1970 - 1995).

period 𝐷𝑚𝑎𝑥 𝑝−value
1695 - 1719 0.137 0.929
1720 - 1744 0.143 0.904
1745 - 1769 0.142 0.910
1770 - 1794 0.228 0.413
1795 - 1819 0.226 0.424
1820 - 1844 0.225 0.436
1845 - 1869 0.143 0.904
1870 - 1894 0.132 0.947
1895 - 1919 0.172 0.766
1920 - 1944 0.305 0.128
1945 - 1969 0.189 0.640

27



Table A.4: Statistical results of the two-sample KS test of 𝛿18𝑂𝑐 variability (𝑚 = 11) for every 25
years compared to the period after the industrial revolution (1970 - 1995).

period 𝐷𝑚𝑎𝑥 𝑝−value
1695 - 1719 0.112 0.985
1720 - 1744 0.143 0.904
1745 - 1769 0.103 0.994
1770 - 1794 0.225 0.436
1795 - 1819 0.223 0.448
1820 - 1844 0.192 0.616
1845 - 1869 0.143 0.904
1870 - 1894 0.209 0.554
1895 - 1919 0.135 0.936
1920 - 1944 0.343 0.063
1945 - 1969 0.228 0.413

28



(a)

(b)

(c)

Figure A.1: Annual 𝛿18𝑂𝑐 variability. Solid green line shows the median, while shaded region
denotes the central 95% quantiles from the 3×1000 random realizations for (a) 𝑚 = 5, (b) 𝑚 = 9, (c)
𝑚 = 11.

29



(a)

(b)

Figure A.2: Empirical cumulative probability density function (CDF) curves of the annual 𝛿18𝑂𝑐

variability (high-pass filtered 𝑚 = 5). The red line shows the period after the industrial revolution,
while the blue lines show the periods before the industrial revolution. We used two different two
different definitions of the anthropogenic era to show that the results are not sensitive to its definition
of the beginning of the exponential increase in the anthropogenic forcings. In (a) the industrial
revolution began in 1850, while in (b) the industrial revolution began in 1970.

30



(a)

(b)

Figure A.3: Empirical cumulative probability density function (CDF) curves of the annual 𝛿18𝑂𝑐

variability (high-pass filtered 𝑚 = 9). The red line shows the period after the industrial revolution,
while the blue lines show the periods before the industrial revolution. We used two different
definitions of the anthropogenic era to show that the results are not sensitive to its definition of the
beginning of the exponential increase in the anthropogenic forcings. In (a) the industrial revolution
began in 1850, while in (b) the industrial revolution began in 1970.

31



(a)

(b)

Figure A.4: Empirical cumulative probability density function (CDF) curves of the annual 𝛿18𝑂𝑐

variability (high-pass filtered 𝑚 = 11). The red line shows the period after the industrial revolution,
while the blue lines show the periods before the industrial revolution. We used two different
definitions of the anthropogenic era to show that the results are not sensitive to its definition of the
beginning of the exponential increase in the anthropogenic forcings. In (a) the industrial revolution
began in 1850, while in (b) the industrial revolution began in 1970.

32



Bibliography

[1] E. Aldrian and R. D. Susanto. Identification of three dominant rainfall regions within
Indonesia and their relationship to sea surface temperature. International Journal of
Climatology, 23(12):1435–1452, 2003.

[2] M. R. Allen and W. J. Ingram. Constraints on future changes in climate and the
hydrologic cycle. Nature, 419(6903):228–232, Sep 2002.

[3] K. J. Anchukaitis, M. N. Evans, T. Lange, D. R. Smith, S. W. Leavitt, and D. P.
Schrag. Consequences of a rapid cellulose extraction technique for oxygen isotope
and radiocarbon analyses. Analytical Chemistry, 80(6):2035–2041, 2008.

[4] M. Barbour, J. Roden, G. Farquhar, and J. Ehleringer. Expressing leaf water and
cellulose oxygen isotope ratios as enrichment above source water reveals evidence of
a Péclet effect. Oecologia, 138:426–35, 03 2004.

[5] E. Bauer, M. Claussen, V. Brovkin, and A. Hünerbein. Assessing climate forcings of
the Earth system for the past millennium. Geophysical Research Letter, 30, 03 2003.

[6] H. Bellenger, E. Guilyardi, J. Leloup, M. Lengaigne, and J. Vialard. ENSO represen-
tation in climate models: From CMIP3 to CMIP5. Climate Dynamics, 42:1999–2018,
2014.

[7] W. Cai, S. Borlace, M. Lengaigne, P. Rensch, M. Collins, G. Vecchi, A. Timmermann,
A. Santoso, M. McPhaden, L. Wu, M. England, G. Wang, E. Guilyardi, and F-F.
Jin. Increasing Frequency of Extreme El Niño Events due to Greenhouse Warming.
Nature Climate Change, 4, 01 2014.

[8] W. Cai, A. Santoso, M. Collins, B. Dewitte, C. Karamperidou, J-S. Kug, M. Lengaigne,
M. J. McPhaden, M. F. Stuecker, A. S. Taschetto, et al. Changing El Niño–Southern
Oscillation in a warming climate. Nature Reviews Earth & Environment, 2(9):628–
644, 2021.

33



[9] W. Cai, G. Wang, A. Santoso, M. McPhaden, L. Wu, F-F. Jin, A. Timmermann,
M. Collins, G. Vecchi, M. Lengaigne, M. England, D. Dommenget, K. Takahashi,
and E. Guilyardi. Increased frequency of extreme La Niña events under greenhouse
warming. Nature Climate Change, 5:132–137, 01 2015.

[10] M. Cane, A. Clement, A. Kaplan, Y. Kushnir, D. Pozdnyakov, R. Seager, S. Zebiak,
and R. Murtugudde. Twentieth-Century Sea Surface Temperature Trends. Science,
275:957, 02 1997.

[11] M. A. Cane. Forecasting El Niño with a Geophysical Model. In M. H. Glantz,
R. W. Katz, and N. Nicholls, editors, Teleconnections Linking Worldwide Climate
Anomalies., chapter 11, pages 345–369. Cambridge U. Press, Cambridge, England,
1991.

[12] D. Chen, M. A. Cane, A. Kaplan, S. E. Zebiak, and D. Huang. Predictability of El
Niño over the past 148 years. Nature, 428(6984):733–736, 2004.

[13] L. Chen, T. Li, Y. Yu, and S. K. Behera. A possible explanation for the divergent
projection of ENSO amplitude change under global warming. Climate Dynamics,
49:3799–3811, 2017.

[14] A. C. Clement, R. Seager, M. A. Cane, and S. E. Zebiak. An Ocean Dynamical
Thermostat. Journal of Climate, 9(9):2190 – 2196, 1996.

[15] Kim M Cobb, Christopher D Charles, Hai Cheng, and R Lawrence Edwards. El
Niño/Southern Oscillation and tropical Pacific climate during the last millennium.
Nature, 424(6946):271–276, 2003.

[16] M. Collins. El Niño or La Niña-like climate change? Clim. Dyn., 24:89–104, 2005.
doi: 10.1007/s00382-004-0478-x.

[17] M. Collins, S. I. An, W. Cai, A. Ganachaud, E. Guilyardi, F.-F. Jin, M. Jochum,
M. Lengaigne, S. Power, A. Timmermann, G. Vecchi, and A. Wittenberg. The impact
of global warming on the tropical Pacific Ocean and ENSO. Nature Geoscience,
3:391–397, 2010. doi: 10.1038/NGEO868.

[18] David Colquhoun. An investigation of the false discovery rate and the misinterpretation
of p-values. Royal Society Open Science, 1(3):140216, 2014.

[19] PAGES2k Consortium et al. A global multiproxy database for temperature reconstruc-
tions of the common era. Scientific data, 4, 2017.

[20] R. D’Arrigo, R. Wilson, J. Palmer, P. Krusic, A. Curtis, J. Sakulich, S. Bijaksana,
S. Zulaikah, L. O. Ngkoimani, and A. Tudhope. The reconstructed Indonesian warm
pool sea surface temperatures from tree rings and corals: Linkages to Asian monsoon
drought and El Niño–Southern Oscillation. Paleoceanography, 21(3), 2006.

34



[21] S. G. Dee, K. M. Cobb, J. Emile-Geay, T. R. Ault, R. L. Edwards, H. Cheng, and C. D.
Charles. No consistent ENSO response to volcanic forcing over the last millennium.
Science, 367(6485):1477–1481, 2020.

[22] H. F. Diaz, M. P. Hoerling, and J. K. Eischeid. ENSO variability, teleconnections and
climate change. International Journal of Climatology, 21(15):1845–1862.

[23] J. Emile-Geay, K. M. Cobb, J. E. Cole, M. Elliot, and F. Zhu. Past ENSO Variability,
chapter 5, pages 87–118. American Geophysical Union (AGU), 2020.

[24] M. N. Evans. Toward forward modeling for paleoclimatic proxy signal calibration:
A case study with oxygen isotopic composition of tropical woods. Geochemistry,
Geophysics, Geosystems, 8(7), 2007.

[25] M. N. Evans, K. J. Selmer, B. T. Breeden III, A. S. Lopatka, and R. E. Plummer.
Correction algorithm for online continuous flow 𝛿13C and 𝛿18O carbonate and cellulose
stable isotope analyses. Geochemistry, Geophysics, Geosystems, 17(9):3580–3588,
2016.

[26] M. N. Evans, S. E. Tolwinski-Ward, D. M. Thompson, and K. J. Anchukaitis. Appli-
cations of proxy system modeling in high resolution paleoclimatology. Quaternary
science reviews, 76:16–28, 2013.

[27] A. Gessler, J. P. Ferrio, R. Hommel, K. Treydte, R. A. Werner, and R. K. Monson. Sta-
ble isotopes in tree rings: towards a mechanistic understanding of isotope fractionation
and mixing processes from the leaves to the wood. Tree Physiology, 34(8):796–818,
06 2014.

[28] M. Ghil, M. R. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, W. Robert-
son, A, A. Saunders, Y. Tian, F. Varadi, et al. Advanced spectral methods for climatic
time series. Reviews of Geophysics, 40(1):3–1, 2002.

[29] A. Giannini, A. Robertson, and J. Qian. A role for tropical tropospheric temperature
adjustment to El Niño-Southern Oscillation in the seasonality of monsoonal Indonesia
precipitation predictability. Journal of Geophysical Research, 112, 08 2007.

[30] N. Golyandina, A. Korobeynikov, and A. Zhigljavsky. Singular Spectrum Analysis
with R. Springer, 2018.

[31] E. Guilyardi, A. Capotondi, M. Lengaigne, S. Thual, and A. Wittenberg. ENSO
Modeling: History, Progress, and Challenges, pages 199–226. 11 2020.

[32] Eric Guilyardi. El Niño–mean state–seasonal cycle interactions in a multi-model
ensemble. Climate Dynamics, 26:329–348, 2006.

35



[33] G. J. Hakim, J. Emile-Geay, E. J. Steig, D. Noone, D. M. Anderson, R. Tardif,
N. Steiger, and W. A. Perkins. The last millennium climate reanalysis project: Frame-
work and first results. Journal of Geophysical Research: Atmospheres, 121(12):6745–
6764, 2016.

[34] J-I. Hamada, M. Yamanaka, J. Matsumoto, S. Fukao, P. Winarso, and T. Sribimawati.
Spatial and Temporal Variations of the Rainy Season over Indonesia and their link to
ENSO. Journal of the Meteorological Society of Japan, 80:285–310, 04 2002.

[35] I. Harris, P. D. Jones, T. J. Osborn, and D. H. Lister. Updated high-resolution grids
of monthly climatic observations – the CRU TS3.10 Dataset. International Journal
of Climatology, 34:623–642, 2014. doi: 10.1002/joc.3711.

[36] H. Hassani and R. Mahmoudvand. Singular spectrum analysis: Using R. Springer,
2018.

[37] P. Hope, B. J. Henley, J. Gergis, J. Brown, and H. Ye. Time-varying spectral char-
acteristics of ENSO over the Last Millennium. Climate Dynamics, 49:1705–1727,
2017.

[38] H-Y Kao and J-Yi Yu. Contrasting Eastern-Pacific and Central-Pacific Types of ENSO.
Journal of Climate, 22(3):615 – 632, 2009.

[39] B-H. Kim and K-J. Ha. Observed changes of global and western Pacific precipitation
associated with global warming SST mode and Mega-ENSO SST mode. Climate
Dynamics, 45, 02 2015.

[40] L. Landrum, B. L. Otto-Bliesner, E. R. Wahl, A. Conley, P. J. Lawrence, N. Rosen-
bloom, and H. Teng. Last Millennium Climate and Its Variability in CCSM4. Journal
of Climate, 26(4):1085 – 1111, 2013.

[41] M. Latif and N. S. Keenlyside. El Niño/Southern Oscillation response to global
warming. Proceedings of the National Academy of Sciences, 106(49):20578–20583,
2009.

[42] S. C. Lewis and A. N. LeGrande. Stability of ENSO and its tropical Pacific telecon-
nections over the last millennium. Climate of the Past, 11(10):1347–1360, 2015.

[43] J. Li, S-P. Xie, E. R. Cook, M. S. Morales, D. A. Christie, N. C. Johnson, F. Chen,
R. D’Arrigo, A. M. Fowler, X. Gou, et al. El Niño modulations over the past seven
centuries. Nature Climate Change, 3(9):822–826, 2013.

[44] J. Liu, B. Wang, M. A. Cane, S-Y. Yim, and J-Y. Lee. Divergent global precipitation
changes induced by natural versus anthropogenic forcing. Nature, 493(7434):656–659,
Jan 2013.

36



[45] Y. Liu, K. M. Cobb, H. Song, Q. Li, C-Y. Li, T. Nakatsuka, Z. An, W. Zhou, Q. Cai,
J. Li, et al. Recent enhancement of central Pacific El Niño variability relative to last
eight centuries. Nature Communications, 8(1):15386, 2017.

[46] Z. Lu, Z. Liu, J. Zhu, and K. M. Cobb. A Review of Paleo El Niño-Southern
Oscillation. Atmosphere, 9(4):130, 2018.

[47] J. L. McBride, M. R. Haylock, and N. Nicholls. Relationships between the Maritime
Continent Heat Source and the El Niño–Southern Oscillation Phenomenon. Journal
of Climate, 16(17):2905 – 2914, 2003.

[48] W. J. Merryfield. Changes to ENSO under CO2 doubling in a Multimodel Ensemble.
Journal of Climate, 19(16):4009 – 4027, 2006.

[49] I. C. Nurhati, M. N. Evans, S. Y. Cahyarini, R. D. D’Arrigo, K. Y. Yoshimura, and
S. H. S. Herho. 𝛿18O of marine-influenced Tectona grandis l. f. from equatorial
Indonesia: A local rainfall amount and remote ENSO indicator. in preparation for
submission to Paleoceanography and Paleoclimatology.

[50] J.P. Peixoto and A.H. Oort. Physics of Climate. American Institute of Physics, New
York, 1991.

[51] N. Pumijumnong, P. Payomrat, S. Buajan, A. Bräuning, C. Muangsong, U. Chareon-
wong, P. Songtrirat, K. Palakit, Y. Liu, and Q. Li. Teak Tree-Ring Cellulose 𝛿13C,
𝛿18O, and Tree-Ring Width from Northwestern Thailand Capture Different Aspects
of Asian Monsoon Variability. Atmosphere, 12(6), 2021.

[52] C. S. Ramage. ROLE OF A TROPICAL “MARITIME CONTINENT” IN THE
ATMOSPHERIC CIRCULATION. Monthly Weather Review, 96(6):365 – 370, 1968.

[53] B. Rein, A. Lückge, L. Reinhardt, F. Sirocko, A. Wolf, and W-C. Dullo. El Niño
variability off Peru during the last 20,000 years. Paleoceanography, 20(4), 2005.

[54] A. Robertson, V. Moron, J. Qian, C-P. Chang, F. Tangang, E. Aldrian, T. Koh, and
L. Juneng. The Maritime Continent monsoon, pages 85–98. 04 2011.

[55] M. Rodriguez-Caton, L. Andreu-Hayles, V. Daux, M. Vuille, A. M. Varuolo-Clarke,
R. Oelkers, D. A. Christie, R. D’Arrigo, M. S. Morales, M. Palat Rao, A. M. Srur,
F. Vimeux, and R. Villalba. Hydroclimate and ENSO Variability Recorded by Oxygen
Isotopes From Tree Rings in the South American Altiplano. Geophysical Research
Letters, 49(4):e2021GL095883, 2022. e2021GL095883 2021GL095883.

[56] G. T. Rustic, A. Koutavas, T. M. Marchitto, and B. K. Linsley. Dynamical excitation
of the tropical Pacific Ocean and ENSO variability by Little Ice Age cooling. Science,
350(6267):1537–1541, 2015.

37



[57] D. H. Schoellhamer. Singular spectrum analysis for time series with missing data.
Geophysical Research Letters, 28(16):3187–3190, 2001.

[58] K. Schollaen, I. Heinrich, B. Neuwirth, P. J. Krusic, R. D. D’Arrigo, O. Karyanto,
and G. Helle. Multiple tree-ring chronologies (ring width, 𝛿13C and 𝛿18O) reveal
dry and rainy season signals of rainfall in Indonesia. Quaternary Science Reviews,
73:170–181, 2013.

[59] S. Stevenson, A. Capotondi, J. Fasullo, and B. Otto-Bliesner. Forced changes to
twentieth century ENSO diversity in a last millennium context. Climate Dynamics,
52:7359–7374, 2019.

[60] A. S. Taschetto, C. C. Ummenhofer, M. F. Stuecker, D. Dommenget, Karumuri A.,
R. R. Rodrigues, and S-W. Yeh. ENSO Atmospheric Teleconnections, chapter 14,
pages 309–335. American Geophysical Union (AGU), 2020.

[61] R. Vautard and M. Ghil. Singular spectrum analysis in nonlinear dynamics, with ap-
plications to paleoclimatic time series. Physica D: Nonlinear Phenomena, 35(3):395–
424, 1989.

[62] G. A. Vecchi and B. J. Soden. Global Warming and the Weakening of the Tropical
Circulation. Journal of Climate, 20(17):4316 – 4340, 2007.

[63] G. A. Vecchi, B. J. Soden, A. T. Wittenberg, I. M. Held, A. Leetmaa, and M. J.
Harrison. Weakening of tropical Pacific atmospheric circulation due to anthropogenic
forcing. Nature, 441(7089):73–76, 2006.

[64] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau,
E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wil-
son, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J.
Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cim-
rman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro,
F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental
Algorithms for Scientific Computing in Python. Nature Methods, 17:261–272, 2020.

[65] R. E. Walpole, R. H. Myers, S. L. Myers, and K. Ye. Probability and Statistics for
Engineers and Scientists, volume 5. Macmillan New York, 1993.

[66] J. Wang, J. Emile-Geay, D. Guillot, J. E. Smerdon, and B. Rajaratnam. Evaluat-
ing climate field reconstruction techniques using improved emulations of real-world
conditions. Climate of the Past, 10(1):1–19, 2014.

[67] D. S. Wilks. Statistical Methods in the Atmospheric Sciences, volume 100. Academic
Press, 2011.

38



[68] C. Xu, M. Sano, and T. Nakatsuka. Tree ring cellulose 𝛿18O of Fokienia hodginsii
in northern Laos: A promising proxy to reconstruct ENSO? Journal of Geophysical
Research: Atmospheres, 116(D24), 2011.

[69] S-W. Yeh, J-S. Kug, B. Dewitte, M-H. Kwon, B. P. Kirtman, and F-F. Jin. El Niño in
a changing climate. Nature, 461(7263):511–514, 2009.

[70] S. Yoden, S. Otsuka, N. J. Trilaksono, and T. W. Hadi. Recent Progress in Research
on the Maritime Continent Monsoon, chapter 6, pages 63–77. 2017.

[71] K. Yoshimura, M. Kanamitsu, D. Noone, and T. Oki. Historical isotope simulation
using Reanalysis atmospheric data. Journal of Geophysical Research: Atmospheres,
113(D19), 2008.

[72] Z. Zhongming, L. Linong, Y. Xiaona, L. Wei, et al. AR6 Climate Change 2021: The
Physical Science Basis. 2021.

[73] F. Zhu, J. Emile-Geay, K. J. Anchukaitis, G. J. Hakim, A. T. Wittenberg, M. S.
Morales, M. Toohey, and J. King. A re-appraisal of the ENSO response to volcanism
with paleoclimate data assimilation. Nature Communications, 13(1):747, 2022.

[74] F. Zhu, J. Emile-Geay, G. J. Hakim, R. Tardif, and A. Perkins. LMR Turbo (LMRt):
a lightweight implementation of the LMR framework, 2021.

[75] M. Zhu, L. Stott, B. Buckley, K. Yoshimura, and K. Ra. Indo-Pacific Warm Pool
convection and ENSO since 1867 derived from Cambodian pine tree cellulose oxygen
isotopes. Journal of Geophysical Research: Atmospheres, 117(D11), 2012.

39


	Acknowledgements
	Introduction
	Motivation
	Is ENSO changing? Evidence from climate models
	Is ENSO changing? Evidence from paleoclimatology
	This study
	Study location
	Data modeling: ^18O of -cellulose
	Goals of this study


	Materials and methods
	Stable isotope analysis
	Imaging, intra-annual sampling, and -cellulose extraction
	18Oc data acquisition

	Statistical analysis
	Time-series compositing
	Nonparametric statistical test


	Results
	Discussion
	Summary of results
	Comparison with prior studies

	Conclusion
	Statistical analyses for other embedding dimensions of the high-pass filtered time series
	Bibliography