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.

Continue reading

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.

Continue reading

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}}


{\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:

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

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.

Continue reading

Crucial recent citations for ENSO

I have found the following research articles vital to formulating a basic model for ENSO.

The first citation finds the disturbance after 1980 leading to the identification of a phase reversal in the ENSO behavior. They apply Takens embedding theorem (which works for linear and non-linear systems such as Mathieu and Hill) to the time series, reconstructing current and future behavior from past behavior.

H. Astudillo, R. Abarca-del-Rio, and F. Borotto, “Long-term non-linear predictability of ENSO events over the 20th century,” arXiv preprint arXiv:1506.04066, 2015.

Continue reading

Tidal Locking and ENSO

As I have been formulating a model for ENSO, I always try to relate it to a purely physical basis. The premise I have had from the beginning is that some external factor is driving the forcing of the equatorial Pacific thermocline. This forcing stimulus essentially causes a sloshing in the ocean volume due to small changes in the angular momentum of the rotating earth. I keep thinking that the origin is lunar as the success of the QBO model in relating lunisolar forcing to the oscillatory behavior of the QBO winds is enough motivation to keep on a lunar path.

Yet, I am finding that the detailed mechanism for ENSO differs from that of QBO. An interesting correlation I found is in the tidal-locking of the Earth to the moon. I think this is a subset of the more general case of spin-orbit resonance, where the rotation rate of a satellite is an integral ratio of the main body. In the case of the moon and the earth, it explains why the same moon face is always directed at the earth -- as they spin at the same rate during their mutual orbit, thus compensating via a kind of counter-rotation as shown in the left figure below.

Fig 1: Tidal Locking (left) results in the spin of the moon

Continue reading

Biennial Mode of SST and ENSO

This recent 2016 paper [1] by Kim is supporting consensus to my model of a modulated biennial forcing to ENSO. I had read some of Kim's earlier papers [2] where he introduced the idea of cyclostationary behavior.

The insight that they and I share is that the strictly biennial oscillation is modulated by longer frequencies such that +/- sideband frequencies are created around the 2-year period.

sin(\pi t) \cdot sin(\omega_m t) = \frac{1}{2} ( cos(\pi t - \omega_m t) - cos(\pi t + \omega_m t) )

This is not aliasing but essentially a non-aliased frequency modulation of the base cycle. The insight is clarified by Kim with respect to Meehl's [3] tropospheric biennial oscillation (TBO).

From [1]

That biennial mode locks the base frequency in place to a seasonal cycle, with the modulation creating what looks like a more chaotic pattern. That's why ENSO has been so stubborn to analysis, in that the number of Fourier frequencies doubles with each modulating term. Yet in reality it's likely half as complicated as most scientists have been lead to believe.

Continue reading

Validating ENSO cyclostationary deterministic behavior

I tend to write a more thorough analysis of research results, but this one is too interesting not to archive in real-time.

First, recall that the behavior of ENSO is a cyclostationary yet metastable standing-wave process, that is forced primarily by angular momentum changes. That describes essentially the physics of liquid sloshing. Setting input forcings to the periods corresponding to the known angular momentum changes from the Chandler wobble and the long-period lunisolar cycles, it appears trivial to capture the seeming quasi-periodic nature of ENSO effectively.

The key to this is identifying the strictly biennial yet metastable modulation that underlies the forcing. The biennial factor arises from the period doubling of the seasonal cycle, and since the biennial alignment (even versus odd years)  is arbitrary, the process is by nature metastable (not ergodic in the strictest sense).  By identifying where a biennial phase reversal occurs, the truly cyclostationary arguments can be isolated.

The results below demonstrate multiple regression training on 30 year intervals, applying only known factors of the Chandler and lunisolar forcing (no filtering applied to the ENSO data, an average of NINO3.4 and SOI indices). The 30-year interval slides across the 1880-2013 time series in 10-year steps, while the out-of-band  fit maintains a significant amount of coherence with the data:

Continue reading