during the last Millennia

coriandercultureΜηχανική

22 Φεβ 2014 (πριν από 3 χρόνια και 7 μήνες)

76 εμφανίσεις

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