Here is an alternate derivation of reducing Laplace's tidal equations along 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,\\[2ex]{\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\qquad {\text{and}}\\[2ex]{\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 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 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,\\[2ex]{\frac {\partial u}{\partial t}}&+{\frac {1}{a}}{\frac {\partial }{\partial \lambda }}\left(g\zeta +U\right)=0\qquad {\text{and}}\\[2ex]{\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 diffEq.

So we can eventually get to a fit that looks like the actual QBO:

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.

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


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