Evaluation of Empirical Mode Decomposition for Quantifying Multi-decadal Variations and Acceleration in Sea Level Records

The ability of empirical mode decomposition (EMD) to extract multi-decadal variability from sea level records is tested using three simulations: one based on a series of purely sinusoidal modes, one based on scaled climate indices of El Niño and the Pacific decadal oscillation (PDO), and the final one including a single month with an extreme sea level event. All simulations include random noise of similar variance to high-frequency variability in the San Fran-cisco tide gauge record. The intrinsic mode functions (IMFs) computed using EMD were compared to the prescribed oscillations. In all cases, the longest-period modes are significantly distorted, with incorrect amplitudes and phases. This affects the estimated acceleration computed from the longest periodic IMF. In these simulations, the acceleration was underestimated in the case with purely sinusoidal modes, and overestimated by nearly 100 % in the case with prescribed climate modes. Additionally, in all cases, extra low-frequency modes uncorrelated with the prescribed variability are found. These experiments suggest that using EMD to identify multi-decadal variability and accelerations in sea level records should be used with caution.


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-periodic fluctuations.IMFs extracted from various tide gauge records have been 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, or nonlinear trend, than simply fitting a quadratic to the original data using ordinary least squares (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, and there is no reason to a priori expect a mode to have peaks and troughs of equal but offsetting magnitude.
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Thus, it is unlikely that a single IMF from the EMD analysis can represent a real, physical climate variation.Because of the assumptions underlying the method, it is more likely that multiple modes will be needed to quantify the physical climate mode.However, without some a priori knowledge of this mode, how can one know which IMFs to add together?In the worst case, the climate signal could be spread among a large number of modes.Already, several authors have performed an EMD on El Niño-Southern Oscillation (ENSO) indices and argued they have extracted distinct modes of interannual to multi-decadal variability (Wu and Huang, 2004;Franzke, 2009;Yang et al., 2010); their argument based solely on the fact that such modes are extracted during the EMD process, but with no physical explanation for them.The same has been done to sea level measurements made by tide gauges, and individual IMFs are interpreted as distinct climatic modes (Ezer and Corlett, 2012;Ezer et al., 2013).
Moreover, Wu and Huang (2004) have previously shown that EMD behaves as a low-pass filter on random noise.If one runs white noise with a normal Gaussian distribution through the process, low-frequency signals will appear in the resulting IMFs.They found there is roughly a doubling of the average period with each subsequent mode.This is a significant issue.It means that any quasi-random, high-frequency signal will propagate into low-frequency signals in the recovered IMFs.Wu and Huang (2009) proposed a method to quantify the uncertainty caused by this behavior by computing an ensemble mean of IMFs, starting from the same time series but with different amounts of added random noise.
However, this method ignores the fact that most geophysical time series have an underlying real signal that has high variance and little serial correlation, i.e., a high-frequency, near-random signal.This signal will also be filtered by the EMD process and will likely appear as a quasi-stationary oscillation in higher-order IMFs that is not real.Although adding multiple realizations of random noise to a time series will account for uncertainty in the IMFs from random error in the measurement, it will not account for the shifting of high-frequency signal to low-frequency signal in the recovered IMFs.One of the assumptions in EMD processing is high-frequency variability is captured in the lowest IMFs, but as far as this author is aware, this assumption has not been evaluated or verified.
Additionally, the assertion that the recovered nonlinear trend from EMD is more accurate than one computed using a parametric model and ordinary least squares has not been evaluated for data that simulates a tide gauge record.Considering the importance of quantifying acceleration in long sea level records to understand ongoing climate change, this is a vital test.Franzke (2011) conducted an experiment of detecting nonlinear trends (i.e., an acceleration) in a small suite of 100 simulated temperature time series, using different statistical estimators, including ordinary least squares and EMD.The results showed no statistically significant improvement using EMD.In fact, in most tests, the nonlinear trend esti-mated using ordinary least squares was closer to the input signal.Whether the same result holds for sea level records is still an open question.
Several questions arise from this discussion.How well can the EMD method recover the acceleration in a long tide gauge record?Is it more accurate than using a linear model and ordinary least squares?Do the individual IMFs reflect distinct climate modes, or do they reflect in part the aliasing of high-frequency variability to the low-frequency because of the EMD low-pass filtering?To answer these questions, we will apply the EMD process to three simulated data sets where known low-frequency modes are prescribed.This is not a novel idea, and should be used to evaluate any new algorithm.However, it has not been used in the application of EMD to sea level records to this author's knowledge.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 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 standard deviations.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 with significant multi-decadal variability 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.These spurious low-frequency IMFs also have a tendency to bias the recovered acceleration either low or high.

Data and methods
The basic idea of EMD is to fit cubic splines to the local maxima and minima of a time series 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 standard deviation 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).
EMD is applied 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.A tolerance value of 0.05 was utilized.
There is 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 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 (e.g., Ezer and Corlett, 2012;Ezer et al., 2013), which shows a continuously increasing sea level, is to fit a quadratic to the final IMF from the EMD process using least squares, and plot the resulting fit.This is conceptually no different than fitting a quadratic to the original time series, other than the fact that it is done to the final mode, which has significantly lower variance than the original data.This should improve the estimate -provided there are no systematic errors or biases in the final IMF that would bias the result.
To demonstrate how EMD shifts some of the random variability to higher IMFs (Huang and Wu, 2008), we ran EMD on a monthly resolution time series that is 150-years long with randomly correlated values that have a standard deviation of 60 mm, using a white noise model with a normal Gaussian distribution.A value of 60 mm was used because this is the standard deviation of the residual monthly sea level at San Francisco after fitting and removing a quadratic function plus annual and semi-annual sinusoids; therefore, it is representative of high-frequency sea level at a typical site, although some sites can have significantly higher variability.Another case was run using a colored noise model that reproduces the autocovariance of the San Francisco tide gauge residuals, based on an AR(5) model, where the coefficients are computed from the autocovariance following the Yule-Walker method.The results were nearly identical to the ones shown with the randomly correlated residuals; subsequently, we choose to use random values as they are faster to compute for the several thousand simulations planned.The EMD finds IMFs that have quasi-periodic fluctuations of nearly 60 years and amplitudes as large as 10 mm (Fig. 1); fluctuations at quasi-30-year periods are the same magnitude.
We use the monthly sea level record from the San Francisco tide gauge for our reference.It was downloaded from the Permanent Service for Mean Sea Level (PSMSL) (Woodworth and Player, 2003;PSMSL, 2012).Annual and semiannual sinusoids are fit to the data along with a trend and an acceleration term using ordinary least squares to obtain the baseline sea level variability for the model (y base ), where and t is the time in years, with dt = t − 1900.0.This baseline model is the same for all experiments.For case 1, three long-period sinusoids (13, 55, and 80 years) are added to the baseline model along with white (random) noise which has a normal distribution (ε(t)): (2) The random noise has a variance to match the variance of the residuals of the real tide gauge data minus the model.In all, 1000 different random noise models were applied to create 1000 different versions of case 1 to quantify how the recovered IMFs change depending on the different high-frequency variability.The periods and amplitudes of the long-period sinusoids were chosen arbitrarily to approximate the level of multi-decadal fluctuations in the San Francisco sea level record (Fig. 2a).The hypothesis being tested is that the highfrequency variations are isolated into the lowest IMFs with little or no distortion of the higher IMFs, and that the higher IMFs will represent the prescribed multi-decadal fluctuations.
Case 2 starts from the same baseline model, but instead of prescribing sinusoidal oscillations, non-stationary climate indices for ENSO variations and the PDO are used.The Southern Oscillation Index (SOI) utilized is 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 is 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 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 the PDO record to fill in the missing data before January 1900.Recall the experiment does not require true ENSO or PDO variability, only a simulation of the type of variability and how well EMD can recover it.
Second, because the two indices are slightly correlated (−0.21, p < 0.001) due to similar interannual (< ten year) variations, the PDO index is low-pass filtered with a 5-year Gaussian, and the SOI is band-pass filter 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, p < 0.01).
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 standard deviation.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 standard deviation of the filtered residuals is calculated.The scaling factor applied to PDO LP (t) is the standard deviation of the low-pass filtered sea level residuals (20.8 mm); the scaling factor applied to SOI BP (t) is the standard deviation of the band-pass filtered sea level residuals scaled by −1 to account for the fact that the 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 standard deviation of the residuals and the model: where 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 standard deviation of 60 mm (representative of the high-frequency variability in the San Francisco sea level), then adds an extra 350 mm for January 1956 to represent the signal of a large anomalous high-water event, such as the effect of a large storm surge event on the monthly average, a large flooding event from sustained rainfall, or climatic variations in winds that can cause sustained high water levels.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 standard deviation).Moreover, most 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.

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 of a much higher frequency and so are not considered, although we 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 IMF 6 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 (IMF 8) 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 (IMF 9 in this case), significantly underestimates the prescribed acceleration.
Although this is a single example, it will reflect the type of the distortion in low-frequency IMFs caused by apply-ing the EMD algorithm to high-frequency variability inherent in a sea level record.Note that the ensemble EMD approach proposed by Wu and Huang (2009) creates the ensemble members from the original time series, differing only by random noise.This process will still filter the inherent, highfrequency, quasi-random signal with EMD, which will likely bias the ensemble mean.Moreover, Wu and Huang (2009) assume that averaging will minimize any residual effect of EMD from the additional random noise on the ensemble mean, although this is not demonstrated.
We test whether this assumption is valid by averaging the 1000 different IMF clusters computed from the simulations.One cannot simply separate the corresponding low- frequency modes based on the IMF number, however, as the total number of IMFs changed from 9 to 11 in the 1000 different simulations.The 13-year signal was found in IMFs ranging from number 4 to 7, while the 55-year mode ranged from IMF 7 to 9. The last mode found ranged from IMF 9 to 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.Figure 4 shows the histogram of computed correlations.Note that the correlations were not the same in every simulation.The 13-year oscillation had a mean correlation of 0.66 (standard deviation = 0.09), the 55-year had a mean of 0.52 (standard deviation = 0.11), while the mean correlation of the 80-year signal was 0.74 (standard deviation = 0.09).
So that the lower correlation in the 55-year test did not bias our results, a minimum bound was set to 0.5.If two or more modes have correlations > 0.5 with one of the input signals, the one with the highest correlation will be chosen.Figure 5 summarizes the results, showing the mean IMF with the standard deviation 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 oscillation only 550, and the 80year oscillation 990.It appears that the extra mode or two that often pops up between 13 years 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. 5).The standard deviation is also quite high relative to the amplitude (80 %), suggesting the actual recovered IMF could be nearly zero for any realization, or 2 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 55and 80-year signal: the phase is only correct at the end of the record, and the amplitude is unrealistically increasing in time (Fig. 5), 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. 5).
We acknowledge this test has its limitations.The final peaks of the 55-year and 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).The fact that EMD creates non-stationary modes where none are present is troubling, and suggests one must be very careful in interpreting the results for a single IMF.
For example, consider the interpretation of the results from this simplistic simulation in terms of longer-term climate change if only the EMD results (Fig. 5) were analyzed.Based only on the returned IMFs, one could easily argue that there was no significant low-frequency variation in the sea level before 1950, then a rather dramatic rise in the 1970s, followed by a return to a normal condition.In fact, there were equally large sea level shifts in the early part of the simulated record that were lost due to the way the EMD method partitions the real signal.
Figure 6 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, in nearly every case (99 %) the EMD computed one to two IMFs with a periodic signal that did not correlate highly with either PDO or ENSO, but had a low-frequency.Because this was not always contained in a single IMF between the two prescribed periodic fluctua- tions, we had to adapt a method to search for it.This signal was isolated by looking at the autocorrelation of the IMFs after removing those correlated with PDO or ENSO, as well as the last mode.To find the mode with the longest-period fluctuation, we examined the autocorrelation at a 1-year lag.Only IMFs with an autocorrelation greater than 0.9 at a lag of 1 year were examined, and if two existed, the one with the higher autocorrelation was selected.
We should note that typically there were several IMFs that correlated significantly with the ENSO index.For the statistics shown in Fig. 6, 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;Yang 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, e.g., on environmental data like sea surface temperature, precipitation, sea level, winds, 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).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. 6), 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 standard deviation 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 and 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 standard deviation.Without knowing a priori what variations were in the data, this mode would be interpreted as a real, physical oscillation in the 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 standard deviation 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, which could be misinterpreted as a sign of climate change; 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. 6).The trend at 1900 is 36 % lower than prescribed, and the overall acceleration is 83 % higher.
Figure 7 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 IMF 6.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 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.Although the random-only case (Fig. 1) also produces low-frequency erroneous oscillations, the amplitudes are significantly less for the longer-period IMFs.With a larger pulse, the magnitude of the error is even higher.It does not affect just this mode.It also shows up in IMF 7 and IMF 8, especially distorting the end of the record (Fig. 7).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 multi-decadal 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 up 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 low-frequency fluctuations than applying a parametric model and linear least squares.Their argument assumes that the high-frequency variations and shorter-period non-stationary signals in the original time series are biasing a quadratic fit to the original data.By distributing these signals in the EMD process to specific IMFs, they believe the final IMF contains the true acceleration plus residual low-frequency variability.Even Fanzke (2011), who demonstrated that EMD was no better than using an ordinary least squares estimator and a parametric model, argued that EMD was still better if the trend was nonlinear, especially exponential.Our experiments, however, show the opposite.The quadratic fit to the last IMF is either no more accurate than one fit with least squares to the full, unfiltered data set, or, in some cases, is significantly biased.In the experiment with ENSO-and 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.This is most likely due to the aliasing behavior of EMD where some of the high-frequency variance is aliased into the low-frequency modes, as we have demonstrated.
EMD is a quick and relatively easy tool to identify possible multi-decadal fluctuations in a sea level record.However, we have demonstrated that real climatic non-stationary signals are generally spread among multiple modes.Analyzing a single IMF for climate variability will likely lead to significantly biased interpretations.Thus, we feel that EMD analy-sis should not be used solely to quantify magnitude and phase of non-stationary climate variations, nor should analysis of climatic signals be based on a single IMF.One should also 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. 5,6,and 7).Where EMD has shown to be useful has been in lowpass filtering data to reduce the impact of high-frequency variability and noise (e.g., Alberti et al., 2014).In that case, the sum of the higher IMFs are used as the low-pass filter.
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 along with EMD to study low-frequency climatic variability.This is 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.

Figure 1 .
Figure 1.Low-frequency IMFs resulting from EMD of random noise with a standard deviation of 60 mm.

Figure 3 .
Figure 3. True oscillations and long-term trend + acceleration (red) for simulation shown in Fig. 2a, along with the closest correlating IMF (blue).

Figure 5 .
Figure 5. Mean (solid blue line) and standard deviation (light blue envelope) of IMFs calculated from the 1000 different case 1 simulations, along with the true signal (red).

Figure 6 .
Figure 6.Mean (solid blue line) and standard deviation (light blue envelope) of IMFs calculated from the 1000 different case 2 simulations, along with the true signal (red).

Figure 7 .
Figure 7. Low-frequency IMFs resulting from EMD of a simulated signal with a trend + acceleration + seasonal variations + random noise + a single large anomalous event in January 1956.