Evaluation of empirical mode decomposition for quantifying multi-decadal variations in sea level records

Introduction Conclusions References Tables Figures


Introduction
Over the last decade, several papers have used the method of empirical mode decomposition (EMD) (Huang et al., 1998;Huang and Wu, 2008) to evaluate non-stationary patterns in time series as disparate as electromyographic signals (Andrade et al., 2006) and sea level (Breaker and Ruzmaikin, 2011;Ezer and Corlett, 2012;Ezer et al., 2013;Lee, 2013;Chen et al., 2014).The use of EMD in sea level records has been motivated in large part by numerous papers discussing the appearance of decadal and longer period fluctuations in tide gauge records and global mean sea level estimates based on tide gauge records (e.g., Feng et al., 2004;Miller and Douglas, 2007;Woodworth et al., 2009;Bromirski et al., 2011;Sturges and Douglas, 2011;Chambers et al., 2012;Calafat and Chambers, 2013;Becker et al., 2014;Dangerdorf et al., 2014).At first glance, EMD appears to be a useful tool to find non-stationary, low-frequency fluctuations in sea level, as it breaks the time-series into a set of Intrinsic Mode Functions (IMFs) that have progressively longer quasi-period fluctuations.IMFs extracted from various tide gauge records have been shown to be correlated with several climate indices (e.g., Ezer and Corlett, 2012;Ezer et al., 2013), which gives some credence to extracted signals.Moreover, authors have argued that the final IMF, representing the continuously increasing sea level mode, is a better representation of an acceleration than simply fitting a quadratic to the original data (Huang and Wu, 2008;Ezer and Corlett, 2012;Ezer et al., 2013).
However, there are some potential pitfalls that we believe have not been fully addressed in previous papers utilizing the method.First and foremost, EMD is a purely mathematical deconstruction of the data, with no regard to intrinsic covariance of the signals or physics.Second, it assumes that IMFs are comprised of fluctuating signals where the magnitude of nearby peaks and troughs are balanced to create a zero mean -an assumption not based on any physical requirement, as real observations can have quite large ranges in magnitudes, especially sea level data affected by climate signals and synoptic storm events.
The basic idea of EMD is to fit cubic splines to the local maxima and minima of a timeseries separately, average the splines, then remove the average from the time-series.
The process is iterated on the residual time-series until the average of the splines converges to have a SD less than some pre-set tolerance.This is the first IMF.This is then subtracted from the original time-series and the process is repeated until only one minimum and one maximum remain.For details of the procedure, readers are referred to the original paper by Huang et al. (1998) or more recent applications (e.g., Huang and Wu, 2008;Ezer et al., 2013).
There is also a subtlety in finding the last IMF that is not discussed in the literature.Since the EMD process requires fitting of cubic splines, the last IMF mode that can Introduction

Conclusions References
Tables Figures

Back Close
Full be calculated has more than one local minima and more than one local maxima, but fewer than four.The only way to get the final IMF shown in most studies, which shows a continuously increasing sea level mode, is to fit a quadratic to the final IMF from the EMD process, and plot the resulting fit.This is conceptually no different than fitting a quadratic to the original time series, although the variance of the final EMD mode is far lower than the original data, which should improve the estimate -provided there are no systematic errors or biases in the final IMF that would bias the result.Moreover, Wu and Huang (2004) have shown that EMD behaves as a low-pass filter.If one runs random noise with a normal Gaussian distribution through the process, low-frequency signals will be seen in the resulting IMFs.They found there is roughly a doubling of the average period with each subsequent mode.In a simple example of running EMD on a monthly-resolution time-series that is 150 years long with random noise that has a SD of 60 mm, one clearly observes IMFs that have quasi-period fluctuations of nearly 60 years and amplitudes as large as 10 mm (Fig. 1); fluctuations at quasi-30 year periods are the same magnitude.Note that we used 60 mm because this is the SD of residual monthly sea level at San Francisco after fitting and removing a quadratic function plus annual and semiannual sinusoids, so is representative of high-frequency sea level at a typical site, although some sites can have significantly higher variability.Wu and Huang (2004) and others (Ezer and Corlett, 2012;Ezer et al., 2013) have attempted to account for this behavior in the error statistics of their recovered IMFs.This is done by creating several pseudo-time-series by adding a small amount of random noise to the original time-series, running a large number of EMDs, and considering the average.The SD of the differences represents the uncertainty.However, we see a problem with this method.Although adding random noise will account for uncertainty in the IMFs from uncertainty in the measurement, it does not consider the effect of the real signal in the tide gauge data that has high variance and little serial correlation.This "signal" will also be filtered by the EMD process and will likely appear as Introduction

Conclusions References
Tables Figures

Back Close
Full a quasi-stationary oscillation in higher order IMFs that is not real (i.e., similar to the EMD filtering of purely random noise shown in Fig. 1).Thus, we propose to test the EMD process not on real observations where one does not know the underlying modes, but on three simulated data sets where the modes are prescribed.We believe the differences between the recovered IMFs and given signals will be a better measure of the accuracy of the EMD method than what has previously been discussed in the literature.Two different simulations will be examined with fluctuating signals and differing random noise to represent high-frequency variability: one using purely sinusoidal oscillations with multiple periods ranging from 13 to 80 years, the second with variations based on band-pass filtered and scaled El Niño-Southern Oscillation (ENSO) and Pacific Decadal Oscillation (PDO) indices, both with additional random noise applied.The third case will examine a situation with only seasonal fluctuations, random noise, and a single month with a variation larger than 3 SDs.This represents an extreme event, typically caused by major storm surge, which is a common feature in many sea level records.We will demonstrate that the EMD method leads to spurious IMFs in all cases, and where the IMFs are correlated with the input signal, their amplitudes and phases are significantly biased in many periods of the record.

Data and methods
We use the monthly sea level record from the San Francisco tide gauge for our reference sea level record.It was downloaded from the Permanent Service for Mean Sea Level (PSMSL) (Woodworth and Player, 2003;PSMSL, 2012).We fit annual and semiannual sinusoids along with a trend and an acceleration term to the data using ordinary least squares to obtain the base-line sea level variability for our model (y base ), where

Conclusions References
Tables Figures

Back Close
Full For Case 1, we add three long-period sinusoids (13, 55, and 80 years) to the baseline model along with random noise (ε(t)): (2) The random noise has a variance to match the variance of the residuals of the real tide gauge data minus the model.1000 different random noise models were applied to create 1000 different versions of Case 1 for statistical testing.The periods and amplitudes of the long-period sinusoids were chosen arbitrarily to approximate the level of multidecadal fluctuations in the San Francisco sea level record (Fig. 2a).
Case 2 starts from the same baseline model, but instead of prescribing sinusoidal oscillations, we use non-stationary climate indices for El Niño-Southern Oscillation (ENSO) variations and the Pacific Decadal Oscillation (PDO).We use the Southern Oscillation Index based on the pressure differences between Tahiti and Darwin, Australia to represent ENSO variability (Trenberth, 1984;Ropelewski and Jones, 1987; downloaded from http://www.cgd.ucar.edu/cas/catalog/climind/soi.html on 5 March 2014), and the PDO index based on the leading principal component of sea surface temperature in the North Pacific (Zhang et al., 1997;Mantua et al., 1997; downloaded from http://jisao.washington.edu/pdo/PDO.lateston 5 March 2014).
Several additional processing steps are required before using these indices for our experiment.First, neither index covers the same period as the tide gauge (January 1856 to December 2010).The SOI index starts in January 1866 while the PDO index begins in January 1900.In order to have a simulated record as long as possible, we start in January 1866 and use values from the end of PDO record to fill in the missing data before January 1900.Recall we are not worried about the "true" ENSO Introduction

Conclusions References
Tables Figures

Back Close
Full or PDO variability, only a simulation of the type of variability and how well EMD can recover it.Secondly, because the two indices are correlated due to similar interannual (< ten year) variations, we low-pass the PDO index with a 5 year Gaussian, and band-pass filter the SOI by first removing the 5 year Gaussian of the SOI, and then filtering the residuals with a 0.5 year Gaussian.After doing this, the correlation between the two filtered indices PDO LP (t) and SOI BP (t) is insignificant (−0.003).
The final step is to determine the scaling factor to apply to both the PDO LP (t) and SOI BP (t) variations.This is done by first normalizing both time-series by their SD.Then, after removing the estimated trend, acceleration and seasonal variations, the sea level data are low-pass and band-pass filtered as the climate indices were, and the SD of the filtered residuals is calculated.The scaling factor applied to PDO LP (t) is the SD of the low-pass filtered sea level residuals (20.8 mm); the scaling factor applied to SOI BP (t) is the SD of the band-pass filtered sea level residuals scaled by −1 to account for the fact El Niño sea level variations at San Francisco are positive when the SOI is negative (−28.7 mm).
The final time-series for Case 2 is assembled as for Case 1, including the random noise term based on the SD of the residuals and the model: noting PDO LP (t) and SOI BP (t) are normalized as described previously.One time-series is shown in Fig. 2b to show the model does a reasonable job of simulating the San Francisco tide gauge record.Case 3 starts as the baseline model, adds random noise with a SD of 60 mm (representative of the high-frequency variability in San Francisco sea level), then adds an extra 350 mm for January 1956 to represent the signal of a large storm surge event.Such a value is possible in sea level records, depending on the size and duration of the storm (e.g., the maximum deviation of monthly sea level residuals after removing a trend for the San Francisco tide gauge record is 4.9 times higher than the SD).Moreover, most Introduction

Conclusions References
Tables Figures

Back Close
Full tide gauge records have numerous events instead of just one; San Francisco has six monthly residuals exceeding 200 mm and two exceeding 300 mm.For this study, however, we consider just one to demonstrate the effect on EMD if authors do not consider this possibility in their analysis.We perform EMD using the EMD Toolkit for SciLab (http://www.scilab.org),based on code documented in Rilling et al. (2003).The specific function utilized was emdc, which stops iterating when a tolerance is reached.We used a tolerance value of 0.05.

Results and analysis
Figure 3 shows the low-frequency IMFs for a single simulation of Case 1, along with the input signals.IMFs 1-5 are all much higher frequency and so are not considered, although we will note that none accurately captures the input seasonal variation.However, we point out that some of the artifacts shown in Fig. 3 for the low-frequency IMFs are a direct result of correcting for errors in the higher frequency IMFs not shown, so that the sum of all matches the original data.
The correlation of IMF6 with the prescribed 13 year sinusoid is significant (> 0.5), but not high.It is clear there are several periods where the EMD method would suggest no variability at a 13 year period  and other periods (∼ 1910) where the variation is significantly faster.Moreover, the amplitude of the recovered IMF is steadily increasing after 1980, although the phase is about correct.The next IMF is an artifact of the method, with no significant correlation with any input signal, yet showing a periodicity of ∼ 20 years with amplitudes as high as 20 mm.The longer period IMFs also have problems (Fig. 3).The one best correlated with the 55 year sinusoid (IMF8) is out of phase with the real signal until about 1940, and the amplitude is increasing in time.The 80 year IMF exhibits a similar behavior of increasing amplitude (Fig. 3).Finally, the estimated quadratic term to the longest oscillatory IMF (IMF9 in this case), significantly underestimates the prescribed acceleration.Introduction

Conclusions References
Tables Figures

Back Close
Full This is just a single example, however.To see if averaging of multiple versions of random noise will better match the prescribed signals, we examine the 1000 different IMF clusters from the simulations.One cannot simply separate the corresponding lowfrequency modes based on the IMF number, however, as the total number of IMFs changed from 9 to 11 in the 1000 different simulations.We found that the 13 year signal was found in IMFs ranging from number 4 to 7, while the 55 year mode ranged from IMF 7-9.The last mode found ranged from IMF 9-11.Thus, we had to rely on correlation with the known oscillation to identify the relevant IMF.This was done by computing the correlation of each recovered IMF from each simulation with the prescribed sinusoids.A minimum bound was set to 0.5.If two or more modes had correlations > 0.5 with one of the input signals, the one with the highest correlation was chosen.
Figure 4 summarizes the results, showing the mean with the SD as a shaded uncertainty band.Not every simulation found an IMF that had a correlation > 0.5 with all the prescribed sinusoids.The 13 year oscillation had 848 matches, the 55 year only 550, and the 80 year 890.It appears that the extra mode or two that often pops up between 13 and 55 years in the EMD distorts the recovery of the 55 year signal (e.g., Fig. 3).
Although the phase of the mean 13 year IMF is consistent with the prescribed signal, the mean amplitude is too small (Fig. 4).The SD is also quite high relative to the amplitude (80 %), suggesting the actual recovered IMF could be nearly zero for any realization, or two times too large, depending on how the high-frequency variability affects it.
For the longer-period oscillations, there is a systematic error in the mean IMF.It is roughly the same in both the 55 and 80 year signal: the phase is only correct at the end of the record, and the amplitude is unrealistically increasing in time (Fig. 4), from almost no fluctuation at the beginning to larger variations than were prescribed at the end.The scatter is again relatively large compared to the largest amplitude (60-80 %).Finally, the acceleration estimated from the final IMF mode is systematically too small (Fig. 4).Introduction

Conclusions References
Tables Figures

Back Close
Full We acknowledge this test has its limitations.The final peaks of the 55 year 80 year sinusoids are very close to each other.However, a rather simplistic harmonic analysis using least squares over ranges of given periods found all three sinusoids precisely with small errors (< 5 mm).If EMD can purportedly uncover non-stationary oscillations in a data set accurately, then it should be able to compute perfectly stationary ones just as well.The level that EMD cannot do this is a reasonable estimate of its accuracy.
Figure 5 summarizes the results of Case 2, the simulation based on the ENSO and PDO indices.As with the experiment in Case 1, 1000 different simulations were run, differing only by the random noise.The IMF with the highest correlation greater than 0.5 with both the prescribed ENSO and PDO index was averaged.In addition, we found in nearly every case (99 %) the EMD computed an IMF with a periodic signal between the ENSO and the PDO signal.We isolated this signal by looking at the autocorrelation of the remaining IMFs uncorrelated with either PDO or ENSO.The IMF with an autocorrelation greater than 0.9 at a lag of 1 year was selected.
We should note that typically there were several IMFs that correlated significantly with the ENSO index.For the statistics shown in Fig. 5, only the one with the highest correlation was chosen.Although we found that by adding the 1-2 additional IMFs to the most significant ENSO mode resulted in a better correlation, we felt this was not a fair evaluation of the EMD process.ENSO is a physical process, and the relationship between the climate indices and the physics of the strength and timing of an ENSO event related to the index has been well established, (e.g., Philander, 1990Philander, , 2006)).
Although some authors have run EMD on ENSO indices and argued they have extracted distinct modes of interannual to multi-decadal variability (Wu and Huang, 2004;Franzke, 2009;Pecai et al., 2010), their conclusions are based solely on the fact such modes are extracted during the EMD process; they have offered no physical explanation for them.We note that other statistical based methods (such as principal component analysis) run on environmental data like sea surface temperature, precipitation, sea level, winds, etc. find modes highly correlated with ENSO and PDO indices (e.g., Mantua et al., 1997;Wolter and Timlin, 1998;Chambers et al., 1999;Bond et al., 2003).

NPGD Introduction Conclusions References
Tables Figures

Back Close
Full We know of none that find multiple modes that add up to correlate with an ENSO index.Thus, we argue it is more relevant to quantify if EMD can extract physically meaningful climate modes than whether it can extract modes with interannual and multi-decadal variability.
The ENSO-mode IMF on average matches the timing of the input ENSO variability (Fig. 5), although the amplitude is smaller; on average it underestimates the size of the El Niño and La Niña events by a factor of 2 to 3.Moreover, the SD is large, ranging from 50 to 250 % of the estimated peak values.This means that no single decomposition exactly matches the simulated ENSO variability.Some may catch a peak or two properly, but other El Niño or La Niña events are not captured at all.
The non-simulated low-frequency IMF has a period of between 25-30 years (Fig. 5), with an average amplitude ranging from 10 to 20 mm.This is the same magnitude of variability as the PDO-related variability, although IMFs extracted from a single simulation could have an amplitude nearly 3 to 4 times higher, based on the SD.Without knowing a priori what variations were in the data, this mode would be interpreted as a real, physical oscillation in sea level, when in fact it is a bogus artifact of the analysis.
As with the ENSO-mode, the mean PDO-mode IMF tracks the general periodicity of the PDO, although the amplitudes are on average too small.Again, the SD suggests any single simulation would give a considerable range of amplitudes.We note that as with the Case 1 results, there is a tendency for an increasing amplitude in time for the mean IMF, inconsistent with the true signal; the first two peaks in the given PDO signal are roughly the same magnitude.Finally, the average long-term rise computed from the last IMF is wrong (Fig. 5).The trend at 1900 is 36 % lower than prescribed, and the overall acceleration is 83 % higher.
Figure 6 shows the results from the EMD of Case 3, with the single extreme event.
Notice the large, non-stationary oscillation with a period of about 10 years in IMF6.The amplitude reaches 25 mm around 1956.Recall that this experiment only included seasonal variations, random noise, and this single large event.Because the EMD method implicitly assumes local highs are balanced perfectly by nearby lows, it cannot handle Introduction

Conclusions References
Tables Figures

Back Close
Full an extreme event like this.By enforcing an unrealistic balance of equal highs and lows, the method creates a low-frequency oscillation that does not exist.With a larger pulse, the magnitude of the error is even higher.It does not affect just this mode.It also shows up in IMF7 and IMF8, especially distorting the end of the record (Fig. 6).We have not tested by adding more extreme events, but we would assume it would cause even more spurious signals like these.

Conclusions
While at first glance empirical mode decomposition appears to be a useful tool for quantifying non-stationary multidecadal oscillations in sea level records, the results of our experiments suggest there are several issues.Probably the biggest one is the fact the EMD process applied to random noise consistent with high-frequency sea level variability and single extreme events will cause relatively large and systematic multidecadal oscillations that are not real.This will distort any underlying true signal.Our results suggest this is especially a problem for the longest period fluctuations; the IMFs are systematically biased away from the true signal, both in amplitude and phase.In some cases the amplitude increases in time, which could lead to incorrect interpretations regarding acceleration.
Moreover, there always appears to be one or more IMFs that are completely spurious fluctuations.These are needed to correct the errors in the other IMFs so they all sum to the original data.With no knowledge of the underlying physical modes, how is one to know which of the signals is spurious?In the articles that have applied EMD to sea level, all long-period IMFs have been assumed real and analyzed in regards to climatic or dynamical fluctuations in sea level.Based on the results of our experiments, we cannot believe that all the analyzed modes are true.
Finally, authors have asserted that the acceleration that comes out of an EMD process is more accurate, as they believe the IMFs better separate the high-and lowfrequency fluctuations than linear least squares.However, our experiments show the Introduction

Conclusions References
Tables Figures

Back Close
Full opposite.The quadratic fit to the last IMF is either no more accurate than one fit with least squares, or, in some cases, is significantly biased.In the experiment with ENSOand PDO-like oscillations, the acceleration estimated from the final IMF was nearly 100 % too large on average.In individual experiments, the error was even more.EMD is a quick and relatively easy tool to identify possible multidecadal fluctuations in a sea level record.However, it should not be used solely to quantify magnitude and phase, and one should be cautious in interpreting acceleration computed from the final IMF, especially in light of the significant errors found in the early and later parts of the low-frequency IMFs (Figs. 4,5,and 6).Instead, we believe other more traditional methods, such as harmonic analysis (e.g., Chambers et al., 2012), linear regression against climatic indices or physical parameters (e.g., Calafat and Chambers, 2013), running means of linear trends evaluated over discrete window-lengths (e.g., Holgate, 2007;Merrifield et al., 2009), or simply low-pass filtering on different timescales should also be utilized in order to find possible spurious signals in the IMFs arising from the way the EMD process filters random noise and extreme events.At the very least, authors should carefully remove extreme events from the sea level records before performing EMD to reduce biasing low-frequency IMFs.Unless other methods are utilized and shown to agree with the EMD results, we remain skeptical of many interpretations of EMD processed sea level data.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and t is the time in years, with dt = t − 1900.0.This baseline model is the same for all experiments.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .Figure 3 .Figure 5 .
Figure 1.Low frequency IMFs resulting from EMD of random noise with a SD of 60 mm.