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 , 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: From , 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
If while not a function of , then the operator can be reduced to:
The next crucial step is to extract the implicit temporal change from via the chain rule. Even though is small, changes in have to be retained to first-order or else the solution would be static along the equator.
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 , then
With that simplification, we reduce the formulation to a temporal second-order differential equation.
where f(t) is a forcing added to the right-hand side (RHS).
Moving the from the denominator, we recover the Mathieu-like wave equation formulation:
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 , 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)
 Wang, Houjun, John P. Boyd, and Rashid A. Akmaev. "On computation of Hough functions." Geoscientific Model Development 9.4 (2016): 1477. PDF