# Solver vs Multiple Linear Regression for ENSO

In the previous ENSO post I referenced the Rajchenbach article on Faraday waves.

There is a telling assertion within that article:

"For instance, to the best of our knowledge, the dispersion relation (relating angular frequency ω and wavenumber k) of parametrically forced water waves has astonishingly not been explicitly established hitherto. Indeed, this relation is often improperly identified with that of free unforced surface waves, despite experimental evidence showing significant deviations"

What they are suggesting is that too much focus has been placed on natural resonances and the dispersion relationships within a free fluid volume. Whereas the forced response is clearly as important — if not more — and that the forcing will show through in the solution of the equations. I have been pursuing this strategy for a while, having started down the Mathieu equation right away and then eventually realizing the importance of the forced response, yet the Rajchenbach article is the first case that I have found made of what I always thought should be a rather obvious assumption. The fact that the peer-reviewers allowed the "astonishingly" adjective in the paper is what makes it telling. It's astonishing in the equivalent sense that Rajchenbach & Clamond are pointing out that a pendulum's motion will be impacted by a periodic push. In other words, astonishing in the sense that this premise should be obvious!

# Is ENSO predominantly tidal?

The model of ENSO is split into two groups of forcing factors, unified by a biennial modulation. The predominant factors have cycles of 6.4 and 14 years, which extends back through historical proxy data.  Other strong factors track the same lunar cycles that appear as QBO forcing factors.

What's somewhat interesting is that after a multiple linear regression, the optimal values are 6.476 and 13.98 years -- with the mean frequency of these two values at 8.852 years, which is very close to the anomalistic tidal long-period of 8.85 years.

# Modeling Red Noise versus ENSO

With respect to the ENSO model I have been thinking about ways of evaluating the statistical significance of the fit to the data. If we train on one 70 year interval and then test on the following 70 year interval, we get the interesting effect of finding a higher correlation coefficient on the test interval. The training interval is just below 0.85 while the test is above 0.86.

This image has been resized to fit in the page. Click to enlarge.

# ENSO model maps to LOD cycles

Elaborating on this comment attached to an LOD post,  noting this recent paper:

Shen, Wenbin, and Cunchao Peng. 2016. “Detection of Different-Time-Scale Signals in the Length of Day Variation Based on EEMD Analysis Technique.” Special Issue: Geodetic and Geophysical Observations and Applications and Implications 7 (3): 180–86.  doi:10.1016/j.geog.2016.05.002.

Because of the law of conservation of momentum sloshing can change the velocity of a container full of liquid, momentarily speeding it up or slowing it down as the liquid sloshes back and forth.  By the same token, suddenly slowing or speeding of that container can also cause the sloshing.   So there is a chicken and egg quality to the analysis of sloshing, making it difficult to ascertain the origin of the effect.

If ENSO is a manifestation of a liquid sloshing in a container and if the length-of-day (LOD) is a measurement of the angular momentum changes of the Earth's rotation, then it is perhaps useful to compare the fundamental time-varying signals in each measurement.

# ENSO Model Final Stretch (maybe)

I recently posted a bog article called QBO Model Final Stretch. The idea with that post was to give an indication that the physics and analytical math model explaining the behavior of the QBO was in decent shape. I would like to do the same thing with the ENSO model but retain the context of the QBO model.  Understanding the QBO was a boon to making progress with ENSO as it provided an excellent training ground for time-series analysis and also provided some insight into the underlying forcing functions.  In the literature, there is a clear indication that ENSO and QBO are somehow related, but the causality chain remains unclear.

# Obscure paper on ENSO determinism

Glenn Brier is on my list of essential climate science researchers. I am sure he has since retired but his status as both a fellow of the American Meteorological Society and a fellow of the American Statistical Association indicated he knew his stuff. One of his research topics was exploring periodicities in climate data. He was eminently qualified for this kind of analysis as his statistics background provided him with the knowledge and skill in being able to distinguish between stochastic versus deterministic behaviors.

One paper Brier co-authored in 1989 is an obscure gem in terms of understanding the determinism of El Nino and ENSO [1]. He essentially analyzed a 463-year historical chronology of strong El Nino events and was able to find significant periodicities of 6.75 and 14 years in the data.

This compares very well with what I have  been able to decipher as the strongest ENSO periodicities (taken from the instrumental record since 1880) of 6.5 and 14 years.

# Alternate simplification of QBO from Laplace's Tidal Equations

Here is an alternate derivation of reducing Laplace's tidal equations along the equator. Recall that the behavior of QBO (and ENSO) is isolated to the equator.

For a fluid sheet of average thickness $D$, the vertical tidal elevation $\zeta$, as well as the horizontal velocity components $u$ and $v$ (in the latitude $\varphi$ and longitude $\lambda$ directions).

This is the set of Laplace's tidal equations (wikipedia). Along the equator, for $\varphi$ at zero we can reduce this.

{\begin{aligned}{\frac {\partial \zeta }{\partial t}}+{\frac {1}{a\cos(\varphi )}}\left[{\frac {\partial }{\partial \lambda }}(uD)+{\frac {\partial }{\partial \varphi }}\left(vD\cos(\varphi )\right)\right]=0,\\{\frac {\partial u}{\partial t}}-v\left(2\Omega \sin(\varphi )\right)+{\frac {1}{a\cos(\varphi )}}{\frac {\partial }{\partial \lambda }}\left(g\zeta +U\right)=0,\\{\frac {\partial v}{\partial t}}+u\left(2\Omega \sin(\varphi )\right)+{\frac {1}{a}}{\frac {\partial }{\partial \varphi }}\left(g\zeta +U\right)=0,\end{aligned}}

where $\Omega$ is the angular frequency of the planet's rotation, $g$ is the planet's gravitational acceleration at the mean ocean surface, $a$ is the planetary radius, and $U$ is the external gravitational tidal-forcing potential.

The main candidates for removal due to the small-angle approximation along the equator are the second terms in the second and third equation. The plan is to then substitute the isolated $u$ and $v$ terms into the first equation, after taking another derivative of that equation with respect to $t$.

{\begin{aligned}{\frac {\partial \zeta }{\partial t}}+{\frac {1}{a}}\left[{\frac {\partial }{\partial \lambda }}(uD)+{\frac {\partial }{\partial \varphi }}\left(vD\right)\right]=0,\\{\frac {\partial u}{\partial t}}+{\frac {1}{a}}{\frac {\partial }{\partial \lambda }}\left(g\zeta +U\right)=0,\\{\frac {\partial v}{\partial t}}+{\frac {1}{a}}{\frac {\partial }{\partial \varphi }}\left(g\zeta +U\right)=0,\end{aligned}}

Taking another derivative of the first equation:

{\begin{aligned}{a\frac {\partial^2 \zeta }{\partial t^2}}&+ \frac {\partial }{\partial t} \left[{\frac {\partial }{\partial \lambda }}(uD)+{\frac {\partial }{\partial \varphi }}\left(vD\right)\right]=0,\end{aligned}}

Next, on the bracketed pair we invert the order of derivatives (and pull out $D$, which may not be right for ENSO where the thermocline depth varies!)

{\begin{aligned}{a\frac {\partial^2 \zeta }{\partial t^2}}&+ D \left[{\frac {\partial }{\partial \lambda }}( \frac {\partial u }{\partial t} )+{\frac {\partial }{\partial \varphi }}(\frac {\partial v }{\partial t} )\right]=0,\end{aligned}}

Notice now that the bracketed terms can be replaced by the 2nd and 3rd of Laplace's equations

{\begin{aligned}{\frac {\partial u}{\partial t}}&=-{\frac {1}{a}}{\frac {\partial }{\partial \lambda }}\left(g\zeta +U\right)\end{aligned}}

{\begin{aligned}{\frac {\partial v}{\partial t}}&=-{\frac {1}{a}}{\frac {\partial }{\partial \varphi }}\left(g\zeta +U\right)\end{aligned}}

as

{\begin{aligned}{a^2\frac {\partial^2 \zeta }{\partial t^2}}&- D \left[{\frac {\partial }{\partial \lambda }}( {{\frac {\partial }{\partial \lambda }}\left(g\zeta +U\right)} )+{\frac {\partial }{\partial \varphi }}({{\frac {\partial }{\partial \varphi }}\left(g\zeta +U\right)})\right]=0\end{aligned}}

The $\lambda$ terms are in longitude so that we can use a separation of variables approach and create a spatial standing wave for QBO, $SW(s)$ where $s$ is a wavenumber.

{\begin{aligned}{\frac {\partial^2 \zeta }{\partial t^2}}&- D \left[( SW(s) \zeta )+{\frac {\partial }{\partial \varphi }}({{\frac {\partial }{\partial \varphi }}\left(g\zeta +U\right)})\right]=0\end{aligned}}

The next bit is the connection between a change in latitudinal forcing with a temporal change:

\begin{aligned} {\frac {\partial \zeta }{\partial \varphi } = \frac {\partial \zeta }{\partial t} \frac {\partial t }{\partial \varphi } } \end{aligned}

so if

\begin{aligned} \frac {\partial \varphi }{\partial t } = \sum_{i=1}^{i=N} k_i \omega_i \cos(\omega_i t) \end{aligned}

to describe the external gravitational forcing terms, then the solution is the following:

\begin{aligned} \zeta(t) = \sin( \sqrt{A} \sum_{i=1}^{i=N} k_i \sin(\omega_i t) ) \end{aligned}

where $A$ is an aggregate of the constants of the differential equation.

So we can eventually get to a fit that looks like the actual QBO, for a fit from 1954 to 2015 for the 30 hPa data.

Besides the excursion-limiting behavior imposed by a $sin(sin(t))$ modulation, this formulation can also show amplitude folding at the positive and negative extremes. In other words, if the amplitude is too large, the outer $sin$ modulation starts to shrink the excursion, instead of just limiting it. Note around the 2016 time frame, the negative peak starts folding inward.

One caveat in this derivation is that the QBO may actually be the $v(t)$ term - the horizontal longitudinal velocity of the fluid, the wind in other words - which can be derived from the above by applying the solution to Laplace's third tidal equation in simplified form above.

EDIT (8/22/2016) — The caveat in the paragraph above is resolved in this post.

What I have derived is a first step, good enough for modeling the QBO for the time being. But the dynamics of ENSO substantially differ in character from the QBO, showing more extreme excursions. How could that come about?

Consider if instead of $U$ being a constant, it varies in a similar sinusoidal fashion, then that term will go to the RHS as a forcing, and the above impulse response function will need to be convolved with that forcing . The same goes for the $D$ factor which I foreshadowed earlier. Having these vary with time changes the character of the solution. So instead of having a relatively constant amplitude envelope characteristic of QBO, the waveform amplitude will become more erratic, and more similar to the dynamics of ENSO, where the peak excursions vary a considerable amount.

Modifying the $D$ factor in fact makes the solution closer to a Mathieu equation formulation. The stratosphere is constant in depth (i.e. stratified) so that's why we set that to a constant. Yet remember that with ENSO, it is the thermocline depth that is sloshing wildly. The density differences of water above and below the thermocline is what is sensitive to the lunar gravitational forcing. That's why with a strict biennial modulation applied to the $D$ term allows an ENSO model that shows a behavior that matches the data so well:

This is a fit for ENSO from 1880 to 1950, with the validation to the right

Note how well the amplitudes vary in prediction.

In contrast for QBO, a training run from 1953 to 1983 shows this (this may be sign inverted, I don't care since the positive direction of QBO is defined arbitrarily):

The validated model to the right of 1983 shows much less variation with amplitude, in keeping with the excursion-limiting behavior of the sin-of-a-sin formulation (i.e. $sin(sin(t))$).

These kinds of approximations are very common in the physics that I was taught, especially in the hairy world of solid-state physics. Simplifying the starting set of equations will not only reduce complexity but it often brings a level of understanding not normally accessible from looking at the original formulation.

# 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:  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)

## References

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

# 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

# Deterministically Locked on the ENSO Model

After several detours and dead-ends, it looks as if I have locked on a plausible ENSO model, parsimonious with recent research.  The sticky widget almost from day 1 was the odd behavior in the ENSO time-series that occurred starting around 1980 and lasting for 16 years. This has turned out to be a good news/bad news/good news opportunity.  Good in the fact that applying a phase inversion during that interval allowed a continuous fit across the entire span.  Yet it's bad in that it gives the impression of applying a patch, without a full justification for that patch. But then again good in that it explains why other researchers never found the deterministic behavior underlying ENSO -- applying conventional tools such as the Fourier transform aren't much help in isolating the phase shift (accepting Astudillo's approach).

Having success with the QBO model, I wasn't completely satisfied with the ENSO model as it stood at the beginning of this year. It didn't quite click in place like the QBO model did. I had been attributing my difficulties to the greater amount of noise in the ENSO data, but I realize now that's just a convenient excuse -- the signal is still in there.  By battling through the inversion interval issue, the  model has improved significantly.  And once the correct forcing and Mathieu modulation is applied, the model locks in place to the data with the potential to work as well as a deterministic tidal prediction algorithm.