1
Climate variability
in the Central Mediterranean
during the last Millennia
Dipartimento di Fisica Generale dell’Università di Torino and Istituto di Fisica dello Spazio Interplanetario
(IFSI), INAF, Torino, Italy
Gianna Vivaldo
Supervisor: Prof. Carla Taricco
Co

supervisor : Dott. Silvia Alessio
DOTTORATO DI RICERCA IN FISICA FONDAMENTALE, APPLICATA ED ASTROFISICA
XXI CICLO, 2005

2008
2
Study of past climatic variations
•
tree rings (dendrocronology)
•
corals (calcification)
•
ice cores (isotopic ratio)
•
sea and lake sediments cores
Very
different temporal scales
of climate variability (1
–
10
4
/10
5
years)
Direct measurements: ~ last 150 y
[1] Grove, J.M. and Switsur R.,
Climatic Change
, 26, 143
–
169, 1994; [2] Pollack, H., et al.,
Science
, 282, 279
–
281, 1998 [3] Jones, P. D.,
Springer

Verlag, Berlin, Germany, 649, 1996; [4]
Luckman, B.H., et al.,
The Holocene
, 7(4), 375
–
389, 1997; [5]
Crowley, T., J.,
Science
,
289(5477), 270
–
277, 2000
To go further back in the past we need indirect climate “proxy”
indicators measured in terrestrial archives
Knowledge of
natural climatic variability
Evaluation of
anthropogenic contribution
to recent climate change
3
SST reconstructions back to 1000 AD
combined with any available instrumental or historical records
[3, 4]
[6]
Hughes, M.K., and Diaz
,
H.F.
,
Climatic Change
, 26
,
109
–
142,
1994
; [7] Bradley, R.S., and Jones, P.D.,
The Holocene
, 3, 367
–
376, 1993; [8]
Jones, P.D. et al,
The Holocene
, 8, 477
–
483, 1998; [9] Mann, M.E., et al.,
Nature
, 392, 779
–
787, 1998
Locations of main proxy records with data back to AD 1000
Main limitations of available temperature reconstructions
merging of different records with different temporal resolutions
only few back to 1000 y ago and very few further back in time
tree archives reduce the amplitude of secular variations
need for long, homogeneous, well dated series with high time resolution
4
Ionian shallow

water sediment cores
downwind to the Campanian area
(historical documentation of major eruptions in the last 2000 y)
GALLIPOLI terrace
favourable site for climatic studies
flat region
no direct river discharge
weak and constant marine flows
5
Carbonate (CaCO
3
) profiles
•
uniformity of CaCO
3
profiles in different cores
uniform platform stratigraphy
•
reproducibility of the results
in the whole Gallipoli terrace
homogene
o
us sedimentation
Radiometric
210
Pb method
•
upper 20 cm of the core
•
constant sedimentation rate:
1 cm of mud deposited in
~ 15.5 y
137
Cs activity
•
maximum expected and found in 1963

64 AD (nuclear

weapon testing)
•
top layers undisturbed (no mixing during drilling)
Tephroanalysis
Using the sedimentation rate determined by
210
Pb method, peaks of volcanic material corresponding to
historical eruptions of the
Campanian area
in the last 2 millennia can be identified
•
confirms radiometric dating
•
extends dating to the last 2 millennia
•
allows determining sedimentation rate with higher accuracy
Dating of the cores
6
Tephrochronology
ejected into the
atmosphere
during an
explosive
eruption
[1] Morimoto et al. (1989)
(i.e.
quartz
,
olivines
,
feldspars
,
micas
)
crystallize
before
lava
erupted
sometimes
very
far
from
the
volcano
vent
Pyroxenes crystals
Event

marker
:
number density of pyroxenes (ubiquitous in these lavas)
[1]
Historical documentation of major eruptions available for the last 2 millenia
allows precise dating
•
GT89/3
core sampled at
2.5 mm
in a continuous sequence
•
0.5 mg of sediment for each sample treated (CaCO
3
removed) and analyzed under a
microscope
(sector zoning and skeletal morphology strongly suggest strong volcanic origin)
Experimental procedure
1 mm
7
Time
–
depth relation
Time

depth
relationship
:
the
depth
at
which
a
volcanic
peak
(
red
point
)
is
found
in
the
sediment
is
plotted
versus
the
historical
date
of
the
corresponding
eruption,
expressed
in
years
counted
backwards
from
1979
A
.
D
.
,
i
.
e
.
the
date
of
the
core
top
(hence,
years

before

top)
.
The
linear
regression
fit
to
the
experimental
data
(
straight
blue
line
)
is
also
shown
;
its
high
correlation
coefficient
(
r
=
0
.
99
)
attests
to
the
goodness
of
the
fit
.
Constant in the last 2 Millennia
Uniform throughout the whole platform
D
z = 2,5 mm
D
t = 3,87
anni
pyroxenes peaks (
22
eruptions)
linear regression fit (
r
2
= 0.99
)
79 AD
1944 AD
core top
1979 AD
•
1638

1944:
catalogue by Arnò et al.
•
ante
1638:
sparse documentation
depth = (0.0645
±
0.0002) y
BT
slope
sedimentation rate
1 cm of mud
deposited in
~
15.5 y
[2] Arnò et al., 1987 (
Somma Vesuvius
–
Eruptive History)
8
Pyroxenes time series : 20 AD
–
1979 AD
POMPEI
(79 AD)
POLLENA
(472 AD)
ISCHIA
(1301 AD)
506 samples
–
D
t
=
3.87 y
9
Objective testing of dating procedure
In
the
frame
of
the
European
Project
on
Extreme
Events
–
Causes
and
Consequences
(E
2
C
2
)
in
which
our
group
is
involved,
I
applied
a
statistical
method,
allowing
the
objective
identification
of
pulse
like

events
in
a
series
,
to
extract
the
volcanic
peaks
from
the
long
pyroxenes
series
measured
in
the
shallow

water
core
GT
89

3
.
The
work
has
recently
been
published
on
Nonlinear
Processes
in
Geophysics
.
Timing and intensity of events and the associated
posterior probability are determined by fitting to
the series a statistical
model including a slowly changing component or
baseline, an abruptly changing component
representing the volcanic signal and a white noise
background.
Sequence of eruptive events in the Vesuvio area recorded in shallow

water Ionian Sea sediments
C. Taricco, S. Alessio, and G. Vivaldo, NPG, in press
10
Automatic extraction of volcanic events
[3] Guo, W., et al.,
Journal of the American Statistical Society,
94(447), 746

756, 1999
irregular abruptly changing component (
pulses
)
slowly changing component (
baseline
)
eruptions
(timing
and
intensity)
objectively
identified
as
significant
at
given
probability
threshold
posterior
probability
not
necessarily
determined
only
by
absolute
spike
intensity
but
linked
mainly
to
signal
local
characteristics
[
3
]
pulses
extracted
from
variable
background
(unlike
in
conventional
thresholding
techniques)
pulse
homogeneity
assumption
Stochastic model
y
(
t
) = input time series
x
(
t
) = volcanic signal (pulses)
f
(
t
) = long

term trend (baseline)
(
t
) = Gaussian white noise,
)
,
0
(
)
(
2
e
N
t
where
a
<1 = constant expressing
expon. decay with time of the effect of
an eruption
v
(
t
) = volcanic input =
random variable modeling the strength
of an eruption
Introduce a Bernouilli variable o(t)
whose parameter
= probability( o(t) =
1) ) denotes the
a priori
probability of
observing an eruption at time t
if o(t) = 1 (eruption), v(t) =
random gaussian variable;
if o(t) = 0 (no eruption), v(t) = 0.
cubic smoothing spline
multi

process dynamic linear model
[8]
The model assumes that the pulses are homogeneous, while the
pyroxenes series contains two peaks much higher than the others
)
(
)
1
(
)
(
t
v
t
ax
t
x
decay
l
exponentia
1
a
0
)
(
if
0
1
)
(
if
,
)
(
2
t
o
t
o
N
t
v
x
x
eruption
no eruption
)
(
)
(
)
(
)
(
t
t
f
t
x
t
y
linear
combination
of
3
independent
components
simultaneous
estimation
of
the
multiple
outputs
–
Multi
Process
State

Space
Model
•
observational

level
equation
–
state

vector
+
random
Gaussian
noise
•
system

level
equation
–
Markov
transition
equation
for
the
state
vector
–
structural
parameters
definition
a dynamic system is not
directly observable
without ERROR
y
(
t
) = input time series
x
(
t
) = volcanic signal (pulses)
f
(
t
) = long

term trend (baseline)
(
t
) = Gaussian white noise
11
Model outputs and parameters estimation
Polinomial smoothing spline (cubic)
–
f(t)
1.
flexibility
–
time dependent trend
2.
straight line + random perturbation (smoothing deviation)
3.
f(t)
and the
m

1
derivatives
–
state

space representation
[4]
(m=2
cubic
)
y(t) = input time series
x(t) = volcanic signal (pulses)
f(t) = long

term trend (baseline)
(t) = Gaussian white noise,
)
,
0
(
)
(
2
e
N
t
)
(
)
(
)
(
)
(
t
t
f
t
x
t
y
[4] Kushler and Brown, 1991
Estimated Parameters:
x(t)
halflife (exponential decay
–
effect of an eruption)
v(t)
mean and variance
(t)
mean
f(t)
smoothing parameter
0
)
(
if
0
1
)
(
if
,
)
(
2
t
o
t
o
N
t
v
x
x
eruption
no eruption
Multi

process dynamic linear model
–
x(t)
)
(
)
1
(
)
(
t
v
t
ax
t
x
decay
l
exponentia
1
a
input
volcanic
v(t)
12
Application to the pyroxenes series
The method does not work well with a series that has such huge peaks as the
ones marking
Pompei
and
Ischia
eruptions
I applied the model to two separate segments of the series:
•
150
–
1250 AD
•
1350
–
1979 AD
13
Good model fit
150
–
1250 AD and 1350
–
1979 AD
14
Detail of model fit to pyroxenes series (last 300 y)
Time AD (years)
15
18 cycles (
13

14 y
of duration on average), each one ending with an
eruption from moderate to high explosive activity
mean dormancy time (
3 y
)
Vesuvius eruptive history documented by a detailed catalogue
[2]
Volcanic activity 1638
–
1944 AD
16
Volcanic activity 1638
–
1944 AD
Results at
99% probability threshold
Historical
[AD]
From
dating
[AD]
1660
1657.79
1682
1684.88
1698
1700.36
1794
1793.24
1822
1824.20
1906
1905.47
1944
1944.17
All the 7 catalogue

documented
events used for dating are
recognized
One extra

eruption documented in the
catalogue but internal to a cycle is found
(1810 AD)
17
Volcanic activity
ante
1638
10/13
eruptions used for dating
found (+, )
(Pompei and Ischia excluded)
3
historical events revealed
at lower probability level
(
+
, )
4
possible undocumented
events revealed ( )
346 AD, 914 AD, 1388 AD,
1472 AD
Results at 99% probability threshold
18
Spectral analysis
TRENDS
:
comparison
between
model
baseline
and
low

frequency
components
VARIABILITY
:
main
cycles
embedded
in
the
series
ADVANCED SPECTRAL METHODS
MultiTaper Method (MTM)
Singular Spectrum Analysis (SSA), Monte

Carlo SSA (MC

SSA)
Wavelet Transform (WT)
19
Multi

Taper Method (MTM)
[5] Thomson, D., J.,
Proceedings of IEEE
, 70, 1055
–
1096, 1982; [6] Ghil, M., et al.,
Reviews of Geophysics
, 40(1), 1
–
41, 2002
Averaging over this
(small) ensemble of
spectra (so K
determines the
stability)
Reshaped
spectrum:
contributions from
harmonic signals
removed
Leakage:artificiall
y high power
estimates at
frequencies away
from the true peak
frequency
resolution
spectral
1
where
,
D
c
R
R
hb
NT
f
pf
f
Power
Spectral
Density
(PSD)
estimate
+
Signal
Reconstruction
a)
harmonic
or
anharmonic
narrow

band
signal
(singular/discrete)
b)
background
noise
(continuous)
MTM
Spectrum
:
weighted
average
of
a
small
number
of
independent
spectra
Objectively optimized small set of
k
orthogonal TAPERS (series pre

filtering
[5]
)
Trade off between spectral resolution and stability (variance reduction)
Two ways of testing the MTM spectrum
1.
Harmonic
analysis
: estimate periodic components and their amplitude
2.
Red

noise test
[6]
: detect non

harmonic outliers (robust estimate)
error bars no more proportional to the peaks amplitude
20
MTM
–
red noise background
Climatic time series: the
intrinsic inertia of the system leads to
greater power at low frequencies,
even in the absence of any signal
[7]
[7] Hasselmann, K., Tellus, 6, 473
–
485, 1976
21
MTM
–
locally

white noise background
Colored noise process that varies
slowly but arbitrarily with
frequency (complex structure)
22
Singular Spectrum Analysis (SSA)
Short,
noisy
and
chaotic
time
series
Non
parametric
method
Empirical
Orthogonal
Functions
(EOFs)
:
data
adaptive
basis
functions

adaptive
filters
Principal Components (PCs)
: series projections on the new basis
Reconstructed Components (RCs)
: from a set of k components
filtered version of the series
SSA spectrum: sum of the PCs power spectra
mode
th
of
variance
(EOF)
r
eigenvecto
th

matrix
covariance

lag
k
k
k
k
C
C
k
k
k
•
a
ssess whether the SSA spectral estimation
can reject a
null hypothesis (normally red noise)
•
gives
statistical significance
of the spectral components
Monte Carlo SSA (MC

SSA)
[6]
23
Model evaluation by Monte

Carlo SSA
Null

hyp: [AR(1) + trend (RCs 1, 2) +
+ 400 y (RCs 3, 4) + 13.8 y (RCs 40, 41)
+ 9.4 y (RCs 53, 54) + 8.4 y (RCs 65, 66)]
Confidence level 90%
Ensemble size 10000
M=170
200
100
M
Similar results with
24
Wavelet Transform (WT)
Scalogram
: square modulus of WT, PSD is a function of time and scale (
Fourier period)
(2D contour plot)
Global Wavelet Spectrum (GWS)
: time average of scalogram values at each Fourier
period
Non stationary time series
Evolutionary spectral analysis
:
dominant periodic signals which can vary both
in amplitude and frequency over record duration
Discretization
of the
Continuous
Wavelet
Transform:
•
all time steps
•
dense set of
scales ( integer
and fractional
powers of two).
For
cmor,
scale
Fourier
period.
as well as
smooth
estimates of
stationary
(time
averaged)
spectra
Morlet wavelet
Window function: wave packet of finite duration and
specific frequency
Correlation between the time series and the chosen wavelet
at each time along the series and at several “scales” (wavelet stretched and shifted)
1.
confidence level 90%
2.
background spectrum of red noise
[8] Torrence, C. and Compo, G., P., Bulletin of American Meteorological Society, 1998
Statistical tests
[8]
25
Scalogram and GWS of pyroxenes series

areas of statistically significant power (90% c.l.)

cone
of influence
:
spectral power outside the cup is reduced with
respect to hypothetical “true” value, due to zero padding of series performed while
computing WT in the frequency domain.
Wavelet location in time
10%
significance
regions (red
noise
background)
temporal resolution
frequency resolution
Series preliminarily detrended subtracting contribution by scales >1260 y
as suggested by spectral break at that period in original series global spectrum (not shown)
26
Global Wavelet Spectrum of pyroxenes series
Good agreement
with MTM and SSA
27
Low frequency components comparison:
trend + 400 y variability and baseline
28
Further developments
1.
assumption
in
statistical
model
:
signals
magnitudes
follow
a
Gaussian
distribution
2.
clear
and
distinct
signature
in
the
sediment
:
largest
events
volcanic
spike
amplitudes
should
follow
a
heavy
skewed
distribution
:
•
Gaussian fit: values outside the tail
•
Exponential fit
[10]
:
not able to capture the largest peaks strength
3.
General
Pareto
Distribution
(GPD)
:
asymptotical
approximation
of
exceedances
distribution
3.
seems to be the better one for characterize this kind of events.
4.
Exp, uniform, cauchy
P
(
X
x
)
1
1
x
t
1
x
,
t
\
x
t
0
1
x
t
0
0
t = threshold
= scale parameter
= shape parameter
Gumbel distribution
–
light tailed
Frèchet distribution
–
heavy tailed
Weibull distribution
–
upper bounded
0
0
0
GPD is able to fit exceedances from any continuous distribution
general methodology
Study of the statistics of volcanic events in the frame of EVT theory
Possible extention of the statistical model for pulse extraction to any amplitude range
29
18
O isotopic analysis
In
the
cores
dated
in
the
way
described
before,
we
measured
and
analyzed
18
O
in
relation
with
temperature
in
the
last
millennia
18
O
is
measured
in
the
carbonatic
shells
of
foraminifera
(
Globigerinoides
ruber
)
present
in
the
sediments
30
Core GT90

3
39
°
45'53''N 17
°
53'33''E, depth 174 m, length 3.57 m
Ionian Sea, Gulf of Tarant
o
1.
Monte dei Cappuccini Laboratory
–
core sampling
2.
Cosmogeophysics Laboratory
–
chemical treatment
3.
ETH Zurich
–
isotopic measurements
31
Monte dei Cappuccini Laboratory

Torino
•
Cores stored in refrigerating
cell at
4
°
C
•
core sampling at 2.5 mm
interval (continuous
sequence)
•
5

6 g of sediment for
analysis
32
DFG Cosmogeophysics Laboratory

Torino
5%
Calgon solution
overnight (esametaphosphate (
NaPO
3
)
6
solution)
10%
H
2
O
2
attack
(organic material)
distilled water jet washing through a
sieve of 150
μm
filtered fraction kept and
oven

dried at 50
°
C
G. Ruber
specimens picked up under a
microscope
(20
–
30/samples)
Cleaning of the
shell carbonate
ensure pure CO
2
to
be analyzed
33
Globigerinoides Rubes selection
Sample after chemical treatment
34
ETH Zurich
Mass spectrometer
3. Ionized CO
2
+
isotopic components selectively deflected by electric and
magnetic field (p =
10

8
mbar
) and revealed by
3 Faraday collectors
4.
18
O isotopic ratio determination
5.
Calibration to the VPDB standard using LSI/SIL carbonate standards (
Carrara
Marmor
)
C
T
PO
Ca
O
H
CO
PO
H
CaCO
90
at
3
3
2
3
2
4
3
2
2
4
3
3
Preparation Device
2.
Successive freezing steps to trap
H
2
O
and
contaminants
(N
2
, O
2
)
CO
2
fed to spectrometer
1.
Reaction with phosphoric acid (pirolysis):
35
Experimental apparatus
36
Calibration of the paleotemperature scale
•
SMOW
= Standard Mean Ocean Water
•
VPDB
= Vienna Pee Dee Belemnite (South Carolina Cretaceous formation)
RATIO of the
heavy to the
light isotope in
molecule of
phase A /(stessa
cosa in phase B)
“expressed as
per mil
difference
relative to
SMOW”
Oxygen isotopic composition expressed as a
‰ difference relative to a standard (VPDB)
3
16
18
/
16
18
/
18
10
1
)
/
(
)
/
(
SMOW
w
c
w
c
O
O
O
O
O
86
.
30
03
.
1
18
18
VPDB
SMOW
O
O
2
)
(
10
.
0
)
(
38
.
4
9
.
16
w
c
w
c
T
Shackleton equation
18
O
w
in the past is impossible to evaluate
18
O in the calcite forming the shells of foraminifera
depends on T
and isotopic oxygen ratio of
ambient water
we resort to the comparison of our
18
O record with
global T reconstructions for NH
in order to extract T information, even in the presence of this unknown parameter
37
18
O time series
D
t = 3.87 y
N = 560
38
18
O
significant spectral components from SSA
125
y
200
y
350
y
600
y
Trend
a)
b)
11
y
6 highly significant
spectral components
(42% of total
variance is described)
In phase with 11 y solar cycle
M=150
39
Sum of significant centennial components
––
raw data
––
sum of components
18
O depends also on
18
O
w
not
all the significant components extracted are
necessarily linked to T
In order to overcome this problem we compared our series with
T reconstructions available over the last 1000 y
40
Comparison with NH SST reconstructions
R
econstructed
long

term
variability
of
NH
temperature
by
Mann
et
al
.
[
6
]
(RCs
1
—
4
,
red
)
and
the
corresponding
18
O
long

term
variability
component
(RC
1
+
RCs
6
—
8
,
black
),
original
temperature
series
of
Mann
et
al
.
[
6
]
(
gray
)
TRENDS
Mann
et
al
.
[
9
]
SST
reconstruction
Jones
et
al
.
[
10
]
SST
time
series
18
O
profile
200 y OSCILLATION
Mann et al.
[9]
SST
reconstruction
18
O profile
Same SSA analysis on
several T reconstructions
Good match between
trend and 200 y components
extracted from
18
O
series and those extracted from T reconstructions
[9]
Mann et al.,
Geoph. Res. Lett.
, 26(6), 759, 1999;
[10]
Jones et al.,
The Holocene
8(4), 455, 1998.
Mediterranean area:
stable region but not adequate for general conclusions about global changes
Global fluctuations:
they may leave a mark also at local scales
41
Long term temperature variability
1.
Linear regression (r =

0.7)
D
T
NH

ㄮ㘠
S
(
S
=
18
O
long

term oscillations)
2.
0
–
1100 AD (
Medieval Optimum
, 900
–
1200 AD)
steep temp. rise (+0.4
°
C
–
D
18
O
=0.23
‰
)
3.
Sp
ö
rer
(1500 AD)
–
Maunder
(1700 AD)
–
Modern
(1900 AD) minima (low solar activity)
Calibration in terms of NH temperature
4.
450 AD and 650 AD: less pronounced
temperature minima
5.
shallow local minimum at the
Little Ice Age
(
LIA
, 1300
–
1850 AD)
6.
Minimum at 0 AD
•
Mann
et
al
.
[
6
]
NH
temperature
reconstruction
(RCs
1
–
4
)
•
18
O
time
series
(RC
1
+
RCs
6
—
8
)
42
Further developments
[] Pirazzoli P., A., Science, 194 (4624), 519
–
521, 1976
Comparison with alkenones time series recently measured in the same cores
(SST proxy)
Extension of the
18
O
series further back in the past.
43
Results
Pyroxenes time series
Pulse

like events automatic extraction
(high posterior probability):
ante

1638: the majority of the eruptions used for dating and a few more possible
events
post

1638: all the eruptions used for dating
Spectral analysis
: trend + 400 y + high frequency components (<15 y period)
18
O time series
High resolution climatic time series over 2200 y
Spectral analysis
: trend + 600 y + 350 y + 200 y + 125 y + 11 y
NH long

term temperature variations well described by trends and 200 y components
this allows obtaining
information about T before 1000 AD
Minimum of T revealed at 0 AD
44
Carbonate profiles
45
Equivalent expressions of isotopic equilibrium
molecule
reactive
more
isotopes
lighter
2
,
1
,
2
1
m
m
1
E
E
m
vibr
Isotopic
fractionation
:
during
a
chemical

physical
process
(thermodynamics)
the
isotopes
are
discriminated
by
their
mass
39
.
3
10
78
.
2
ln
1000
2
6
T
w
c
O
O
O
O
O
H
O
H
CaCO
CaCO
O
H
CaCO
O
H
CaCO
k
16
18
16
18
16
2
18
2
3
1
16
3
18
3
18
2
3
1
16
3
16
2
3
1
18
3
)
]
[
]
[
(
)
]
[
]
[
(
]
[
]
[
]
[
]
[
1000
1000
)
(
18
18
w
c
O
O
T
k
2
)
(
10
.
0
)
(
38
.
4
9
.
16
w
c
w
c
T
2
)
(
13
.
0
)
(
2
.
4
9
.
16
w
c
w
c
T
Isotopic equilibrium in Benthonic foraminifera
Uvigerina
(C
–
Shakleton)
(A
–
Craig)
)
(
0
.
4
9
.
16
w
c
T
Commentaires 0
Connectezvous pour poster un commentaire