Simplifying Laplace's Tidal Equations for QBO

I have been playing the models of ENSO and QBO off of each other. So for whatever insights I glean from one model, I have tried to apply it to the other.  My last ENSO analysis involved simplifying Laplace's Tidal Equations and finding similarity to the Mathieu-modulation used to describe sloshing of water volumes.  However, I had left one equation in a rough state, and stated that I would get back to it.  That's what I will do in this post, and I can flatly assert that the simplification is an eye-opener for understanding the behavior of QBO! (Which shouldn't be surprising since Richard Lindzen tried vainly to apply Laplace's equations to QBO, with the outcome of creating a monster of complexity, see  [1])

Continue reading

Ridiculous Paper on Ocean Oscillation

These scientists thought that they could increase interest in their research by converting a very low frequency ocean oscillation into a human audible tone.  They then called it a "Rossby Whistle".  As they write "this resonant mode is dynamically equivalent to the operation of a whistle."

Hughes, C. W., J. Williams, A. Hibbert, C. Boening, and J. Oram (2016), A Rossby Whistle: A resonant basin mode observed in the Caribbean Sea, Geophys. Res. Lett., 43,  doi:10.1002/2016GL069573.

This is really silly because converting it to a human audible tone is meaningless, as the original frequency has a period of months, which is absurdly below the human audible detection limit of 10 to 20 Hz at the low-frequency end of the audio spectrum. And so it is only intended to illicit a response in people that might find it spooky sounding.

Here is my take on comparing a tone with harmonics of the fundamental frequency to that of a converted Quasi-Biennial Oscillation (QBO) waveform.

This is the spectrum of the harmonic tone; note the octave doubling in the peaks starting with the harmonic at 270 Hz.

harmonic_spectrumThis is what the QBO looks like. These are not harmonics but are aliased frequencies of the original tidal frequency.

qbo_spectrumThis is a scope trace of the harmonic tone. The envelope is the fundamental and the harmonics provide the finer structure, forming a regular pattern.

harmonicAnd this is a scope trace of the equivalent QBO tone, showing the more erratic nature of combining the aliased frequencies.

qboAnd this is what they sound like in comparison. The first tone is the harmonic tone, and then separated by a pause is the QBO tone.  The first sounds pleasing and the second converted QBO waveform sounds out-of-tune because the frequencies aren't harmonics.

As John Gielgud once said, "I'll alert the media". Sigh. In our current media circus, rossbyWhistlereporters will lap up the inconsequential aspect of the strange sound and ignore the fact that Richard Lindzen couldn't figure out that the QBO comes about from seasonal aliasing of the lunar nodal cycle.  Even though he knew that a lunar tidal mechanism could form a valid premise for atmospheric oscillation.

Example of the second-derivative of the QBO, which exaggerates the aliased frequencies.

Still can't figure out why Lindzen never saw this.

Recent research say QBO frequency is emergent property

Gavin Schmidt of NASA tweeted a cite to a recent paper on QBO by a group of 12 scientists

Geller, M. A., Zhou, T., Shindell, D., Ruedy, R., Aleinov, I., Nazarenko, L., Tausnev, N.L., Kelley, M., Sun, S., Cheng, Y., Field, R.D. and Faluvegi, G. (2016), Modeling the QBO – Improvements Resulting from Higher Model Vertical Resolution. J. Adv. Model. Earth Syst.. doi:10.1002/2016MS000699

They claim that the fundamental frequency of QBO relates inversely proportional to the pressure:

"The mean zonal winds from the ERA-Interim reanalysis are shown between 4.5 and 4.5 N between heights of 100 and 1 hPa for the 20 years 1991-2010. Note that ERA-Interim shows that approximately 8 1/2 QBO cycles occurred during this period, implying a QBO period averaging about 28 months. The model results in figure 1b show that no coherent QBO resembling observations exists for the gravity wave momentum flux forcing of 1.5 mPa, which is consistent with the steady jets at low wave forcing demonstrated by Yoden and Holton (1988). It also shows that the QBO-like oscillation for a gravity wave momentum flux forcing of 2.0 mPa has a period of about 8 years, while a forcing of 2.5 mPa gives a period of about 37 months, and a forcing of 3.0 mPa gives a period of about 25 months, and a forcing of 3.5 mPa gives a period of about 21 months. In fact, we find that the best fit to observed QBO periods is for a gravity wave momentum flux forcing of 2.9 mPa, as will be shown in the next section. "

I disagree with this interpretation. I know that this might sound like I am using the argument of personal astonishment, but I can't see how a physical model could operate this way. If anything, pressure may modify the amplitude of the oscillations, i.e. the peak speed, but asserting that the pressure sets the oscillation period is likely the result of a model that no one understands. If it was indeed caused by the pressure, they should be able to come up with a simplified formula, like deriving the resonant frequency of a Rijke tube, instead of suggesting that it is an emergent property of a complicated model. Certainly, pitch can change as a result of say overblowing a woodwind instrument, but even there, the shifts are more than likely harmonics or overtones of the original fundamental frequency, and not a continuous change in frequency.

RHEED oscillations more here

And I realize that the QBO is obviously not an acoustic oscillation yet the physical model for any oscillation does need some sort of mathematically intuitive explanation (see figure to the right), which is greatly lacking in this paper.  To say it emerges from the model is an extremely weak argument.

The 28 month period is ubiquitous in a number of geophysical oscillations. It shows up in (1) the seasonally aliased Draconic lunar cycle, (2) the strongest factor in the model of ENSO, (3) exactly half that period in the Chandler Wobble, (4) its defined as the mean annual flood recurrence interval by the USGS (first reported in 1960), and (5) the fundamental QBO period.

Perhaps no one wants to admit that seeing this number show up in all these measures is most likely a forced response from the nodal lunar cycle factor, as opposed to a natural resonance. The reality is that a forced response is not as "sexy" as finding a natural frequency eigenvalue.

But of course there is no way to verify either way, since a controlled experiment is not available to test a model against.  There is only one earth and it can't be duplicated in the lab. All we have to go on is the empirical data and any model that we can dream up.

A bunch of meteorologists are freaking out that the current QBO data is not meeting up to expectations

I don't understand the rarity of this since the higher altitude QBO aligns more with semiannual seasonal cycles, while the lower altitude QBO aligns more with the nodal lunar cycles. These will constructively add according to a beat frequency and so the alignment will occur occasionally. This summer solstice there is a full moon and that only happens once every 19 years.  I can see this with my QBO model as well, so what is the big deal?

Pukite's Model of ENSO

Putting the finishing touches on the ENSO model, I decided to derive the Mathieu modulation of the wave equation. It's actually a rather obvious simplification of Laplace's Tidal Equations or the Shallow-water Wave Equation applied along the equator. That last assumption is important as the Coriolis forces are known to cancel at the equator (added: perhaps more precisely Isaac Held explains that it vanishes at the equator -- it is both zero and a reversal) and so we can use a small angle approximation to reduce the set of equations to a Mathieu or Hill formulation.

To make the job even easier, the three equations that comprise the the Laplace formulation can be formulated initially in terms of the Laplace tidal operator. Per the recent reference [1], the authors separate the temporal component from the standing-wave spatial component(s) and express the operator concisely via implicit time and wavenumber (see Figure 1).

Fig 1:

Fig 1:  From [1], formulation of the Laplace Tidal Equation

Although the paper proceeds to detail the computational approach to solving the tidal equation, we will only piggy-back off the formulation and arrive at a gross simplification.

The first step is to use the small angle approximation to \phi

\mu = sin \phi \cong\phi

If  \sigma \gg \mu while \sigma not a function of \mu , then the operator can be reduced to:

F(\Theta) = \frac{1}{\sigma^2} \left( \frac{d^2\Theta}{d\mu^2} - \left[\frac{s}{\sigma} + s^2\right] \Theta \right)

The next crucial step is to extract the implicit temporal change from \mu via the chain rule.  Even though \mu is small, changes in \mu have to be retained to first-order or else the solution would be static along the equator.

\frac{d^2\Theta}{d\mu^2} = \frac{d^2\Theta}{dt^2} \frac{d^2t}{d\mu^2}
added 6-17-2016: In spirit, the above is close, but not a correct expansion.
The actual chain-rule expansion is a bit more complicated and features first-derivative cross-terms.
Will get to that later, not a big deal.

The second term has the dimensions of a reciprocal latitudinal acceleration. If we assume that  \mu = \mu_0 sin(\omega t) , then
F(\Theta) = - \frac{1}{\sigma^2} \left( \frac{d^2\Theta}{d t^2} \frac{1}{\mu_0 \omega^2 sin(\omega t)} + \left[\frac{s}{\sigma} + s^2\right] \Theta \right)

With that simplification, we reduce the formulation to a temporal second-order differential equation.

\frac{1}{\sigma^2} \left( \frac{d^2\Theta}{d t^2} \frac{1}{\mu_0 \omega^2 sin(\omega t)} + \left[\frac{s}{\sigma} + s^2\right] \Theta \right) + \gamma \Theta = f(t)

where f(t) is a forcing added to the right-hand side (RHS).

Moving the sin(\omega t) from the denominator, we recover the Mathieu-like wave equation formulation:

\frac{d^2\Theta}{dt^2} + \left[\frac{s}{\sigma} + s^2 + \sigma^2 \gamma \right] \Theta \mu_0 \omega^2 sin(\omega t) = \sigma^2 \mu_0 \omega^2 sin(\omega t) f(t)

For now, the primary objective is to verify the sinusoidal behavior in terms of correlating periodicities. Previously, we determined that the f(t) term was composed of seasonally aliased lunar periods, such as the Draconic 27.2 day cycle. (Recently I also found this term buried in the LOD data.)

And the Mathieu modulation likely includes an annual contribution, and, importantly, a strict biennial cycle \omega=\pi , which is both observed and implied based on other climatic measures.

The seasonally aliased 27.2 day cycle generates harmonics in ENSO at the same locations as for QBO
1/T = 1/P + n, where P is the tidal period in years, and n is ... -3, -2, -1, 0

so for Draconic, 1/T = 365.242/27.2 = 13.43. Stripping off the integer, the value of 0.43 is apparent (and 0.43-1=-0.57) as the strongest low-frequency factor in the power spectrum of Figure 2. Similarly for the fortnightly Draconic (13.6 day), 1/T = 365.242/13.6 = 27.85, and stripping off the integer, one gets 0.85 (and 0.85-1=-0.15). This also corresponds to the Chandler wobble. And with the multiplicative biennial modulation factor adding a +/- 0.5 to each aliased value, the other candidates are 0.85-0.5 = 0.35, 0.43+0.5 = 0.93 (and -0.57+0.5=-0.07, -0.15-0.5=-0.65). All these values are observed as among the strongest low-frequency components in the power spectrum below.

Fig. 2:  The strongest components correspond to the seasonally aliased Draconic period. The arrows separate the biennial symmetric terms.

The model fit to the ENSO spectrum is shown in Figure 3. This includes more than the Draconic terms, including potential Tropical fortnight (27.32/2 day) and Anomalistic (27.55 day) terms. Apart from the 1/f filtering envelope, the symmetry about the f=0.5/yr is evident. Note that the power spectrum calculated is essentially the LHS response (data) against the RHS forcing (model). The second derivatives and modulation are included in the LHS, so it is not completely "raw" data.

Fig. 3:  ENSO power spectrum comparison to model

And in temporal space, the multiple regression fit is shown in Figure 4. Interesting that the correlation coefficient is higher outside the training interval than within the training interval!

Fig. 4: Fit using the training interval from 1881 to 1950.

This is a very elegant theory in terms of deriving from first-principles; the fact that it is a simplification direct from Pierre Laplace's original tidal formulation means that we have bypassed sophisticated GCM approaches to potentially solving ENSO .  The only caveats to mention are the phase inversion from 1981 to 1996, and the imprecision of the  Draconic term, where 27.199 days fits better than 27.212 days.  The former aligns with a 7 year interval as there are precisely 94 cycles of 27.199 day periods in 7 years.  As ENSO is known to be seasonally aligned, this makes qualitative sense. However, the phase error from that cycle to the actual Draconic cycle will build up over time, and it will eventually go completely out of phase after 70+ years.  This actually may be related to the 1981 phase inversion, as that would be a means of re-aligning the phase after a number of years (Another phase inversion may occur prior to 1880, see comment).

The white paper is located on arXiv.

I have to get to work so will leave this post subject to fixing any Latex markup typos later. (Note: fixed some sign-carryover errors due to the 2nd-deriv, same day as posting)


[1] Wang, Houjun, John P. Boyd, and Rashid A. Akmaev. "On computation of Hough functions." Geoscientific Model Development 9.4 (2016): 1477. PDF

QBO Aliasing Graph

Before getting to the needless complexity discussion hinted at in the last post, here is a simple graphical analysis of finding aliased frequencies. Based on this analysis one can determine aliased signals by looking at the frequency power spectrum. The rules for aliasing in frequency space involve peeling back from the assumed frequency until the aliased harmonics appear in the spectrum.

1/T = 1/P + n, where P is the tidal period in years, and n is ... -3, -2, -1, 0

So applying the Draconic period aliasing to QBO, the harmonics are clearly seen in the following Figure 1.  The strongest factor is the approximately 2.33 year period shown on the left, but since this is a second derivative of the QBO time series, higher frequency aliased harmonics are also observed to the right.  Note that the factors are approximately 0.43, 1.43, and 2.43, which are separated by Δn=1. And note that 1/0.43 ~ 2.33 years, which is the approximate observed QBO frequency..

Fig. 1:  Spectrum of the second derivative of QBO, showing seasonally aliased factors

Other potential periods are shown in the figure, including the unaliased annual and semi-annual signals separated by the violet colored arrow. The green arrow separates the harmonics of the 13.633 day lunar tide periodicity. Clustered around the strong Draconic period are also satellite frequencies. At least one of these satellites along with the strong aliased Draconic also appear in the LOD spectrum. That provides the physical basis for QBO (and also ENSO) forcing.

Putting these factors together via multiple regression results in marvelous fits of the various QBO time series, as seen in Figure 2 below.

Fig. 2:  Multiple regression fit of the 30 hPa QBO time series.

And to highlight the second derivative composition and how the higher frequency aliased harmonics contribute to the fine detail, see Figure 3.

Fig. 3 : Second derivative of QBO highlights the higher frequency aliased harmonic detail.

Only in a few places do the aliased harmonics not appear well aligned. Those years, such as the interval 1983-1984 may be telling us that other climate events may be playing a role in temporarily impacting the QBO dynamics.

Seasonal Aliasing of Long-Period Tides Found in Length-Of-Day Data

From a careful analysis of Length-of-Day (LOD) measurements, we can glean lots of information about phenomenon such as tides, earthquakes, core/mantle motion, and ENSO [1].  In Figure 1 below, one can see how this is manifested at different time scales.

Fig. 1: Differing temporal scales of LOD analysis.  Short period tides have scale of diurnal and semidiurnal.

I revisited the LOD data to see if I could detect the angular momentum changes that appear to force the QBO and ENSO oscillations, similar to what Wang et al do in [2].  In each of these natural behaviors, there is a clear indication that the 2.369 year seasonally-aliased long period draconic tide is a primary forcing mechanism.

There are two sets of data to look at. The original set from Gross at JPL called LUNAR97 [3] is not online in a machine-readable format, but it is included in the paper. So it was a simple matter to OCR a GIF of the paper's table. Another more recent set is available from IERS here.

Fig. 2: Compariscon between the two LOD data sets.  The Gross set has less filtering

The issue with the more recent set is that the site claims that they do a "5-point quadratic convolute", at least up to 1955.  I don't want any filtering, so chose the Gross LUNAR97 set. You can see the difference between the two data sets in Figure 2 above.

To draw out the fine features even more, I performed a straightforward second-derivative on the LUNAR97 time-series (same approach as used on the QBO see result). This was then fed into a symbolic regression machine learning tool to determine the principle components. The results are shown in Figure 3 below.

Fig. 3:  Result of machine learning on 2nd-derivative of LOD data.

After about an hour of machine learning, the solution included a strong frequency at 22.481 rads/year.  In Figure 3, that solution was one of the simplest along the Pareto front, shown by the red dot. This frequency is aliased to the slower frequency of 2.3694 rads/year, which is correspondingly seasonally aliased to an approximate lunar month of 27.2121 days. Incredibly, but perhaps not surprisingly, this compares to the known draconic or nodal lunar month length of 27.21222 days.   The other somewhat weaker term is very close to aliasing the anomalistic lunar month of 27.55 days, in the fortnight mode of 13.777 days.

That was the first run I attempted and subsequent runs are picking out some other potential long-period lunar periods (and also the 6.45 year Chandler wobble), but this is really enough to demonstrate a salient result -- that the same draconic period that the QBO model uses and that the ENSO model uses as the primary angular momentum forcing term also appears in the most direct and practical measurement of angular momentum known in Earth geophysics.   And I do not think this is a instance of artificial sampling aliasing, since each of the LOD measurements was averaged over at least a month.  It is likely a real seasonal aliasing behavior and serves to further substantiate the theorized mechanisms forcing the QBO and ENSO behavior.   Read the short paper that I placed into arXiv to see how the other parts fall into place.

And a continuing lesson to be learned is to not filter data unless you can be sure that what you are filtering out is nuisance noise. In this case, it appears the filtering performed by IERS was too aggressive and removed a very important signal in the time series! Contrary to what the statistician Tamino has been preaching, these time-series do not have a lot of statistical content. I donated money to his lessons on time-series analysis, but largely disagree on the applicability of his approaches to earth data such as we see in QBO, ENSO, Chandler wobble, and now LOD. If you followed his instructions, you might not find any of the periodic cycles buried in the data, and you might dismiss it all as red noise. Yet, there is likely nothing statistical or complex about these behaviors, as they are largely driven by known deterministic forces. Stay tuned, as the next post will be titled "Needless Complexity".


[1] B. Fong Chao, “Excitation of Earth’s polar motion by atmospheric angular momentum variations, 1980–1990,” Geophysical research letters, vol. 20, no. 3, pp. 253–256, 1993.

[2] Shen, Wen-Jun Wang Wen-Bin, and Han-Wei Zhang. "Verifications for Multiple Solutions of Triaxial Earth Rotation."  PDF

Note: the authors in [2] also find a 14-year signal in the LUNAR97 data (see figure below). This is likely a result of a wavelet filter that isolated periods between 8 and 16 years as they described the windowing procedure.  This 14 year period is critical in the ENSO model.

From Wang et al, Figure 6 "Spectrum of wavelet for data series LUNAR97 with evident peak at period 14 yr ".   The reciprocal of 0.0715 cycles/year is 14 years.

[3] Gross, Richard S. "A combined length-of-day series spanning 1832–1997: LUNAR97." Physics of the Earth and Planetary Interiors 123.1 (2001): 65-76.  PDF



Biennial Stratosphere Mode

In addition to the biennial mode found in ENSO and in GPS readings, a stratospheric biennial mode also exists. This is different than the quasi-biennial oscillation (QBO) modeled previously, as it shows a more strict two year cycle.

From Dunkerton [1]:

Second, concerning their time dependence, subbiennial variations should not be viewed in isolation from other modes. The seasonal dependence of QBO anomalies cannot be described by a single harmonic, whether quasi-biennial or subbiennial; rather, their superposition provides for the seasonality.


Only three parameters were varied to obtain the least squares fit: the QBO period, QBO phase, and subbiennial phase. This calculation was performed on the three tracers independently, but the same QBO period was obtained in each case, namely, 26 months. The corresponding subbiennial period is 22.3 months

What Dunkerton is talking about essentially is akin to a 2.166 year quasi-biennial signal mixed with a 1.858 year sub-biennial signal. These combine in frequency space as Dunkerton states by their superposition, and also couple together via a ~26-year modulation on the biennial signal -- as per the trig identity:

cos(2\pi t/2.166) + cos(2\pi t/1.858) = 2 cos(\pi t) cos(2\pi t/26)


This is substantiated by Remsberg in his methane study [2] where in addition he sees a smaller amplitude factor which is closer to the exact biennial signal.

Three interannual terms are almost always present and are 19 included in the model: an 838-dy (27.5-mo) or QBO1 term; a small amplitude, 718-dy (23.5-mo) biennial or QBO2 term; and a 690-dy (22.6-mo) sub-biennial (SB) term, whose period arises from the difference interaction between QBO and annual terms. The relative amplitudes of these three interannual terms vary with latitude and altitude.

More recently (in 2014) Remsberg reports [5]:

Instead, Fourier analysis of the time series of the residuals after removing the seasonal terms almost always indicates that there are two significant, interannual terms having periods of order of 28 (QBO-like) and 21 months (subbiennial term denoted as IA).

This is akin to a 2.333 year quasi-biennial signal mixed with a 1.75 year sub-biennial signal, leading to a 14-year signal modulating the biennial signal:

cos(2\pi t/2.333) + cos(2\pi t/1.75) = 2 cos(\pi t) cos(2\pi t/14)


Which is strikingly similar to the behavior observed in the ENSO model's 14-year modulation of the biennial cycle. The idea was that this 14 year modulation is close to the additional triaxial wobble suggested by Wang.  Slightly weaker is an 18.6 year modulation which is close to the lunar nodal period.

Fig. 1:  Biennial modulation in the ENSO power spectrum

The similarity between the stratospheric measures and the ENSO observations is based on the common behavior of a symmetrical balance about the biennial cycle. This is starkly observed in the power spectrum of ENSO shown above.

In a closely related finding, the single-minded AGW-denier Murry Salby also discusses this modulation in [6], but ascribes the modulation to the sunspot cycle, as he finds a value closer to an 11-year modulation. That assumption has to be taken with some circumspection, as it is well known that Salby favors a solar-based rationale to global temperature variation, which has gotten himself into some bizarre predicaments.


1. Dunkerton, Timothy J. "Quasi-biennial and subbiennial variations of stratospheric trace constituents derived from HALOE observations." Journal of the atmospheric sciences 58.1 (2001): 7-25.
2. Remsberg, Ellis E. "Methane as a diagnostic tracer of changes in the Brewer–Dobson circulation of the stratosphere." Atmospheric Chemistry and Physics 15.7 (2015): 3739-3754.
3. Remsberg, E. E., and L. E. Deaver. "Interannual, solar cycle, and trend terms in middle atmospheric temperature time series from HALOE." Journal of Geophysical Research: Atmospheres 110.D6 (2005).
4. Remsberg, Ellis E. "On the observed changes in upper stratospheric and mesospheric temperatures from UARS HALOE." NASA Langley Report (2006).
5. Remsberg, E. E. "Decadal-scale responses in middle and upper stratospheric ozone from SAGE II version 7 data." Atmospheric Chemistry and Physics 14.2 (2014): 1039-1053.
6. Salby, Murry, Patrick Callaghan, and Dennis Shea. "Interdependence of the tropical and extratropical QBO: Relationship to the solar cycle versus a biennial oscillation in the stratosphere." Journal of Geophysical Research: Atmospheres 102.D25 (1997): 29789-29798.

\frac{1}{2.333} + \frac{1}{1.75} = 1


Biennial Signals from GPS

This paper is from an open-access journal, yet it's a mind blower.

Pan, Yuanjin, et al. "The Quasi-Biennial Vertical Oscillations at Global GPS Stations: Identification by Ensemble Empirical Mode Decomposition." Sensors 15.10 (2015): 26096-26114.

GPS signals have enough information content that they can pick up sensitive measures related to earth deformation, so the paper describes the results.

That strict biennial signal is readily apparent along even years and it likely has everything to do with the biennial-modulated ENSO behavior. Have to remember that a thermocline separating two-slightly different density liquids acts as a highly sensitive sensor -- any slight forcing will get the sloshing in motion, and the GPS is likely isolating the biennial component.

Also some indication that the measure is picking up the 28 month QBO signal. See Table 1 in the paper, which shows an average period of 2.3 years for the quasi-biennial oscillation.

Biennial Connection to Seasonal Aliasing

I found an interesting mathematical simplification relating seasonal aliasing of short-period cycles with a biennial signal.

If we start with a signal of an arbitrary frequency \omega_L = 2\pi/T_L

L(t) = k \cdot sin(\omega_L t + \phi)

and then modulate it with a delta function array of one spike per year

s(t) = \sum\limits_{i=1}^n a_i sin(2 \pi t i +\theta_i)

This is enough to create a new aliased cycle that is simply the original frequency \omega_L summed with an infinite series of that frequency shifted by multiples of 2\pi .

f(t) = k/2 \sum\limits_{i=1}^n a_i sin((\omega_L - 2 \pi i)t +\psi_i) + ...

This was derived in a previous post. So as a concrete example, the following figure is a summed series of f(t) -- specifically what we would theoretically see for an anomalistic lunar month cycle of 27.5545 days aliased against a yearly delta.

Fig 1: Series expansion of aliasing

If you count the number of cycles in the span of 100 years, it comes out to a little less than 26 cycles, or an approximately 3.9 year aliased period. If a low-pass filter is applied to this time series, which is likely what would happen in the lagged real world, a sinusoid of period 3.9 years would emerge.

The interesting simplification is that the series above can also be expressed exactly as a biennial + odd-harmonic expansion

[cos(\pi t) + cos(3 \pi t) + cos(5 \pi t) + ...] \cdot sin(\frac{2\pi}{T} t)

where T is given by

1/T = 1/2 + int(1/T_L) - 1/T_L

where the function "int" truncates to the integer part of the period reciprocal. Since T_L is shorter than the yearly period of 1, then the reciprocal is guaranteed to be greater than one.

As a check, if we take only the first biennial term of this representation (ignoring the higher frequency harmonics) we essentially recover the same time series.

Fig 2: Biennial modulation term

What is interesting about this factoring is that a biennial modulation may naturally emerge as a result of yearly aliasing, which is possibly related to what we are seeing with the ENSO model in its biennial mode. Independent of how many lunar gravitational terms are involved, the biennial modulation would remain as an invariant multiplicative factor.

For the ENSO signal, an anomalistic term corresponding to a biennial modulation operating on a 4.085 year sinusoid appears significant in the model

1/4.085 = 1/2 + int(365.242/27.5545) - 365.242/27.5545

which is expanded as this pair of biennially split factors (see ingredient #5 in the ENSO model)

1/4.085 ~ 2/3.91 - 2/1.34

ENSO Phase Reversal

As I indicated in this post, lots of good information is available in "unpublished" NASA memorandum.

Returning to the mystery of an inferred phase reversal in the biennial forcing, the best description of the interaction of the Chandler wobble with ENSO is this NASA technical memorandum.

Fig 1: From B.F.Chao, NASA Technical Memorandum 86231

So Chao essentially reports that the Chandler wobble is in-phase with ENSO from 1900-1979, but it flips phase starting around 1980. This is precisely in line with what I am finding.

Here's an example of how you can gain confidence that you are on the right track. I configured the multiple linear regression to train on the ENSO time series applying the wave equation transform from 1880-1940, and then looked at the correlation coefficient over the interval 1940-2013 -- which is completely outside the training interval.

Fig 2: Training from 1880-1940, correlation outside. The ENSO signal is phase reversed from 1980 to 1996.

The in=inside cc is 0.725 and the out=outside cc is 0.758. It is rare that you can get a correlation higher outside of the training interval, but because of the stronger noise in the earliest ENSO measurements this is not out of the question. In terms of variance contributions, the reduced noise promotes the real underlying signal.

Next, we generate a set of several synthesized red-noise random walks that have similar variance to the actual ENSO data. This is essentially the Ornstein-Uhlenbeck algorithm, describing a random walk with a reversion-to-the-mean potential well factor. Ten runs were combined into an animated GIF shown below:

Fig 3: A set of O-U random walk profiles.

Note that although the correlation coefficient is moderate within the training interval (anywhere from 0.4 to 0.6), it is insignificant outside of the training interval, and a few times it actually goes negative. In other words, there is no phase coherence outside of the interval, and any apparent agreement within the training interval is likely the result of over-fitting.

Alas, research is still being published [1] arguing whether or not ENSO is red noise. This is actually a moot point since what the GCMs generate (as far as I can tell) are invariably reported as sets of time series with random outcomes.  If ENSO is actually deterministic -- apart from the rare metastable phase reversal and volcanic activity -- the GCMs should be slight variations of Figure 2 and not Figure 3.

[1] Chen, Xianyao, and John M. Wallace. "Orthogonal PDO and ENSO indexes." Journal of Climate 2016 (2016).

You don't need a weatherman to know which way the wind blows.