LOD Revisited for CSALT

One of the most questioned aspects of the CSALT model of global temperature is the LOD to Temperature factor. This creates a multi-decadal variation in temperature useful for optimizing a multiple-linear regression AGW model dependent on CO2 and other factors.

Lunisolar tides impact variations in Length-of-Day (LOD). So does ENSO and QBO. There is a recursive aspect to these relationships as well, since both LOD and ENSO have the same Chandler wobble match in apparent forcing periodicity. This is what I believe generates a 6-year signal that gets identified routinely in the LOD time-series, such as the latest finding in ref [1] below.

From ref [1], the 6-year signal in LOD series. Counting cycles this is close to an average 6.25 year period, close to the 6.4 year Chandler wobble angular momentum variation.

Its a mystery why the LOD may be considered deterministic/periodic based on how well the lunisolar tides resolve the features, but ENSO is only considered quasi-periodic or nearly chaotic (the latter according to Tsonis), even though they likely arise from common mechanisms. Above all, these phenomena all have that curious tie-in to the seasonally aliased Draconic-monthly lunar cycle.

Now we can add this paper by Marcus [2] to the mix. This is a very detailed look at the correlation between long-range LOD variations (longer than the 6-year variation of Ref [1]) and global surface temperature. His application of a 5-year running mean is essentially similar to a 5-year lag that works as a best fit in the CSALT model. Marcus stops short of assigning a source cause for the LOD-to-Temperature correlation, but the general idea is that angular momentum variations are the forcing terms that slosh the sources of heat to the surface -- "via core-induced rotational and/or related global-scale processes".  (I also have to note that Marcus is an independent researcher, who at one time had an affiliation with NASA JPL.)

From Marcus [2], correlation of LOD against various temperature indices.

All these observations of LOD, ENSO, QBO, Chandler Wobble, Flood Return periods have a strong sense of self-consistency (IMO ultimately tied to lunisolar forcing), but the problem is that the discussions are scattered among different research groups. And even on this blog, the discussions reside in scattered postings (and over at Azimuth Project, also see another lunar connection).  Eventually I will write a longer manuscript to tie it all together much like I did with The Oil Conundrum and my old fossil-fuel depletion blog.


[1] Duan, Pengshuo, Genyou Liu, Lintao Liu, Xiaogang Hu, Xiaoguang Hao, Yong Huang, Zhimin Zhang, and Binbin Wang. 2015. “Recovery of the 6-Year Signal in Length of Day and Its Long-Term Decreasing Trend.” Earth, Planets and Space 67 (1): 1.

[2] Marcus, Steven L. 2015. “Does an Intrinsic Source Generate a Shared Low-Frequency Signature in Earth’s Climate and Rotation Rate?” Earth Interactions 20 (4): 1–14. doi:10.1175/EI-D-15-0014.1.

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.

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


Inferring forced response from QBO wave equation

If the wave equation is to be effective in modeling a forced behavior such as QBO then certain invariant properties would need to be observed.

f''(t) + \omega_0^2 f(t) = F(t)

So if we split this into two parts

f''(t) \approx k_1 F(t)


\omega_0^2 f(t) \approx k_2 F(t)

I have been applying multiple linear regression to fit the wave equation via either of the above two forms

So if we combine these two and see how well the Fourier series coefficients align -- keeping care to compensate for the \omega^2 scaling factor after taking the second derivative of a sinusoid -- the linear result should be:

\frac{f''(t)}{f(t)} = \omega_0^2 k_1/k_2

and indeed ! (see Figure 1 below)

Fig 1: Correlation between direct f(t) and second-derivative f''(t) coefficient strengths. It should be linear with an intercept at the origin. The strongest term in the upper right is the aliased Draconic tidal period.  The rest are other lunar periods, solar periods, and aliased lunar harmonics.

Continue reading

Forced versus Natural Response, not a secret

This is a concept that is beaten into the head of electrical engineers.

Engineers tend to assume the forced response first, as being pessimists they always assume the worst case.  For example: Why is that annoying 60 Hz hum coming out of that circuit?  hmmm .. probably because the AC power somehow leaked into the input and the hum resulted from the forced response of the internal transfer function?  Issue solved.

Physicists tend to look for a natural response first. Being optimists, a research physicist first and foremost wants to make a discovery. A natural response would tell us something interesting about the behavior ... and perhaps a new physical law!  Look and See! The eigenvalues!   Often to them, the forced response is a nuisance that needs to be removed to observe the natural response.

But what happens when the forced response is actually the interesting physical behavior?

A prototypical example of a forced response is the phenomenon of cyclic ocean tides. The periods of tides, though complex in breadth, are merely a reflection of the orbital rates of the moon and sun.

Pictures of (left) high tide (right) low tide, linked from Science Buddies

Obviously the ocean has a natural response to a forcing stimulus, obeying some variation of the wave equation, but the overriding output directly links to the observed orbits. Even the second-order effects such as the 9-day tides are simply beats arising from nonlinear interactions of the fortnightly and monthly long period tides.  Everything from the semi-diurnal periods to the long-period monthly periods are forced responses to the gravitational tug of the lunisolar orbits.  Any resonant frequencies in the ocean's natural response are over-ridden by the forced response.

And so it likely goes for the quasi-biennial oscillations (QBO) of the equatorial stratospheric winds, as these most definitely are a forced response to the inputs ... in spite of what the current literature infers.

Continue reading