[mathjax]Is there such a thing as too simple a model? Take a look at this fit to ENSO
Fig 1: The SOI fit explained in this post. The time axis is in months from 1880.
My original approach to reconstructing the SOI time series using basis Mathieu functions is useful but not the only way to do the modeling. I believe that I took that as far as I could go and so decided to bite the bullet and solve the SOIM (Southern Oscillation Index Model) differential equations numerically. What pops out is a gloriously simple climate model that fits the data remarkably well — and BTW, one of those models deemed not to exist by the climate science deniers.
What follows is the description of the model and how it connects to the previously described Chandler Wobble along with a crucial link to the quasi-biennial oscillation (QBO). Some type of forcing is responsible for providing the necessary energy to get the ENSO sloshing behavior going, and these are as good of candidates as any.
Read on to see how it can’t possible be done 🙂
The SOIM differential equation that we want to solve numerically was introduced here. If we arrange it so that the left hand side is the perturbed wave equation, and the right-hand side is a forcing function, F(t), this gives us a non-linear diffEq (shorthand for differential equation) that will require numerical techniques to solve for anything more sophisticated than a delta impulse forcing stimulus, F(t) = δ(t).
$$ frac{d^2x(t)}{dt^2}+[a-2qcos(2omega t)]x(t)=F(t) $$
This equation also describes a known variation of sloshing behavior in a fluid-filled tank [1]. I will at some later point reproduce the derivation, but for now read Frandsen’s thorough treatise, and in particular pay attention to how the forcings, vertical and horizontal,are introduced. (Note that I really don’t care that the math has only been previously applied to tanks built on a human scale. Just because we are dealing with a “tank” the size of the Pacific Ocean doesn’t change the underlying math and physics)
In trying to solve this equation we are definitely entering the computational physics realm, where industrial strength computer science techniques are required for finding solutions [2]. What I lean on right now is a combination of the NonLinearModelFit and ParametricNDSolve capability available in Mathematica. After that, I use the Eureqa tool exploratory model explanation tool to identify that my nonlinear equation is on the Pareto frontier of plausible solutions.
The basic premise is that we generate a forcing function that is a combination of a Chandler Wobble-related inertial mechanism with radial frequency W and a QBO-related forcing of frequency Q. We know that the Chandler Wobble will be around a 6 year period and that the QBO frequency will be around 2.3 to 2.4 years. The periodic perturbation in the wave equation, ω, is also assumed to be approximately 6 years based on the earlier analysis assuming an impulse response only. So we have this equation that we want to solve numerically
$$ frac{d^2x(t)}{dt^2}+[a-2qcos(2omega t)]x(t)= CW sin(Qt+psi) + QB sin(Wt + phi) $$
The initial conditions of x(0) and x'(0) are also included but are weaker terms than the long-term cyclic forcings. This has to do with the fact that climate is based on boundary conditions whereas weather is based on initial conditions as described by Andrew Lacis [3] . In this case, the boundary conditions are set by the ENSO forcing functions, which will synchronize over time, thereby over-riding the initial conditions of any point in time.
Characterization of the SOI using the SOIM differential equation
The initial temptation is to apply this equation to the entire SOI time series and see if we can discover any kind of correlation at all. That will work, but we want to reduce the state search space beforehand. So the first tactic is to model the QBO dataset. The QBO data starting from the year 1953 is available here. Figure 2 below shows the result of a best fit to the QBO data using a few Fourier series terms. We will use this as an initial seed to the solver and then later come back to verify any adjustments made.
Fig 2: The QBO data starting in 1953 was approximated by a few terms in the Fourier series and a constant term that is factored out. Note that this is likely a varying term of around 2.4 years and a weaker term at 6+ years.
The next step is to generate a trial equation and zero in on the Mathieu equation modulation by hand (painful but instructive) or by an automated search mechanism (painless but risky if it does not converge or gets stuck in a local minimum). As a first try, we opt for the manual search mechanism.
Using the QBO model fit above as a seed, we start from 1932 and create a time series model that runs up to 2005. The fit was to the SOI data maintained by NCAR here. Since zero crossings are important in characterizing oscillatory behavior and long-term oscillations are a nuisance, the Tahiti and Darwin sets have the multidecadal variations removed via a LOESS filter (this is critical for a cross check performed later)
The result is shown in Figure 3 along with the Mathematica code used to generate the model. The fit evaluation was done by visual inspection and maximizing the correlation coefficient (where absolute scaling is not a vital factor)
Fig 3 : The SOIM fit to the SOI data starting in 1932 and ending in 2005. Most peaks and valleys are reproduced.
This converges remarkably quickly to an acceptable correlation coefficient, largely because the initial QBO seed function was effective in narrowing the parameter search space.
Using this seed, the data pre-1932 was then fit to the same SOIM differential equation using the post-1932 parameters as seeding.
Fig 4 : The SOIM fit covering the pre-1932 SOI time series.
Again the fit is remarkable considering the model’s simplicity, with only slight changes needed in the SOIM parameters.
Next we can reconstruct the modeled QBO that was parameter tweaked against the original QBO data. This is shown below in Figure 5. Since only a few periods were used to approximate the slightly varying QBO time series, the fit would not be expected to be perfect. Yet it does align along long stretches where the data and model are in phase, and is only poor in the regions highlighted in yellow. Unsurprisingly, these yellow regions correspond to the intervals where the DiffEq solution in Figure 3 also gives a poor SOIM fit with respect to the SOI data. For example, the interval starting in the late 1960’s shows a phase reversal artifact in the QBO model and this corresponds to a bad fit in the same interval starting in month #1050 = late 1960’s for the SOI model. Also note the interval starting in 1990. We can expect these disagreements to exist, as the initial model was purposely made to be as simple as possible so that we could not be accused of over-fitting.
Fig 5: The backfit of the optimized QBO model (in dashed RED) giving the best possible SOI agreement against the original QBO data (in BLUE). Due to the simple approximation of a slowly varying QBO wavelength, beat effects will introduce artifacts that reduce the model’s fidelity in certain regions (shown as yellow highlights).
Next, we reconstruct the Mathieu equation modulation over the entire interval. In a previous post, using the Mathieu basis functions alone, we identified coherence with a 6+ year cycle — this was postulated to align with the Chandler Wobble beat frequency. But as Figure 6 below shows, the numerically computed model appears as an 8+ year cycle, with a subtle period change prior to the year 1932. This may not contradict the previous findings, as the QBO period of ~2.4 years interacting with an 8+ year cycle could create a strong difference harmonic at ~3.3 years — which is the Mathieu period that will double due to the nonlinearity to a physical effect at 6.6 years. Also, in a previous post I found a Mathieu basis function with a semidiurnal tidal period of 8.85 years that showed some promise in fitting the SOI data. So this could also be a mis-identification of the observed frequencies.
Fig 6 : The Mathieu equation modulation term is a steady oscillation over time, with a perceptible change occurring around 1932.
The Chandler wobble also showed significant effects in the years from 1920 to 1930, which is the interval at which the Mathieu equation modulation changes perceptibly and we lose coherence over longer time scales. This loss of phase coherence is shown in the figure below
Fig 7: This figure illustrates what happens if an early period fit is extrapolated to later years without considering possible changes in the Mathieu modulation frequency or QBO frequency.
Validation
A very simple validation to the meat of the analysis above exists. Two possible alternatives are to estimate (1) the ratio of the divergence (second derivative) over the signal, and (2) to machine learn on the divergence. The first alternative is brute force as shown in a prior post and is subject to noise; as whenever the signal goes through a zero crossing, the ratio gives rise to a singularity. This makes it difficult to get anything but a rough guide to the periodic perturbation. The second much more effective approach involves using an unbiased model discovery tool called Eureqa. What we do is take the SOI data, with as much of the long-term multi-scale oscillation removed (recall the earlier warning), and then apply the Eureqa search target as follows
D(sma(soi, 12), Time, 2) = f(sma(soi, 12), Time)
This is essentially asking Eureqa to map the second derivative of SOI (smoothed at the 12 month level by a sma filter) as a function of unknown parameters with respect to time and the SOI value. The search will evaluate trig functions such as sine and cosine and any other mathematical expressions that exhibit low complexity. The addition of the SOI term is a slight coercion or hint to Eureqa to guide it to find solutions of the Mathiue equation at the top of this post, i.e.
$$ frac{d^2x(t)}{dt^2} = f(x, t) $$
After Eureqa grinds away for several hours, it astoundingly discovers our Mathieu equation with very good error reduction. This is in spite of the ridiculously noisy second derivative computation. See Figure 8 below, and in particular look at the yellow highlighted solution located at a moderate complexity point along the Pareto frontier.
Fig 8 : The Eureqa machine fit to the second derivative of SOI finds the same Mathieu equation formulation postulated earlier.
The two higher-frequency sine waves have nearly the same time period as we identified in the QBO analysis ! And the Mathieu modulation of approximately 8 years exists as well ! This in effect provides a largely unbiased substantiation that the Mathieu model postulated earlier indeed exists in the data. One of Eureqa’s strong suits is that it’s machine fit has no conscience and will punish anyone with the truth if they dare to ask it to validate results.
Discussion
The mystery of what causes the observed periods of 433 days in the Chandler Wobble and the origin of the 8 year Mathieu modulation and whether that is connected to the Chandler wobble is a bit unsettling. I think that interactions exist in the 1 year seasonal cycle, the QBO cycle, the Mathieu equation characteristic frequency set by √a, the Mathieu equation modulation period set by ω, and the Chandler wobble period. I wanted to characterize these relationships more in this post, but will have to save it for later.
References
[1] Frandsen, Jannette B. “Sloshing Motions in Excited Tanks.” Journal of Computational Physics 196, no. 1 (2004): 53–87.
[2] Kawale, J., Liess, S., Kumar, V., Lall, U., & Ganguly, A. (2012, October). Mining time-lagged relationships in spatio-temporal climate data. In Intelligent Data Understanding (CIDU), 2012 Conference on (pp. 130-135). IEEE.
[4] http://judithcurry.com/2014/05/13/climate-dialogue-on-climate-sensitivity-and-transient-climate-response/#comment-552844
“Some people are under the illusion that ∆T = λ ∆F governs something. But it doesn’t. No need to be a skeptic to write the real and different governing equations – there are 5 of them. The first one is : div(V) = 0 where V is the velocity field.”
GeoEnergy Math 
This analysis has elicited a response from a well-regarded climatologist
A deeper analysis with a higher fidelity resolution in the QBO model should be able to pull out any slight seasonal details that the good professor is pointing out.
As far as not showing longer-range coherence in the time series, the analysis in this blog post contradicts that assertion.
LikeLike
WHT: Nice writeup. I’ve used Eureqa myself and it is a beautiful tool, especially for the math-challenged such as myself, but a lover of data and analysis. I particularly like its completely unbiased agnostic approach to picking out signals from multiple noisy series of data, which has very much helped me get comfortable with “the pause” and other nasty dissonance thrown out by the hordes of cherry-pickers out there. I was, however, frustrated by trying to do forecasting/predictions with seemingly predictable oscillations that never quite worked out. Good for me to have gained further insight into the need for much more complex physical models as tantalizing as the simpler highly correlated ones Eureqa so often cooks up.
LikeLike
Thanks,
I have found that Eureqa can be used to either explore data in a completely random fashion where you have no idea of the underlying model, or as a verification tool where you use it to cross-check what you have postulated by other means. The problem with the former is that it can lead you down some dead-ends if the correlation is spurious. That can happen with the latter in that it can reinforce a spurious model to begin with, but that provides a sanity check that you didn’t screw up on your math, even if it is invalid for the problem at hand.
LikeLike
This is an interesting recent paper using the same math but in a different discipline:
Shivamoggi, Bhimsen K. “Vortex motion in superfluid 4He: effects of normal fluid flow.” The European Physical Journal B 86.6 (2013): 1-7.
“… Mathieu’s equation describing the parametric amplification of Kelvin waves by the friction force. Equation (52) reveals again that the friction term associated with α influences the vortex motion in a qualitative way.”
The thinking is that the QBO is parametrically amplifying the ocean’s Kelvin waves in a similar manner, leading to ENSO. Remember that the physics and math do not change just because of the scale of the system.
LikeLike
Thanks for the tips about QBO and Kelvin waves. I wish they’d figure it out quick, I’m so tired of getting hammered by engineers carping about failed predictions as if the planet is something that fits on their lab bench and all you need is a good oscilloscope and a soldering iron to fix it.
Speaking of, I’ve been scratching around globally averaged model outputs for surface temps, CMIP5 and/or EMIC, anything else. I’m looking at fig 9.8 in AR5 ch 9, pg 768. I’d like some bleeding-edge numbers. Ideas?
LikeLike
The strongest correlation of the ENSO indices is with the Quasi-Biennial Oscillation (QBO) of upper atmospheric winds. The QBO is much more clearly periodic than ENSO, yet there seems to be some mutual interaction between the two.
As with many of these phenomenon, it is often difficult to determine which is the originating factor, or whether there is some type of self-organization going on between the factors.
So what I am struggling with in the ENSO / QBO correlation is how one behavior is so erratic (ENSO) while the other is much more predictable (QBO). The arrow of entropy would suggest that QBO would be the driver since it is much more ordered — yet what causes the QBO to create its more regular oscillation? The literature suggests that it is underlying oceanic waves, which points it back to ENSO.
So there is likely an underlying periodic nature to ENSO, which is only revealing its erratic nature via the sloshing dynamics of constructively or destructively interacting waves, simulated either by a Mathieu equation or a delayed action oscillator. Bottomline is that the order is there, but it is difficult to extract.
And once we can extract that order, prediction may become easier.
That’s the reason I started looking at ENSO. It is the biggest natural variability contributing factor to the global temperature — yet our inability to predict it becomes a tool for global warming skeptics to use as a bludgeon to claim uncertainty of outcomes. And therefore a way to marginalize climate scientists claims of continued warming.
LikeLike
A few years ago I conducted my first R “analysis” of ARGO data. It was simply fitting a sine wave to monthly temperature data at different depths for 45ºS. It was interesting to see a QBO-like pattern in the phase. That is, the phase pattern was long, short, long… relative to climatology.
I doubt if this of interest to you, but here it is anyway:
https://sites.google.com/site/climateadj/argo-sine-fitting
and an updated version:
https://sites.google.com/site/climateadj/argo-sine-fitting-2012
LikeLike
AJ, That is an interesting data reduction. The alternating oscillations could be a downwelling and upwelling effect tied to the same mechanism (of sloshing) that describes ENSO.
A sloshing of cool water into warm water will delay annual solar warming while a sloshing of warm water into already warm water will hasten the trend. This could lead to an erratic lag from year-to-year.
LikeLike
I used to think it was a variation in upwelling, but I’m now inclined to believe it is due to a small oscillation in the velocity of heat carrying currents. The slower the “river of heat” flows, the more time for diffusion and the cooler the river becomes (and vise-versa).
Here’s my explanation for the term “river of heat”. Heat uptake at the surface is greatest in the tropics and is pushed westward by the tradewinds until it runs into a continent. The heat is then pushed poleward and begins to slowly sink. Eventually the flow escapes the tropics, gets caught up in the eastward flow and continues to slowly sink. I think that this heat river is basically diffused somewhere after the 1000m mark and then the conversation shifts to intermediate/deep water formations.
If this is the case then a QBO-like signal in the surface tradewinds would likely also get transmitted to the equatorial current velocities and then to mid-latitudinal waters. The following link is relevant:
https://sites.google.com/site/climateadj/argo-animations
Figure 3 shows that the heat from the tropics is carried to the mid-latitudinal depths.
One can deduce from Figure 4 that the tropical heat slowly sinks and is carried poleward when it hits the western margin (i.e. a continent). Eventually it gets caught up in the westerlies (eastward flow) and continues to slowly sink.
So in my view, thanks to Dr. Bernoulli, an oscillation in stratospheric winds will correlate with a phase adjusted oscillation of mid-latitudinal ocean temperature at depth. This might imply cause and effect, but in a coupled system it’s hard to figure who’s zooming who 🙂
Then again, this could be shit. Best ask an Oceanographer.
LikeLike
What these engineers are trying to tell you is that if you really don’t understand something then you should stop pretending that you do. A good start is having an open mind and noticing that the ‘noise’ actually means something – it’s natural variation and it’s more important than too many people initially assumed. If people had left their prejudices outside the lab door and actually listened to critics (who in the main turned out to be correct) then climate science might have progressed somewhat quicker. As it stands it lurches from dogma to excuses to eventual acceptance that nature is not just standing still waiting for mankind to nudge it.
LikeLike
We use mathematical abstractions and simplified models to aid in understanding.
Not much more to it than that.
LikeLike
PW – “ENSO is a nonlinear property of climate, naturally varying but the onset of a phase is unpredictable.”
WHT – “As far as not showing longer-range coherence in the time series, the analysis in this blog post contradicts that assertion.”
It looks intriguing, but where is your ENSO prediction?
LikeLike
I am sorry Kevin, you are the one that is supposed to predict when I will be ready to make my own prediction.
So please make your prediction right now so that we can compare your skill in judging when I am good and ready.
LikeLike
WHT – why so snippy? Why would I try to guess when you’ll be ready? Just throwing out a random calendar date between now and never would be as insightful.
I am genuinely intrigued, but not sure why you dismiss PW’s assertion without having actually predicted ENSO correctly. If it’s wrong, it’s wrong – big deal. Just means more work needs to be done.
Of course Einstein spent 40 years trying to get rid of ‘spooky action at a distance’ – but quantum entanglement remains. Physics doesn’t care about our likes or dislikes. I hope you are onto something – but there are enough people with shitty attitudes around the web that I don’t need anymore of them.
LikeLike
What exactly is a prediction of future behavior of ENSO at this point in time going to prove?
If I make a prediction right now, it won’t prove anything other than I can take the equations and extrapolate some numbers. I don’t even know what it means when you say “I am genuinely intrigued, but not sure why you dismiss PW’s assertion without having actually predicted ENSO correctly. ” How can I have predicted ENSO correctly when the future has not even arrived yet?
It won’t prove anything because the future time span needed to prove this model out won’t be here for decades.
I wasn’t being snippy, just incredulous that people would think that a prediction would be useful, other than as a fallacious argument to raise the bar or move the goalposts.
Given that, I would also be very careful about making any predictions because you only get one shot, as people only remember the first prediction.
LikeLike
Thought you might like this site:
http://waterdata.usgs.gov/mn/nwis/uv?cb_00065=on&cb_00045=on&format=gif_mult_parms&site_no=05289800&period=&begin_date=2014-04-18&end_date=2014-06-05
The creek is sensitive to changes. It amplifies them.
LikeLike
Ragnaar,
Runoff is a good example of dispersion in action.
I have written on the topic in the past:
http://mobjectivist.blogspot.com/2010/09/hydrogeology-for-dummies.html
This is part of the dispersive math I developed when writing The Oil Conundrum book. The blog post has some broken charts but the book has everything in place.
Like Kevin said, runoff into Minnehaha Creek is a pretty trivial problem. That’s a minnie Ha Ha in terms of stumping the chump.
LikeLike
ragnaar – I think you are forgetting what happens to a creek – just about any creek – when rain falls; it collects not just the rain that falls directly on the creek, but also collects runoff from the surrounding area.
The creek is a watershed for an area much larger than that of the creek itself. As such it would be surprising if the creek did not rise by more than the amount of precipitation.
LikeLike
WHT, I really like this kind of systems analysis approach and you seem to handle the tools well.
“So what I am struggling with in the ENSO / QBO correlation is how one behavior is so erratic (ENSO) while the other is much more predictable (QBO). The arrow of entropy would suggest that QBO would be the driver since it is much more ordered — yet what causes the QBO to create its more regular oscillation? The literature suggests that it is underlying oceanic waves, which points it back to ENSO.”
The irregularity is probably just an interference pattern with other cycles.
You’re onto likely significance of the perigee cycle, how it is related in phase to lunar declination (latitude) will lead to variability. Tides are principally horizontal movement of water. When close perigees coincide with zero declination there will be a drawing of water away from high latitudes towards the tropics. Inversely as this happens near max declination.
Because of the push-pull nature of the tidal forces, declination either side of the equator has the same effect and the 18.6 period acts as a folded (fullwave rectified) cosine. This produces 9.3 in addition to 8.85. There is enough noise and averaging going on that these are often not resolvable. The mean frequency is 9.04. I suspect this is the 9.1 +/-… reported in several papers, eg recent BEST/Curry paper on European SAT found it in cross-correlation of AMO/PDO.
I found the same thing using SST directely , rather than the detended “indices”.
http://climategrog.wordpress.com/?attachment_id=754
Keeling & Whorf 1997 also pointed how 18.605 and 8.847 produce 6.00y. This interacting with 1y would give something rather close to Chandler. I don’t recall having see the latter pointed out anywhere.
That would seem to link your key elements together and provides for “oscillations” like ENSO etc appearing ‘erratic’.
I also found what may be evidence of QBO and 8.87 in spectral analysis of trade wind data and Dengler et al’s 4.5 years.
http://climategrog.wordpress.com/?attachment_id=283
The latter is speculative attribution but may provide clues to how all these systems interact and why everything looks erratic.
I like the hydrodynamic approach, water is the key to climate on Earth. Hopefully some of this will help provide the right drivers for your model.
LikeLike
Pingback: The QBOM | context/Earth
Pingback: Proxy confirmation of SOIM | context/Earth
Pingback: Sloshing Animation | context/Earth
Pingback: Correlation of time series | context/Earth
Pingback: The Southern Oscillation Index Model | context/Earth
Pingback: The Chandler Wobble and the SOIM | context/Earth
Pingback: The SOIM: substantiating the Chandler Wobble and tidal connection to ENSO | context/Earth
Pingback: Azimuth Project on El Ninos | context/Earth
Pingback: Demodulation and the SOIM | context/Earth
Pingback: Paper on Sloshing Model for ENSO | context/Earth
Pingback: Scaling El Nino | context/Earth
Pingback: Using Tidal Gauges to Estimate ENSO | context/Earth
Pingback: The QBOM | context/Earth
Pingback: Canonical Solution of Mathieu Equation for ENSO | context/Earth
Pingback: The Southern Oscillation Index Model | context/Earth