Transcript Presentazione di PowerPoint - Istituto Nazionale di Fisica

Dipartimento di Fisica Generale dell’Università di Torino and Istituto di Fisica dello Spazio Interplanetario
(IFSI), INAF, Torino, Italy
DOTTORATO DI RICERCA IN FISICA FONDAMENTALE, APPLICATA ED ASTROFISICA
XXI CICLO, 2005-2008
Climate variability
in the Central Mediterranean
during the last Millennia
Gianna Vivaldo
Supervisor: Prof. Carla Taricco
Co-supervisor : Dott. Silvia Alessio
1
Study of past climatic variations
 Knowledge of natural climatic variability
 Evaluation of anthropogenic contribution to recent climate change
 Very different temporal scales of climate variability (1– 104/105 years)
 Direct measurements: ~ last 150 y
To go further back in the past we need indirect climate “proxy”
indicators measured in terrestrial archives
• tree rings (dendrocronology)
• corals (calcification)
• ice cores (isotopic ratio)
• sea and lake sediments cores
2
SST reconstructions back to 1000 AD
Locations of main proxy records with data back to AD 1000
Main limitations of available temperature reconstructions
 needformerging
long, homogeneous,
well
dated
series
with high
time resolution
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
3
Ionian shallow-water sediment cores
GALLIPOLI terrace
favourable site for climatic studies
 flat region

no direct river discharge

weak and constant marine flows
downwind to the Campanian area
(historical documentation of major eruptions in the last 2000 y)
4
Dating of the cores
 Carbonate (CaCO3) profiles
• uniformity of CaCO3 profiles in different cores uniform platform stratigraphy
• reproducibility of the results in the whole Gallipoli terrace homogeneous sedimentation
 Radiometric 210Pb method
• upper 20 cm of the core
• constant sedimentation rate: 1 cm of mud deposited in ~ 15.5 y
 137Cs 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 210Pb 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
5
Tephrochronology
Event-marker: number density of pyroxenes (ubiquitous in these lavas)[1]
1 mm
Pyroxenes crystals
Experimental procedure
• GT89/3 core sampled at 2.5 mm in a continuous sequence
• 0.5 mg of sediment for each sample treated (CaCO3 removed) and analyzed under a
microscope (sector zoning and skeletal morphology strongly suggest strong volcanic origin)
Historical documentation of major eruptions available for the last 2 millenia
allows precise dating
[1] Morimoto et al. (1989)
6
Time – depth relation
pyroxenes peaks (22 eruptions)
linear regression fit (r2 = 0.99)
79 AD
• 1638-1944:
catalogue by Arnò et al.
depth = (0.0645 ± 0.0002) yBT
• ante 1638:
sparse documentation
1944 AD
slope  sedimentation rate
core top  1979 AD
[2] Arnò et al., 1987 (Somma Vesuvius – Eruptive History)
1 cm of mud
deposited in ~
15.5 y
7
Pyroxenes time series : 20 AD – 1979 AD
506 samples – t = 3.87 y
ISCHIA
(1301 AD)
POMPEI
(79 AD)
POLLENA
(472 AD)
8
Objective testing of dating procedure
In the frame of the European Project on Extreme Events – Causes and
Consequences (E2C2) 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 GT89-3.
The work has recently been published on Nonlinear Processes in Geophysics.
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
9
Automatic extraction of volcanic events
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 )  x(t )  f (t )   (t )
y(t) = input time series
x(t) = volcanic signal (pulses)
f(t) = long-term trend (baseline)
(t) = Gaussian white noise
 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
[3] Guo, W., et al.,Journal of the American Statistical Society, 94(447), 746-756, 1999
10
Model outputs and parameters estimation
y(t )  x(t )  f (t )   (t )
Estimated Parameters:
y(t) = input time series
x(t) = volcanic signal (pulses)
f(t) = long-term trend (baseline)
(t) = Gaussian white noise,  (t )  N (0, e2 )
x(t)  halflife (exponential decay – effect of an eruption)
v(t)  mean and variance
(t)  mean
f(t)  smoothing parameter
Polinomial smoothing spline (cubic) – f(t)
1.
2.
3.
flexibility – time dependent trend
straight line + random perturbation (smoothing deviation)
f(t) and the m-1 derivatives – state-space representation [4] (m=2  cubic)
Multi-process dynamic linear model – x(t)
x(t )  ax(t 1)  v(t )
a  1 exponential decay



 N  x , x2 if o(t )  1
v(t )  

if o(t )  0
0
eruption
no eruption
v(t)  volcanic input
[4] Kushler and Brown, 1991
11
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
12
Good model fit
150 – 1250 AD and 1350 – 1979 AD
13
Detail of model fit to pyroxenes series (last 300 y)
Time AD (years)
14
Volcanic activity 1638 – 1944 AD
Vesuvius eruptive history documented by a detailed catalogue [2]
 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)
15
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 extra-eruption
the 7 catalogue-documented
One
documented in the
events
for dating
areis found
catalogue
butused
internal
to a cycle
(1810 AD)
recognized
16
Volcanic activity ante 1638
Results at 99% probability threshold
 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
17
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)
18
Multi-Taper Method (MTM)

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
 error bars no more proportional to the peaks amplitude
2.
Red-noise test [6]: detect non-harmonic outliers (robust estimate)
[5] Thomson, D., J., Proceedings of IEEE, 70, 1055 – 1096, 1982; [6] Ghil, M., et al., Reviews of Geophysics, 40(1), 1 –41, 2002
19
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
20
MTM – locally-white noise background
Colored noise process that varies
slowly but arbitrarily with
frequency (complex structure)
21
Singular Spectrum Analysis (SSA)

Short, noisy and chaotic time series

Non parametric method

Empirical Orthogonal Functions (EOFs): data adaptive basis functions - adaptive filters
C k  k  k
 C  lag - covariancematrix

 k  k - th eigenvector (EOF)
   varianceof k  th mode
 k
 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
Monte Carlo SSA (MC-SSA)[6]
• assess whether the SSA spectral estimation can reject a null hypothesis (normally red noise)
• gives statistical significance of the spectral components
22
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)]
M=170
Similar results with
100  M  200
Confidence level 90%
Ensemble size 10000
23
Wavelet Transform (WT)
 Non stationary time series
 Evolutionary spectral analysis:
dominant periodic signals which can vary both
in amplitude and frequency over record duration
 Window function: wave packet of finite duration and
specific frequency
Morlet wavelet
 Correlation between the time series and the chosen wavelet
at each time along the series and at several “scales” (wavelet stretched and shifted)
 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
Statistical tests [8]
1.
2.
confidence level 90%
background spectrum of red noise
[8] Torrence, C. and Compo, G., P., Bulletin of American Meteorological Society, 1998
24
Scalogram and GWS of pyroxenes series
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)
25
Global Wavelet Spectrum of pyroxenes series
Good agreement
with MTM and SSA
26
Low frequency components comparison:
trend + 400 y variability and baseline
27
Further developments
1.
2.
assumption in statistical model: signals magnitudes follow a Gaussian distribution
clear and distinct signature in the sediment: largest events
 volcanic spike amplitudes should follow a heavy skewed distribution:
3.
General Pareto Distribution (GPD): asymptotical approximation of exceedances
distribution
1


x  t  
P(X  x)  1 1 


 

t = threshold
 = scale parameter
 = shape parameter
 0
   0
 0


x  t 
x,t   \ x  t   0  1 
 0


 
  0

Gumbel distribution – light tailed
Frèchet distribution – heavy tailed
Weibull distribution – upper bounded
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
28
18
 O
isotopic analysis
In the cores dated in the way described before, we measured and
analyzed 18O in relation with temperature in the last millennia
18O is measured in the carbonatic shells of foraminifera
(Globigerinoides ruber) present in the sediments
29
Core GT90-3
39°45'53''N 17°53'33''E, depth 174 m, length 3.57 m
Ionian Sea, Gulf of Taranto
1.
Monte dei Cappuccini Laboratory – core sampling
2.
Cosmogeophysics Laboratory – chemical treatment
3.
ETH Zurich – isotopic measurements
30
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
31
DFG Cosmogeophysics Laboratory - Torino





5% Calgon solution overnight (esametaphosphate (NaPO3)6 solution)
10% H2O2 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)
32
Globigerinoides Rubes selection
Sample after chemical treatment
33
ETH Zurich
Preparation Device
1.
Reaction with phosphoric acid (pirolysis):
3CaCO3  2H3PO4  3CO2  3H 2O  Ca3 PO4 2
2.
at T  90C
Successive freezing steps to trap H2O and contaminants (N2, O2)
 CO2 fed to spectrometer
Mass spectrometer
3.
4.
5.
Ionized CO2+  isotopic components selectively deflected by electric and
magnetic field (p = 10-8 mbar ) and revealed by 3 Faraday collectors
 18O isotopic ratio determination
Calibration to the VPDB standard using LSI/SIL carbonate standards (Carrara
Marmor)
34
Experimental apparatus
35
Calibration of the paleotemperature scale
Oxygen isotopic composition expressed as a ‰ difference relative to a standard (VPDB)
 Oc / w
18
 (18O /16O)c / w

  18 16
 1  103
 ( O / O) SMOW

 18OSMOW  1.03 18OVPDB  30.86
Shackleton equation T  16.9  4.38(c   w )  0.10(c   w )2
18O in the calcite forming the shells of foraminifera depends on T and isotopic oxygen ratio of
ambient water
18Ow in the past is impossible to evaluate
we resort to the comparison of our 18O record with global T reconstructions for NH
in order to extract T information, even in the presence of this unknown parameter
• SMOW = Standard Mean Ocean Water
• VPDB = Vienna Pee Dee Belemnite (South Carolina Cretaceous formation)
36
18O time series
t = 3.87 y
N = 560
37
18O significant spectral components from SSA
a)
125 y
M=150
200 y
350 y
600 y
6 highly significant
spectral components
(42% of total
variance is described)
Trend
b)
11 y
In phase with 11 y solar cycle
38
Sum of significant centennial components
–– raw data
–– sum of components
18O depends also on 18Ow  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
39
Comparison with NH SST reconstructions
 Same SSA analysis on several T reconstructions
 Good match between trend and 200 y components extracted from 18O
series and those extracted from T reconstructions
TRENDS
Mann
et
al.[9]
SST
reconstruction
Jones et al.[10] SST time series
18O profile
200 y OSCILLATION
Mann et al.[9] SST
reconstruction
18O profile
[9] Mann et al., Geoph. Res. Lett., 26(6), 759, 1999; [10] Jones et al., The Holocene 8(4), 455, 1998.
40
Long term temperature variability
1.
Linear regression (r = - 0.7)  TNH  - 1.6 S
4.
(S = 18O long-term oscillations)
2.
0 – 1100 AD (Medieval Optimum, 900 – 1200 AD)
5.
 steep temp. rise (+0.4 °C – 18O =0.23‰)
3.
Spörer (1500 AD) – Maunder (1700 AD) – Modern
6.
450 AD and 650 AD: less pronounced
temperature minima
shallow local minimum at the Little Ice Age
(LIA, 1300 – 1850 AD)
Minimum at 0 AD
(1900 AD) minima (low solar activity)
41
Further developments
 Comparison with alkenones time series recently measured in the same cores
(SST proxy)
 Extension of the 18O series further back in the past.
[] Pirazzoli P., A., Science, 194 (4624), 519 – 521, 1976
42
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)
18O 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
43
Carbonate profiles
44
Equivalent expressions of isotopic equilibrium
T  16.9  4.0( c   w )
(C – Shakleton)
T  16.9  4.38( c   w )  0.10( c   w ) 2
(A – Craig)
T  16.9  4.2( c   w )  0.13( c   w ) 2
Isotopic equilibrium in Benthonic foraminifera Uvigerina
45