geophysics

Earthquakes, tides, and tsunami prediction

/ geoenergymath / 6 Comments

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.

q

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

qf

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

/ geoenergymath / Leave a comment

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.
Continue reading →

Second-Order Effects in the ENSO Model

/ geoenergymath / 4 Comments

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

/ geoenergymath / 2 Comments

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).

Continue reading →

The QBO anomaly of 2016 revisited

/ geoenergymath / 6 Comments

Remember the concern over the QBO anomaly/disruption during 2016?

Quite a few papers were written on the topic

  1. Newman, P. A., et al. “The anomalous change in the QBO in 2015–2016.” Geophysical Research Letters 43.16 (2016): 8791-8797.
    Newman, P. A., et al. “The Anomalous Change in the QBO in 2015-16.” AGU Fall Meeting Abstracts. 2016.
  2. Randel, W. J., and M. Park. “Anomalous QBO Behavior in 2016 Observed in Tropical Stratospheric Temperatures and Ozone.” AGU Fall Meeting Abstracts. 2016.
  3. Dunkerton, Timothy J. “The quasi‐biennial oscillation of 2015–2016: Hiccup or death spiral?.” Geophysical Research Letters 43.19 (2016).
  4. Tweedy, O., et al. “Analysis of Trace Gases Response on the Anomalous Change in the QBO in 2015-2016.” AGU Fall Meeting Abstracts. 2016.
  5. Osprey, Scott M., et al. “An unexpected disruption of the atmospheric quasi-biennial oscillation.” Science 353.6306 (2016): 1424-1427.
According to the lunar forcing model of QBO, which was also presented at AGU last year, the peak in acceleration should have occurred at the time pointed to by the BLACK downward arrow in the figure below. This was in April of this year. The GREEN is the QBO 30 hPa acceleration data and the RED is the QBO model.

Note that the training region for the model is highlighted in YELLOW and is in the interval from 1978 to 1990. This was well in the past, yet it was able to pinpoint the sharp peak 27 years later.

The disruption in 2015-2016 shown with shaded black may have been a temporary forcing stimulus.  You can see that it obviously flipped the polarity with respect to the model. This will provoke a transient response in the DiffEq solution, which will then eventually die off.

The bottom-line is that the climate scientists who pointed out the anomaly were correct in that it was indeed a disruption, but this wasn’t necessarily because they understood why it occurred — but only that it didn’t fit a past pattern. It was good observational science, and so the papers were appropriate for publishing.  However, if you look at the QBO model against the data, you will see many similar temporary disruptions in the historical record. So it was definitely not some cataclysmic event as some had suggested. I think most scientists took a less hysterical view and simply pointed out the reversal in stratospheric winds was unusual.

I like to use this next figure as an example of how this may occur (found in the comment from last year). A local hurricane will temporarily impact the tidal displacement via a sea swell. You can see that in the middle of the trace below. On both sides of this spike, the tidal model is still in phase and so the stimulus is indeed transient while the underlying forcing remains invariant. For QBO, instead of a hurricane, the disruption could be caused by a SSW event. It also could be an unaccounted-for lunar forcing pulse not captured in the model. That’s probably worth more research.

As the QBO is still on a 28 month alignment, that means that the external stimulus — as with ENSO, likely the lunar tidal force — is providing the boundary condition synchronization.

Search for El Nino

/ geoenergymath / 4 Comments

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.

Continue reading →

Variation in the Length of the Anomalistic Month

/ geoenergymath / 1 Comment

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.

The Hawkmoth Effect

/ geoenergymath / 5 Comments

Contrasting to the well-known Butterfly Effect, there is another scientific modeling limitation known as the Hawkmoth Effect.  Instead of simulation results being sensitive to initial conditions, which is the Butterfly Effect, the Hawkmoth Effect is sensitive to model structure.  It’s a more subtle argument for explaining why climate behavioral modeling is difficult to get right, and named after the hawkmoth because hawkmoths are “better camouflaged and less photogenic than butterflies”.

Not everyone agrees that this is a real effect, or it just reveals shortcomings in correctly being able to model the behavior under study. So, if you have the wrong model or wrong parameters for the model, of course it may diverge from the data rather sharply.

In the context of the ENSO model, we already provided parameters for two orthogonal intervals of the data.  Since there is some noise in the ENSO data — perfectly illustrated by the fact that SOI and NINO34 only have a correlation coefficient of 0.79 — it is difficult to determine how much of the parameter differences are due to over-fitting of that noise.

In the figure below, the middle panel shows the difference between the SOI and NINO34 data, with yellow showing where the main discrepancies or uncertainties in the true ENSO value lie. Above and below are the model fits for the earlier (1880-1950 shaded in a yellow background) and later (1950-2016) training intervals. In certain cases, a poorer model fit may be able to be ascribed to uncertainty in the ENSO measurement, such as near ~1909., ~1932, and ~1948, where the dotted red lines align with trained and/or tested model regions. The question mark at 1985 is a curiosity, as the SOI remains neutral, while the model fits to more La Nina conditions of NINO34.

There is certainly nothing related to the Butterfly Effect in any of this, since the ENSO model is not forced by initial conditions, but by the guiding influence of the lunisolar cycles. So we are left to determine how much of the slight divergence we see is due to non-stationary variation of the model parameters over time, or whether it is due to missing some other vital structural model parameters. In other words, the Hawkmoth Effect is our only concern.

In the model shown below, we employ significant over-fitting of the model parameters. The ENSO model only has two forcing parameters — the Draconic (D) and Anomalistic (A) lunar periods, but like in conventional ocean tidal analysis, to make accurate predictions many more of the nonlinear harmonics need to be considered [see Footnote 1]. So we start with A and D, and then create all combinations up to order 5, resulting in the set [ A, D, AD, A2, D2, A2D, AD2, A3, D3, A2D2, A3D, AD3, A4, D4, A2D3, A3D2, A4D1, A1D4, A5, D5 ].

This looks like it has the potential for all the negative consequence of massive over-fitting, such as fast divergence in amplitude outside the training interval, yet the results don’t show this at all.  Harmonics in general will not cause a divergence, because they remain in phase with the fundamental frequencies both inside and outside the training interval. Besides that, the higher order harmonics start having a diminished impact, so this set is apparently about right to create an excellent correlation outside the training interval.  The two other important constraints in the fit, are (1) the characteristic frequency modulation of the anomalistic period due to the synodic period (shown in the middle left inset) and (2) the calibrated lunar forcing based on LOD measurements (shown in the lower panel).

The resulting correlation of model to data is 0.75 inside the training interval (1880-1980) and 0.69 in the test interval (1980-2016).  So this gets close to the best agreement we can expect given that SOI and NINO34 only reaches 0.79.  Read this post for the structural model parameter variations for a reduced harmonic set to order 3 only.

Welcome to the stage of ENSO analysis where getting the rest of the details correct will provide only marginal benefits;  yet these are still important, since as with tidal analysis and eclipse models, the details are important for fine-tuning predictions.

Footnote:

  1. For conventional tidal analysis, hundreds of resulting terms are the norm, so that commercial tidal prediction programs allow an unlimited number of components.

 

 

 

Switching between two models

Recipe for ENSO model in one tweet

/ geoenergymath / Leave a comment

//platform.twitter.com/widgets.js


and for QBO

//platform.twitter.com/widgets.js

The common feature of the two is the application of Laplace’s tidal equation and its closed-form solution.