One thing to remember about 2nd-order DiffEq formulations is that the results are attenuated based on the frequency of the RHS forcing. So that if a natural resonant frequency of $\omega_0$ =1.5 rads/year is excited with equal magnitudes of frequencies $\omega$ = 1, $\pi$, $2\pi$, and $4\pi$ rads\year, the amount of attenuation will increase in the order, according to this fraction :

$ \frac{1}{1+(\omega/\omega_0)^2} $

And a Mathieu-type non-linearity in the DiffEq can further selectively amplify or attenuate frequencies.

This helps explain how lower frequencies of relative small magnitudes can emerge, even though they aren't the strongest perturbation. And why the higher frequencies don't become erratic as well, since they are not as close to the threshold of resonance as the lower frequencies are.

The 2-year delay differential is a convenient way of notching each of the $\pi$, $2\pi$, and $4\pi$ forcing factors from contributing to the lower frequency ENSO dynamics. They may still be there but whenever anyone plots SOI, typically a yearly smoothing filter is applied beforehand so they don't appear in the chart. Otherwise , the SOI will get obscured by this seasonal noise.

I can stand to be corrected on this interpretation. The main objective I have is to understand the rationale for being able to isolate the ENSO features in the tidal gauge data. I am concerned that the $\pi$, $2\pi$, and $4\pi$ forcing factors may all have an impact and that by throwing this data away, and in particular, when we potentially have a non-linear dynamic occurring, may not be wise.