%%%\documentclass[cpd, online, hvmath]{copernicus_discussions}

\documentclass[cpd, online, hvmath, natbib]{copernicus}

\usepackage{color}

\begin{document}\sloppy

\title{Potential analysis reveals changing number of climate states during
the last 60\,kyr}

\author[1]{V.~N.~Livina}
\author[2]{F.~Kwasniok}
\author[1]{T.~M.~Lenton}

\affil[1]{School of Environmental Sciences, University of
East Anglia, Norwich, UK}

\affil[2]{School of Engineering, Mathematics and Physical Sciences,
University of Exeter, Exeter, UK}

\correspondence{V.~N.~Livina (v.livina@uea.ac.uk)}

\runningtitle{Potential analysis}
\runningauthor{V.~N.~Livina et al.}

\received{28 August 2009}
\accepted{9 September 2009}
\published{}

\firstpage{1}

\maketitle

\begin{abstract}
We develop and apply a new statistical method of potential analysis
for detecting the
number of states of a geophysical
system, from its recorded time series. Estimation of the degree of a
polynomial potential allows us to derive the number of potential
wells in a system. The method correctly detects changes in the number of wells
in artificial data. In ice-core proxy records of Greenland
paleotemperature, a reduction in the number of climate states from
two to one is
detected sometime prior to the last glacial maximum (LGM),
23--19\,kyr\,BP.
This result is also found in analysis of Greenland Ca data.
The bifurcation can be interpreted as loss of
stability of the warm interstadial state of the Dansgaard-Oeschger~(DO)
events. The proposed method can be applied to a wide
range of geophysical time series exhibiting bifurcations.
\end{abstract}

\introduction
Knowing the number of states in a geophysical system is
crucial for understanding its underlying dynamics and for modeling
its behavior.  Of particular interest are changes in the number of states which
represent bifurcations of a system.
Commonly the number of system states is inferred
either by simply looking at time series data, or on the basis of
theoretical models. Yet when there are two or more states present in
a system the precise number is not always obvious from examining
time series by eye (see e.g.~Fig.~\ref{combined}b). In particular,
what looks like a system with two states may actually contain more states.
Consequently we have developed a more formal (though still
necessarily approximate) method of estimating the number of system
states from time series data.

In studies of paleoclimate changes, it is common to describe the
glacial-interglacial cycles as shifts between different climate
states. For example, one simple model \citep{paillard} assumes three
(mild glacial, full glacial, interglacial)
states for Northern Hemisphere ice volume.
Within the last ice age, the Dansgaard-Oeschger events
that are recorded in Greenland ice cores \citep{dansgaard}, North
Atlantic palaeoclimate records \citep{bond}, and are now recognized
further afield \citep{wang}, are commonly treated as a two-state
system (cold stadial, warm interstadial). The DO events are widely
explained in terms of a model \citep{stommel,rahmstorf1,ganopolski}
for bistability of the Atlantic ocean thermohaline circulation,
coupled to changes in sea-ice cover and atmospheric circulation.
Potentially, different combinations of Northern Hemisphere ice
volume and ocean-atmosphere-sea-ice state could give rise to
multiple climate states as recorded in Greenland. Ice volume clearly
affects the stability regime of the ocean-atmosphere-sea-ice as the
Holocene interglacial has not exhibited the extreme DO variability
seen in the last ice age.

Here we apply our new method to examine how many climate states are
present in high-resolution Greenland Ice-Core Project (GRIP)
\citep{dansgaard} and North Greenland Ice-Core Project (NGRIP)
\citep{ngrip} records, as a function of time through the last 60\,kyr
(part of the last ice age and the Holocene). 
We also consider GRIP calcium data, which has annual resolution, 
and although it is available since 91~kyr~BP, we use data from 60~kyr~BP
to make comparison on the same time interval with the above two temperature proxies.

The precise mechanisms
causing transitions between climate states do not concern us here.
We make the simplest
assumption that multiple climate states can be approximated by a
non-oscillatory polynomial potential, albeit one that changes
through time. Furthermore, we assume that transitions between
climate states are triggered purely by stochastic noise. This is
consistent with recent analysis of the DO events (especially in the
GRIP and NGRIP records), which does not support any underlying
periodicity~\citep{ditlevsen}.

\section{Dynamical framework}

We treat the climate system as a nonlinear
dynamical system which can possess multiple states, with shifts
between these different climate states induced by stochastic
forcing. A one-dimensional conceptual model for this is given by the
Langevin equation
\begin{equation}
\label{eq_pot}
dz=-U^\prime (z)\,dt + \sigma dW,
\end{equation}
where $U(z)$ is a potential function, $\sigma$ is the noise level and 
$W$ denotes a Wiener process.
The state variable
$z$ represents a large-scale climate variable like Greenland %Northern Hemisphere
palaeotemperature and is here identified with ice-core proxy
($\delta^{18}$O stable water isotope) records.
In~\citep{ditlevsen99}, the Fokker-Planck equation was used to derive the potential
of the  GRIP calcium record, with white non-Gaussian noise and two noise terms in the Langevin
equation.

The shape of the potential is given by a polynomial
\begin{equation}
U(z)=\sum_{i=1}^L a_i z^i,
\end{equation}
where the order $L$ is even and the leading coefficient $a_L$ is
positive for Eq.~(\ref{eq_pot}) to possess a stationary solution.
The order of the polynomial controls the complexity of the
potential. Increasing values of $L$ allow more states to be
accommodated; for example, a fourth-order polynomial can capture a
system with two states (double-well potential)
\citep{kwasniok1,kwasniok2}.

The number of system states is estimated by means of a polynomial
fit of the probability density function of the data.
Suppose the system is governed by Eq.~(\ref{eq_pot}). The corresponding
Fokker-Planck equation for the probability density function $p(z,t)$
\begin{equation}
\partial_t p(z,t)=\partial_z [U'(z)p(z,t)]
+\frac{1}{2}\sigma^2\partial_z^2 p(z,t)
\end{equation}
has a stationary solution given by \citep{gardiner}
\begin{equation}
p(z)\sim \exp[-2U(z)/\sigma^2].
\end{equation}
Given this one-to-one correspondence between the potential and the stationary
probability density of the system, the potential can be reconstructed from
time series data of the system as
\begin{equation}
U=-\frac{\sigma^2}{2}\log p_{\rm d},
\end{equation}
where $p_{\rm d}$ is the empirical probability density of the data.
This is estimated using
a standard Gaussian kernel estimator~\citep{silverman}.
We applied Matlab code with built-in {\tt kdensity} function and the
NAG Fortran library (computing the Gaussian kernel density estimator using 
a fast Fourier transform).
The estimator is
\begin{equation}
\hat f(z) = \frac 1 {nh} \sum\limits_{i=1}^n K\left(\frac {z-z_i} h\right),
\end{equation}
where $K$ denotes the Gaussian kernel, $n$ is the length of the data set and 
$h$ is the bandwidth controlling the smoothness of the estimator.
Following~\citep{silverman}, we chose $h=1.06 s/n^{1/5}$, where $s$ is 
the standard deviation of the data set.
Then least-square fits of $-\log p_{\rm d}$ weighted with the probability
density of the data \citep{kwasniok1} with
polynomials of increasing even order $L$ are
calculated, starting with $L{=}2$, until a negative leading coefficient $a_L$
is encountered. The polynomial of highest degree before first obtaining a
negative leading coefficient is considered the most appropriate representation
of the probability density of the time series, both locally and globally,
avoiding overfitting of sampling fluctuations in the probability density.

The number of states $S$ in the system is then determined as
\begin{equation}
\label{inflections}
S= 1 + \frac{I}{2},
\end{equation}
where $I$ is the number of inflection points of the fitted polynomial
potential of appropriate degree $L$ as described above. This definition takes
into account not only the degree of the polynomial but its actual shape. We
only look at even-order potentials with positive leading coefficient. These
have positive curvature both at minus and plus infinity. Thus, inflection
points can only occur in pairs (if any). Any potential has at least one state
(with no inflection points). Then we count one further state for each pair
of inflection points. This can be either a real minimum (well) or just a
flattening in the potential corresponding to a degeneracy in the potential;
definition (\ref{inflections}) accommodates
both possibilities. The number of inflection points is numerically given as
the number of sign changes in the second derivative on a fine mesh.

Examples of the polynomial fits for simulated data are given in
Fig.~\ref{pots}. Given
double-well potential data ($U(z){=}z^4{-}2z^2$),
the probability density of the time series
has bimodal shape. The empirical potential given as $U{=}{-}\frac{\sigma^2}{2}
\log p_{\rm d}$ can
be fitted by polynomials of various orders, but only those of degree
2 and 4 have a positive leading coefficient. Higher-order polynomials have
negative leading
coefficients --- although fitting the data better and better locally,
they are clearly inadequate globally. The fourth-order polynomial is
the appropriate representation of the  probability density and the method
identifies the time series as double-well potential.

A simplified version of definition (\ref{inflections}) would be to only count
the relative minima (real wells) in the potential. We opt here for the more
comprehensive definition of system states as degenerate or nearly degenerate
potentials have been
shown to occur in the context of ice-core records \citep{kwasniok1,kwasniok2}.
Moreover, even if the system possesses real potential wells (for example, the
artificial data shown in Fig.~2) changes in the number of states (bifurcations)
tend to be picked up more quickly when using the inflection point definition
as this is a more sensitive feature.

%For all algorithms we have developed Matlab routines that can be
%provided upon enquiry by the corresponding author.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Data}

To test the method for estimating the number of system states, we
generated four sets of artificial data for four different
potentials: one-well ($U(z)=z^2$), double-well
($U(z)=z^4-2z^2$), triple-well ($U(z)=z^6-4.5z^4+5z^2$)
and four-well potential ($U(z)=z^8-6.5z^6+13z^4-8z^2$)
(Fig.~\ref{combined}a). Simulated data of these systems were
obtained by numerically integrating eq.~(\ref{eq_pot}) using the
Euler scheme. The noise level $\sigma=1.5$ was adopted, resulting in a mean 
state transition time (Kramers waiting time) similar to that occurring in 
the palaeoclimatic records when interpreting system time units as kyr. The sampling
interval is 0.05 system time units. For each of the potentials, 150 time
units worth of data were generated, resulting in 3000 data points.
The subsets are combined to produce a record with 12\,000 data
points, corresponding to time running from 0 to 600
(Fig.~\ref{combined}b).

Having tested the methods on artificial data, we examined ice-core
proxy records of palaeotemperature at two different sites 325\,km
apart in Greenland; GRIP~\citep{dansgaard} (Fig.~\ref{NGRIP}a) and
NGRIP~\citep{ngrip} (Fig.~\ref{NGRIP}c), on the most recent GICC05
time scale \citep{svensson2}.
The data are $\delta^{18}$O stable water isotope records, which are
a proxy for past air temperature at the ice-core sites. The records
can also be influenced by changing water source temperatures and
snowfall seasonality. There are significant differences between the
two records during the last ice age, with NGRIP systematically
colder (depleted in $^{18}$O), and more variable, when Northern
Hemisphere ice sheets are more extensive~\citep{ngrip}. This may be
because NGRIP received a greater fraction of colder air coming over
the northern side of the Laurentide ice sheet.
Importantly, the difference between the two records does not show a
component of millennial variability. Here
we consider the last 60 kyr NGRIP series with resolution 20\,yr and a
part of the GRIP series (which is 112 kyr long) in the same interval
(60--0\,kyr\,BP) with the same temporal resolution of 20\,yr. For
the detection of changing numbers of states we use sliding windows
of varying length through each dataset.

To further test initial results obtained on the $\delta^{18}O$ data, 
we also examined Ca data from the GRIP project (Fig.~\ref{calcium}) that has 
annual resolution and spans the interval 91--11~kyr~BP. In our analysis, we consider 
a subset of this dataset starting at 60~kyr~BP to make it comparable with the other two datasets.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Results}

We visualize the estimated number of system states as a colour
contour plot, expressed as a function of the time at the middle of the
data window ($x$-axis, aligned with the data time scale) and the time
length of the data window ($y$-axis). As the results are plotted at
the middle of their corresponding time windows, each point in the
contour plot should be compared with the segment of the time series
data centered at that particular point, and having the length of the
corresponding time window.

In the artificial data, the number of states is generally identified correctly; 
changes are picked up quickly and reliably (Fig.~\ref{combined}c). For the smallest time windows, 
there are sporadic misidentifications of
the number of states, due to poor statistics.
Generally time windows of order 400 data points (20 time units) are
sufficient to get reliable results. 

The GRIP (Fig.~\ref{NGRIP}a) and NGRIP (Fig.~\ref{NGRIP}c)
palaeotemperature data are highly variable and non-stationary.
Nevertheless, the algorithm detects a number of interesting common
features in both records (Fig.~\ref{NGRIP}b,d). Over
60--25\,kyr\,BP, two climate states are most commonly detected,
consistent with the conventional view of the Dansgaard-Oeschger  events.
Exceptions of 1-state detection appear to be caused when long
intervals of steady cooling, e.g.~circa 55--50\,kyr\,BP, occupy a
significant part of the analysis window.

The most pronounced feature, common to both records, is a change
from 2-well to 1-well potential, generally detected by 25~kyr~BP
and inferred to
have occurred somewhat prior to its detection. The transition is most sharply
detected in the GRIP data.
In the NGRIP data, the transition is less pronounced, consistent
with the greater noise level in this dataset obscuring the detection
of states.
In both records, the shift to a 1-well potential is persistent,
indicating a bifurcation in the climate system that occurred late in
the last ice age but prior to the Last Glacial Maximum (LGM)
23--19\,kyr\,BP \citep{yokoyama}.

The return of a second state in the climate is only detected in both
records, for most window sizes, around 12\,kyr\,BP (i.e.~at the end
of the Younger Dryas), although it is sporadically detected earlier,
from the time of the Bolling warming, especially in the NGRIP data.
This reflects the growing influences of changes during the
deglaciation on the algorithm. The detection of up to four states
as the data window spans the interval of the last
glacial termination is probably an artefact of the nonstationarity in the time series.
Although these four climate states plausibly represent the LGM, Bolling-Allerod warm
interval, Younger Dryas cold interval, and the Holocene, it is perhaps more accurate to interpret the system
as having two states about a moving trend.

As the end of the analysis window advances into the Holocene, the
estimated number of system states declines. However, the number of
climate states only reduces to one for both datasets when the analysis
window comprises only the data since around 10\,kyr\,BP, i.e.~the
start of the Holocene.
For the GRIP data, there is some indication that the 8.2\,kyr\,BP
event can influence the algorithm sufficiently that two climate states
are detected when using small windows that span the event. The
8.2\,kyr\,BP event is also visible in the NGRIP data but has poorer
signal-to-noise ratio, so does not exert the same influence on the
algorithm.

Given that the ice-core time series are much shorter (in terms of the
number of recorded transitions between states) than the simulated time
series used above, the question arises how robust
our results on the ice-core data are statistically. The mean recurrence
time (Kramers waiting time) between DO events is 2.8\,kyr which is about
the same as in the artificial data in nondimensional system time units. Hence
the maximum window size of 20\,kyr in Fig.~\ref{NGRIP}b,d corresponds to 20 units
on the vertical axis in Fig.~\ref{combined}c. This is just enough to enable a correct
detection of the number of states with very high probability. Together
with the fact that our results are robust for a range of window sizes
between 10 and 20\,kyr this gives us some confidence that our conclusions
are reliable.

As a further test, we examined GRIP Ca data, and the obtained result (Fig.~4) confirms the detected 
bifurcation in GRIP and NGRIP $\delta^{18}O$, although the datasets have different origin and 
temporal resolution.

\conclusions

The results support the widely held view of the Dansgaard-Oeschger
events that they represent switches between 2 states of the climate
system, which were present from before the start of our analysis
60~kyr~BP.  Most interesting is the detection of a bifurcation in
the climate system that reduced the number of states from 2 to 1,
sometime prior to the last glacial maximum
(LGM) 23--19\,kyr\,BP. This can be interpreted as the loss of a stable
warm interstadial state.
This result is not inconsistent with the conventional labelling of
DO warming event 2 occurring at 23.4\,kyr\,BP timing according
to~\citep{rahmstorf3} or the earlier DO events 3 at 27.8\,kyr\,BP and
4 at 29.0\,kyr\,BP. Rather we suggest that although transitions to a
warm state were still triggered up to the LGM, the resulting state
became not even marginally stable. In the case of DO event 2 at
least (and possibly 3 and even 4), the system was always destined to
revert back to a cold state at a rate determined by its internal
dynamics and/or the removal of some forcing factor. Indeed
shortening of the duration of the warm interstadial state
can be seen in the time series (Fig.~\ref{NGRIP}a,c)
as time progresses toward the LGM.

Clearly a warm state reappears during the deglaciation, in the form
of the Bolling warming event (conventionally labelled DO event 1 at
14.6\,kyr\,BP~\citep{rahmstorf3}), but whether it is appropriate to
liken this to the earlier interstadials or the Holocene is unclear
from our analysis. Likewise whether the Younger Dryas cool event
should be likened to the earlier cold stadial state is ambiguous.
Clearly Northern Hemisphere ice volume had declined significantly
prior to either of these events, and was undergoing a transition
from glacial to interglacial states.
For the last 10\,kyr of the Holocene, under interglacial ice
volume only 1 stable Atlantic ocean-atmosphere-sea-ice
state is detected. We interpret the 8.2\,kyr\,BP cooling event as a meltwater
induced transient change in ocean-atmosphere-sea-ice from which the
system was always destined to recover.

In this view the stability regime
of the Atlantic ocean-atmosphere-sea-ice system is reshaped and
sometimes bifurcated by longer timescale changes in ice volume,
as suggested by a number of models,
e.g.~\citep{ganopolski,colindeverdiere}.
Bistability of the Atlantic ocean-atmosphere-sea-ice system is
present for intermediate Northern Hemisphere ice volumes, whereas
under glacial maximum conditions there is only one stable cold state
for Greenland temperature, and under interglacial conditions there
is only a single (different) stable warm state.

The Langevin equation (eq.~(1)) used in the present paper as a conceptual model is a 
fairly natural starting point which is already quite general but admittedly has limitations. 
The model is of first order; a higher-order model could accommodate a richer dynamics. 
Sometimes, palaeoclimatic cycles are treated as oscillations, e.g.~\citep{ghil}. 
Our method can be readily extended to a second-order stochastically driven potential 
model supporting oscillatory behaviour \citep{kwasniok2}. Moreover, eq.~(1) does not 
allow for non-random external forcing, e.g.~orbital forcing. Also the possibility of 
state-dependent noise rather than purely additive noise as well as coloured noise would be 
interesting to include. The proposed approach is initial, and further generalization 
of the method is to be developed.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{acknowledgements}
We would like to thank J.~Kantelhardt and Yu.~Sapronov for
discussions and P.~Ditlevsen for useful comments on the manuscript.
We are grateful to P.~Ditlevsen for providing the GRIP Ca dataset for our analysis.
The research was supported by NERC through project NE/F005474/1.
\end{acknowledgements}

%%%%%%%%%%%%%%%%%%%%%%
%\providecommand{\natexlab}[1]{#1} \providecommand{\url}[1]{{\tt #1}}
%\providecommand{\urlprefix}{URL } \expandafter\ifx\csname
%urlstyle\endcsname\relax
%\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
%\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
%\urlstyle{rm}\Url}\fi
%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{}

\bibitem[Bond et al.(1997)]{bond}
Bond, G. Showers, W., Cheseby, M., Lotti, R., Almasi, P., de Menocal, P.,
Priore, P., Cullen, H., Hajdas, I., and Bonani, G.: A Pervasive Millennial-Scale Cycle in
North Atlantic Holocene and Glacial Climates, Science,
278, 1257--1266, 1997.

\bibitem[Colin de Verdiere et al.(2006)]{colindeverdiere}
Colin de Verdiere, A., Ben Jelloul, M., and Sevellec, F.: Bifurcation
structure of thermohaline millennial oscillations,  J. Climate, 19, 5777--5795, 2006.

\bibitem[Dansgaard et al.(1993)]{dansgaard}
Dansgaard, W., Johnson, S. J., Clausen, H. B., Dahljensen, D., Gundestrup, N. S.,
Hammer, C. U., Hvidberg, C. S., Steffensen, J. P., Sveinbjornsdottir, A. E.,
Jouzel, J., and Bond, G.: Evidence for general instability of
past climate from a 250\,kyr ice-core record,  Nature,
364, 218--220, 1993.

\bibitem[Ditlevsen(1999)]{ditlevsen99}
Ditlevsen, P.~D.: Observation of $\alpha$-stable noise induced millennial climate changes from an ice-core
record, Geophys.~Res.~Lett., 26(10), 1441--1444, 1999.

\bibitem[Ditlevsen et al.(2007)]{ditlevsen}
Ditlevsen, P.~D., Andersen, K.~K., and Svensson, A.: The DO-climate
events are probably noise induced: statistical investigation of the
claimed 1470 years cycle, Climate of the Past, 3(1), 129--134, 2007.

\bibitem[Ganopolski and Rahmstorf(2001)]{ganopolski}
Ganopolski, A. and Rahmstorf, S.: Rapid changes of glacial climate
simulated in a coupled climate model, Nature, 409, 153--158, 2001.

\bibitem[Gardiner(2004)]{gardiner}
Gardiner, C. W.: Handbook of Stochastic Methods for Physics, Chemistry  and
the Natural Sciences (Springer Series in Synergetics), Springer, 3rd ed.,
415~pp., 2004.

\bibitem[Ghil et al.(1987)]{ghil}
Ghil, M., Mullhaupt, A., and Pestiaux, P.:
Deep water formation and Quaternary glaciations, Climate Dyn., 2, 1--10, 1987.

\bibitem[Kwasniok and Lohmann(2009a)]{kwasniok1}
Kwasniok, F. and Lohmann, G.:
Deriving dynamical models from paleoclimatic records: Application to glacial
millennial-scale climate transitions,
Phys.~Rev.~E, 80, 066104, 2009.

\bibitem[Kwasniok and Lohmann(2009b)]{kwasniok2}
Kwasniok, F. and Lohmann, G.:
A stochastic nonlinear oscillator model for glacial millennial-scale climate
variability derived from ice-core data,
Nonl.~Proc.~Geophys., submitted, 2009b.

\bibitem[NGRIP project(2004)]{ngrip}
North Greenland Ice Core Project members:
High-resolution record of Northern Hemisphere climate extending
into the last interglacial period, Nature, 431, 147--151, 2004.
Data are available at:
\url{http://www.iceandclimate.nbi.ku.dk/data/}.

\bibitem[Paillard(2001)]{paillard}
Paillard, D.:
Glacial cycles: toward a new paradigm,
 Reviews of Geophys., 39(3), 325--346, 2001.

\bibitem[Rahmstorf(1995)]{rahmstorf1}
Rahmstorf, S.: Bifurcation of the Atlantic thermohaline
circulation in response to changes in the hydrological cycle,
Nature, 378, 145--149, 1995.

\bibitem[Rahmstorf(2003)]{rahmstorf3}
Rahmstorf, S.: Timing of abrupt climate change: a precise clock,
 Geophys.~Res.~Lett., 30(10), 1510, 2003.

\bibitem[Silverman(1986)]{silverman}
Silverman, B. W.: Density estimation for statistics and data analysis,
Chapman and Hall, 1986.

\bibitem[{Stommel}(1961)]{stommel}
Stommel, H.: Thermohaline convection with two stable regimes of flow,
Tellus, XIII, 224--230, 1961.

\bibitem[Svensson et al.(2008)]{svensson2}
Svensson, A., Andersen, K. K., Bigler, M., Clausen, H. B., Dahl-Jensen, D.,
Davies, S. M., Johnsen, S. J., Muscheler, R., Parrenin, F., Rasmussen, S. O.,
Roethlisberger, R., Seierstad, I, Steffensen, J. P., and Vinther, B. M.:
A 60000 year Greenland stratigraphic ice core chronology,
Climate of the Past, 4, 47--57, 2008.

\bibitem[Wang et al.(2001)]{wang}
Wang, Y. J., Cheng, H., Edwards, R. L., An, Z. S., Wu, J.Y., Shen, C. C., and Dorale, J. A.: A High-Resolution Absolute-Dated Late
Pleistocene Monsoon Record from Hulu Cave, China, Science,
294, 2345--2348, 2001.

\bibitem[Yokoyama et al.(2000)]{yokoyama}
Yokoyama Y., Lambeck, K., DeDekker, P., Johnston, P., and Fifield, L.~K.: Timing of the Last Glacial Maximum from
observed sea-level minima,  Nature, 406, 713--716, 2000.

\end{thebibliography}

%%%%%%%%%%%%%%%%%%%% FIGURES
\begin{figure}[h!]
\includegraphics[width=0.45\textwidth]{fig1.eps}
\caption{Empirical potential (magenta curve) of double-well
potential artifical data and its polynomial fits of 2nd, 4th and 6th
orders. The parabolic fit (2nd order) is not as accurate as the
4th-order polynomial, whereas the fit of 6th order approximates the
potential only locally, becoming negative for large $|z|$.
The numerical algorithm chooses the highest
order of polynomial that approximates the potential globally with
positive leading coefficient; in this case it is the 4th-order
polynomial. The method correctly classifies the data as double-well
potential.}
\label{pots}
\end{figure}

\begin{figure}[h!]
\includegraphics[width=0.45\textwidth]{fig2.eps}
\caption{Detection of the number of states in artificial data:
a)~four different potentials with increasing number of wells from
one to four; b)~the resulting time series; c)~contour plot of the
number of detected states versus time and size of sliding window
(red -- one, green -- two, cyan -- three, purple --
four). The number of detected states is mapped at the middle of the
sliding time windows.}
\label{combined}
\end{figure}

%%%%\begin{figure*}[p]
\begin{figure}[h!]
%\begin{landscapefigure*}
\includegraphics[width=0.45\textwidth]{fig3.eps}
\caption{GRIP ice-core data: a)~time series; b)~contour plot of the
number of detected states versus time and size of sliding window,
mapped at the middle of the sliding time windows (red -- one,
green -- two, cyan -- three, purple -- four).
NGRIP ice-core data: c)~time series; d)~contour plot of the
number of detected states versus time and size of sliding window,
mapped at the middle of the sliding time windows (red -- one,
green -- two, cyan -- three, purple -- four).}
\label{NGRIP}
%\end{landscapefigure*}
\end{figure}
%%%%%%%%%%%%5\end{figure*}

%%%%\begin{figure*}[p]
\begin{figure}[h!]
%\begin{landscapefigure*}
\includegraphics[width=0.45\textwidth]{fig4.eps}
\caption{GRIP calcium data: a)~time series; b)~contour plot of the
number of detected states versus time and size of sliding window,
mapped at the middle of the sliding time windows (red -- one,
green -- two, cyan -- three, purple -- four).}
\label{calcium}
%\end{landscapefigure*}
\end{figure}
%%%%%%%%%%%%5\end{figure*}


\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CP FIGURES PROCESSING

\clearpage

%f01
\begin{figure}
\vspace*{2mm}
\center\includegraphics[width=12cm]{cp-2009-53-discussions-f01.pdf}
\caption{Empirical potential (magenta curve) of double-well
potential artifical data and its polynomial fits of 2nd, 4th and 6th
orders. The parabolic fit (2nd order) is not as accurate as the
4th-order polynomial, whereas the fit of 6th order approximates the
potential only locally, becoming negative for large $|z|$.
The numerical algorithm chooses the highest
order of polynomial that approximates the potential globally with
positive leading coefficient; in this case it is the 4th-order
polynomial. The method correctly classifies the data as double-well
potential.
\label{pots}}
\end{figure}

\clearpage

%f02
\begin{figure}
\vspace*{2mm}
\center\includegraphics[width=9cm]{cp-2009-53-discussions-f02.pdf}
\caption{Detection of the number of states in artificial data:
{\bf (a)} four different potentials with increasing number of wells from
one to four; {\bf (b)} the resulting time series; {\bf (c)} contour plot of the
number of detected states versus time and size of sliding window
(red -- one well, green -- two, cyan -- three, purple --
four). The number of detected wells is mapped at the end of the
sliding time windows.
\label{combined}}
\end{figure}

\clearpage

%f03
\begin{figure}
\vspace*{2mm}
\center\includegraphics[width=7cm]{cp-2009-53-discussions-f03.pdf}
\caption{Potential analysis of Greenland ice-core data: {\bf (a)} GRIP time series; {\bf (b)} contour plot of the
number of detected states versus time and size of sliding window,
mapped at the end of the sliding time windows (red -- one well,
green -- two, cyan -- three, purple -- four).
{\bf (c)} NGRIP time series; {\bf (d)} contour plot of the
number of detected states versus time and size of sliding window,
mapped at the end of the sliding time windows, as in (b).
\label{NGRIP}}
\end{figure}

\end{document}

