Identification of Lunar Parameters and Noise

For the ENSO model, there is an ambiguity in simultaneously identifying the lunar month duration (draconic, anomalistic, and tropical) and the duration of a year. The physical aliasing is such that the following f = frequency will give approximately equivalent fits for a range of year/month pairs (see this as well).

$f = Year/LunarMonth - 13$

So that during the fitting process, if you allow the duration of the individual months and the year to co-vary, then the two should scale approximately by the number of lunar months in a year ~13.3 = 1/0.075. And sure enough, that's what is found, a set of year/month pairs that provide a maximized fit along a ridge line of possible solutions, but only one that is ultimately correct for the average year duration over the entire range:

By regressing on the combination of linear slopes, the value of the year that minimizes the error to each of the known lunar month values is 365.244 days. This lies within the interval defined by the value of the calendar year = 365.25 days — which includes a leap day every 4 years, and the more refined leap year calculation = 365.242 days —  which includes the 100 and 400 year corrections (there are additional leap second corrections).

This analysis provides further confidence that the ENSO model is approaching the status of a metrology tool for gauging lunisolar cycles.  The tropical month is estimated slow by about 1/2 a minute, while the draconic month is fast by a 1/2 a minute, and the average anomalistic month is spot on to within a second.

This is what the fit looks like for a 365.242 day long calendar year trained over the entire interval.  It is the accumulation of the sharply matching peaks and valleys which allow the solver function to zone in so precisely to the known tidal factors.

About the only issue that hobbles our ability to achieve fits as good as ocean tidal analysis is the amount of noise near neutral ENSO conditions in the time-series data. The highlighted yellow regions in the comparison between NINO34 and SOI time-series data shown below indicate intervals whereby a sliding correlation coefficient drops closer to zero. (The only odd comparison is the blue highlighted region around 1985, where SOI is extremely neutral while NINO34 appears La Nina-like. Is SOI pressure related to a second derivative of NINO34 temperature?).

Those same yellow regions are also observed as discrepancies between the NINO34 data and the ENSO best model fit.

Yellow shading at intervals around 1930, 1936, 1948  indicate discrepancies between the NINO34 data in green and the ENSO model in red.

.

ENSO tidal forcing validated by LOD data

This is a straightforward validation of the forcing used on the lunar-driven ENSO model.

The paper by Chao et al [1] provides a comprehensive spectral analysis of the earth's length of day (LOD) variations using both a wavelet analysis and a power spectrum analysis. The wavelet analysis provides insight into the richness of the LOD cyclic variations (c.f. the Chao ref 6 in a recent post) :

Both the wavelet and the power spectrum (below) show the 6-year Fourier component that appears in the ENSO model as a mixed tidal forcing.

The original premise is that the change in LOD via the equivalent angular momentum change will impart a forcing on the Pacific ocean thermocline as per a reduced-gravity model:

Taken from presentation for ref [2]

Calculating a spectral analysis of the best fit ENSO model forcing, note that all of the model peaks (in RED) match those found by Chao et al in their ΔLOD analysis :

There are additional peaks not found by Chao but those are reduced in magnitude, as can be inferred from the log (i.e. dB) scale. If these actually exist in the Chao spectrum, they may be buried in the background noise.  Also, the missing Sa and Ssa peaks are the seasonal LOD variations that are taken into account separately by the model, as most ENSO data sets are typically filtered to remove seasonal data.

The tidal constituents shown above in the Chao power spectra are defined in the following Doodson table [3]. Chao likely is unable to discriminate the tropical values from the draconic and anomalistic values, being so close in value. On the other hand, the ENSO model needs to know these values precisely.  Each of the primary Mm, Mf, Mtm, and Mqm and satellite Msm, Msf, Mstm, Msqm factors align with the first 4 harmonics of the mixed nonlinear ENSO model with the 2nd order satellites arising from the anomalistic correction.

Tidal constituent coefficients taken from ref [3]

This is an excellent validation test because this particular LOD power spectrum has not been used previously in the ENSO model fitting process. If the peaks did not match up, then the original premise for LOD forcing would need to be reconsidered.

References

[1] B. F. Chao, W. Chung, Z. Shih, and Y. Hsieh, “Earth’s rotation variations: a wavelet analysis,” Terra Nova, vol. 26, no. 4, pp. 260–264, 2014.

[2] A. Capotondi, “El Niño–Southern Oscillation ocean dynamics: Simulation by coupled general circulation models,” Climate Dynamics: Why Does Climate Vary?, pp. 105–122, 2013.

[3] D. D. McCarthy (ed.): IERS Conventions (1996) (IERS Technical Note No. 21) :
Chapter 6

Earthquakes, tides, and tsunami prediction

I've been wanting to try this for awhile — to see if the solver setup used for fitting to ENSO would work for conventional tidal analysis.  The following post demonstrates that if you give it the recommended tidal parameters and let the solver will grind away, it will eventually find the best fitting amplitudes and phases for each parameter.

The context for this analysis is an excellent survey paper on tsunami detection and how it relates to tidal prediction:

S. Consoli, D. R. Recupero, and V. Zavarella, “A survey on tidal analysis and forecasting methods for Tsunami detection,” J. Tsunami Soc. Int.
33 (1), 1–56.

The question the survey paper addresses is whether we can one use knowledge of tides to deconvolute and isolate a tsunami signal from the underlying tidal sea-level-height (SLH) signal. The practical application paper cites the survey paper above:

Percival, Donald B., et al. "Detiding DART® buoy data for real-time extraction of source coefficients for operational tsunami forecasting." Pure and Applied Geophysics 172.6 (2015): 1653-1678.

This is what the raw buoy SLH signal looks like, with the tsunami impulse shown at the end as a bolded line:

After removing the tidal signals with various approaches described by Consoli et al, the isolated tsunami impulse response (due to the 2011 Tohoku earthquake) appears as:

As noted in the caption, the simplest harmonic analysis was done with 6 constituent tidal parameters.

As a comparison, the ENSO solver was loaded with the same tidal waveform (after digitizing the plot) along with 4 major tidal parameters and 4 minor parameters to be optimized. The solver's goal was to maximize the correlation coefficient between the model and the tidal data.

The yellow region is training which reached almost a 0.99 correlation coefficient, with the validation region to the right reaching 0.92.

This is the complex Fourier spectrum (which is much less busy than the ENSO spectra):

The set of constituent coefficients we use is from the Wikipedia page where we need the periods only. Of the following 5 principal tidal constituents, only N2 is a minor factor in this case study.

In practice, multiple linear regression would provide faster results for tidal analysis as the constituents add linearly (see the CSALT model). In contrast, for ENSO there are several nonlinear steps required that precludes a quick regression solution.  Yet, this tidal analysis test shows how effective and precise a solution the solver supplies.

The entire analysis only took an evening of work, in comparison to the difficulty of dealing with ENSO data, which is much more noisy than the clean tidal data.  Moreover, the parameters for conventional tidal analysis stress the synodic/tropical/sidereal constituents — unfortunately, these are of little consequence for ENSO analysis, which requires the draconic and anomalistic parameters and the essential correction factors. The synodic tide is the red herring for the unwary when dealing with global phenomena such as ENSO, QBO, and LOD. The best way to think about it is that the synodic cycle impacts locally and most immediately, whereas the anomalistic and draconic cycles have global and more cumulative impacts.

The 6-year oscillation in Length-of-Day

A somewhat hidden cyclic variation in the length-of-day (LOD) in the earth's rotation, of between 6 and 7 years, was first reported in Ref [1] and analyzed in Ref [2]. Later studies further refined this period [3,4,5] closer to 6 years.

 Change in detected LOD follows a ~6-yr cycle, from Ref [3]

It's well known that the moon's gravitational pull contributes to changes in LOD [6]. Here is the set of lunar cycles that are applied as a forcing to the ENSO model using LOD as calibration.

Limits to Goodness of Fit

Based on a comparison of local interval correlations between the NINO34 and SOI indices, there probably is a limit to how well a model can be fit to ENSO.  The lower chart displays a 4-year-windowed correlation coefficient (in RED) between the two indices (shown in upper chart):

Note that in the interval starting at 1930, the correlation is poor for about 7 years.

Next note that the ENSO model fit shows a poor correlation to the NINO34 data in nearly the same intervals (shown as dotted GREEN). This is an odd situation but potentially revealing. The fact that both the ENSO model and SOI don't match the NINO34 index over the same intervals, suggests that the model may match SOI better than it does NINO34.  Yet, because of the excessive noise in SOI, this is difficult to verify.

But more fundamentally, why would NINO34 not match SOI in these particular intervals? These regions do seem to be ENSO-neutral, not close to El Nino or La Nina episodes.  Some also seem to occupy regions of faster, noisy fluctuations in the index.

It could be that the ENSO lunar tidal model is revealing the true nature of the ENSO dynamics, and these noisier, neutral regions are reflecting some other behavior (such as amplitude folding) — but since they also appear to be obscured by noise, it makes it difficult to unearth.

The paper by Zajączkowska[1] also applies a local correlation to compare the lunar tidal cycles to plant growth dynamics. There's a treasure trove of recent research on this topic.

Second-Order Effects in the ENSO Model

For ocean tidal predictions, once an agreement is reached on the essential lunisolar terms, then the second-order terms are refined. Early in the last century Doodson catalogued most of these terms:

"Since the mid-twentieth century further analysis has generated many more terms than Doodson's 388. About 62 constituents are of sufficient size to be considered for possible use in marine tide prediction, but sometimes many fewer can predict tides to useful accuracy."

That's possibly the stage we have reached in the ENSO model.  There are two primary terms for lunar forcing (the Draconic and Anomalistic) cycles, that when mixed with the annual and biannual cycles, will reproduce the essential ENSO behavior.  The second-order effects are the  modulation of these two lunar cycles with the Tropical/Synodic cycle.  This is most apparent in the modification of the Anomalistic cycle. Although not as important as in the calculation of the Total Solar Eclipse times, the perturbation is critical to validating the ENSO model and to eventually using it to make predictions.

The variation in the Anomalistic period is described at the NASA Goddard eclipse page. They provide two views of the variation, a time-domain view and a histogram view.

 Time domain view Histogram view

Since NASA Goddard doesn't provide an analytical form for this variation, we can see if the ENSO Model solver can effectively match it via a best-fit search to the ENSO data. This is truly an indirect method.

First we start with a parametric approximation to the variation, described by a pair of successive frequency modulated (and full-wave rectified) terms that incorporate the Tropical-modified term, wm. The Anomalistic term is wa.

COS(wa*t+pa+c_1*ABS(SIN(wm*t+k_1*ABS(SIN(wm*t+k_2))+c_2)))

$\cos(\omega_a t+\phi_a+c_1 \cdot |\sin(\omega_m t+k_1 \cdot |\sin(\omega_m t+k_2)|+c_2)|)$

This can generate the cusped behavior observed, but the terms pa, c_1, c_2, k_1, and k_2 need to be adjusted to align to the NASA model. The solver will try to do this indirectly by fitting to the 1880-1950 ENSO interval.

Plotting in RED the Anomalistic time series and the histogram of frequencies embedded in the ENSO waveform, we get:

 Time domain view of model Histogram view of model

This captures the histogram view quite well, and the time-domain view roughly (in other cases it gives a better cusped fit).  The histogram view is arguably more important as it describes the frequency variation over a much wider interval than the 3-year interval shown.

What would be even more effective is to find the correct analytical representation of the Anomalistic frequency variation and then plug that directly into the ENSO model. That would provide another constraint to the solver, as it wouldn't need to spend time optimizing for a known effect.

Yet as a validation step, the fact that the solver detects the shape required to match the variation is remarkable in itself. The solver is obviously searching for the forcing needed to produce the ENSO waveform observed, and happens to use the precise parameters that also describe the second-order Anomalistic behavior.  That could happen by accident but in that case there have been too many happy accidents already, i.e. period match, LOD match, Eclipse match, QBO match, etc.

Using Solar Eclipses to calibrate the ENSO Model

This is the forcing for the ENSO model, focusing on the non-mixed Draconic and Anomalistic cycles:

Note that the maximum excursions (perigee and declination excursion) align with the occurrence of total solar eclipses. These are the first three that I looked at, which includes the latest August 21 eclipse in the center chart.

There are about 90 more of these stretching back to 1880. The best way to fit the calibration is to take the negative excursions of the two lunar forcings and multiply these together, i.e. use the effective Draconic*Anomalistic amplitudes (also only take the fortnightly cycle of the Draconic, as eclipses occur during both the ascending and descending node crossings). The main fitting factors are the phases of the two lunar months.  To get the maximum alignment from the search solver, we maximize the sum of the effective amplitudes across the entire interval. This results in a phase difference between the two of about 0.74 radians based at the starting year of 1880 (i.e. year 0).

Search for El Nino

The model for ENSO includes a nonlinear search feature that finds the best-fit tidal forcing parameters.  This is similar to what a conventional ocean tidal analysis program performs — finding the best-fitting lunar tidal parameters based on a measured historic interval of hundreds of cycles. Since tidal cycles are abundant — occurring at least once per day — it doesn't take much data collected over a course of time to do an analysis.  In contrast, the ENSO model cycles over the course of years, so we have to use as much data as we can, yet still allow test intervals.

What follows is the recipe (more involved than the short recipe) that will guarantee a deterministic best-fit from a clean slate each time. Very little initial condition information is needed to start with, so that the final result can be confidently recovered each time, independent of training interval.

Variation in the Length of the Anomalistic Month

For the ENSO model, we use two constraints for the fitting process. One of the constraints is to maximize the correlation coefficient for the model over the ENSO training interval selected. The other constraint is to maximize the correlation of the selected lunar tidal forces over a measured Length-of-Day (LOD) interval. The latter constrains the lunar tidal forcing to known values that will actually change the angular momentum of the earth's rotation. This in turn drives the sloshing of the Pacific ocean's thermocline leading to the ENSO cycle. The two constraints are simultaneously met by heuristically maximizing the average of the correlation coefficients.

In addition, there are the fixed constraints of the primary lunar periods corresponding to the Draconic/nodal cycle and the Anomalistic cycle.

This combination gives a fairly effective fit over the entire training cycle, but there is an important additional constraint that needs to be applied to the Anomalistic cycle. The NASA eclispse and moon's orbit page describes the situation :

"The anomalistic month is defined as the revolution of the Moon around its elliptical orbit as measured from perigee to perigee. The length of this period can vary by several days from its mean value of 27.55455 days (27d 13h 18m 33s). Figure 4-4 plots the difference of the anomalistic month from the mean value for the 3-year interval 2008 through 2010. ... the eccentricity reaches a maximum when the major axis of the lunar orbit is pointed directly towards or directly away from the Sun (angles of 0° and 180°, respectively). This occurs at a mean interval of 205.9 days, which is somewhat longer than half a year because of the eastward shift of the major axis. "

This is a significant variation in the anomalistic cycle over the course of a year. We don't use this variation as a constraint but we can use it as a fitting parameter and then compare the variation obtained over that shown above.

Using the 205.9 day value ~365/(2-2/8.85), we break this into Fourier components of half this value and twice this value. The mean Anomalistic period of 27.5545 days is then frequency modulated by the slower periods by the standard engineering procedure. We then allow the amplitude and phase of each factor to vary during the training to obtain the best fit (this is slightly different than the concise form used previously).

If we zoom in on the anomalistic period variation, we get this match to the NASA Goddard model:

There is no reason to believe that this match would spontaneously occur given that there are 3 amplitudes and 3 phase factors involved. Yet it matches precisely to the (1) peak positions, (2) relative amplitudes, and to the (3) cusped shape via the Fourier series summation. An even better fit is obtained if we use abs(sin(π time/205.9+Φ)) as the fitting function as it naturally creates more of a cusp shape due to the full-wave rectification of the sine wave.

Conventional tidal analysis is renowned for being an exacting procedure [1], where the known tidal periods are broken down into equivalently similar sets of harmonic factors, yet applied on a diurnal or semidurnal basis. The only difference here is that ENSO responds to the monthly and fortnightly long-period tides and not the short-period ones.

→ This model fit gives further validation to the lunar tidal mechanism for forcing ENSO.  The exacting process of generating the correct lunar tidal variations (along with the subtle biennial modulation and the tricky aliasing) have likely contributed to the fact that the pattern has remained hidden for so long.  This is actually not so different a situation as the long hidden connection between triggering of earthquakes and the dynamic  moon-sun-earth alignment. That pattern is also hidden, only exposed recently. Alas, not everything can be quite as obvious as the pattern matching of ocean surface tides to the lunisolar cycles.

References

[1] S. Consoli, D. R. Recupero, and V. Zavarella, “A survey on tidal analysis and forecasting methods for Tsunami detection,” arXiv preprint arXiv:1403.0135, 2014.

UPDATE

The reason for the peculiar shape of the Anomalistic frequency variations is due to a different slope (i.e. velocity) on one lobe of the elliptical orbit than on the other. You can get this by generating another sinusoidal modulation on top of the average elliptical sinusoid. This generates an asymmetric sawtooth in the phase angle (see blue line below) and the characteristic spiked or cusped profile in the effective frequency or derivative of this value (see red line below).
I am starting to use this formulation in the ENSO model as it is quite concise.