Asymmetries in the ENSO phase space

El Niño Southern Oscillation (ENSO) dynamics are best described by the recharge oscillator model, in which the eastern tropical Pacific sea surface temperatures (T) and subsurface heat content (thermocline depth; h) have an out-of-phase relationship. This defines a 2-dimensional phase space diagram between T and h. In an idealized, stochastically forced damped oscillator, the mean phase space diagram should be a perfectly symmetrical circle with a clockwise propagation over time. However, the observed phase space shows strong asymmetries. In this study we illustrate how the ENSO phase space can be used to discuss the phase-dependency of ENSO dynamics. A normalized spherical coordinate system allows the definition of phase-depending ENSO growth rates and phase transition speeds. Based on these we discuss the implications of the observed asymmetries with regards to the dynamics and predictability of ENSO; with a particular focus on the variations in the growth rate and coupling of ENSO along the oscillation cycle. Using linear and non-linear recharge oscillator models we will show how dynamics and noise are driving ENSO at different phases of the ENSO cycles. The results illustrate that the ENSO cycle with positive phase transitions is present in all phases but has strong variations in its strength. Much of these variations result from presenting the ENSO phase space with estimates of h based on the iso-thermal depth, that is not ideal as it is not out-of-phase with T. Future work should address how h can be estimated better, including aspects such as the vertical temperature gradients and the meridional or zonal range. We further illustrated that a non-linear growth rate of T can explain most of the observed non-linear phase space characteristics.


Introduction
The most widely used theoretical, conceptual model of the El Nino Southern Oscillation (ENSO) mode is the linear recharge oscillator (ReOsc) model [Burgers et al., 2005;Jin, 1997;Timmermann et al., 2018]. In this model ENSO is described by a cycle between subsurface upper ocean heat content and the sea surface temperature (SST), see sketch Fig. 1. Here increased (recharge) upper ocean heat content, which is measured by a deepening of the thermocline depth (h) leads to the development of El Nino SST anomalies in the eastern equatorial Pacific (T). These SST anomalies an analytical discussion of how non-linear dynamics affect the probabilities in the ENSO phase space. Takaheshi et al.
[2019] used the ENSO phase space to illustrate differences between a linear and non-linear model of ENSO.
Non-linear aspects of ENSO have been documented in the past in many different studies; including non-linearities in the amplitude, time evolution and patterns [Burgers and Stephenson 1999;Dewitte et al. 2013;Su et al. 2009;Ohba et al. 2010;Okumura and Deser 2010;Takahashi et al. 2011;Dommenget et al. 2013]. The ENSO phase space should be able to reflect the non-linearities in the amplitude and time evolution of ENSO and could potentially help to better understand the underlying dynamics of these two characteristics.
Several studies have tried to model ENSO non-linearities with the help of a non-linear variation of the ReOsc or other models [e.g., Choi et al. 2013;Levine et al. 2016;Frauen and Dommenget 2010]. They have been able to explain a number of different non-linear aspects of ENSO, but it is unclear how these approaches capture the asymmetries observed in the ENSO phase space.
Previous studies suggest that the predictability of ENSO is likely to be phase-depending [e.g., Dommenget et al. 2013;Timmermann et al. 2018]. Dommenget et al. [2013] found that strong La Nina events are likely to be more predictable than strong El Nino events of lead times of 7-11 months, due to the non-linear wind-SST relation. In contrast, Timmermann et al. [2018] argue that the transition from a recharge to an El Nino state is more predictable and that La Nina conditions are generally less predictable.
The aim of this study is to take a closer look at the ENSO phase space and present a detailed analysis of its observed characteristics. We aim to combine this analysis with a comparison of the observed phase space and the observed ReOsc model fits. By doing so, we would like to illustrate the extent a linear and non-linear ReOsc model can describe Fig. 1 Sketch of the ENSO recharge oscillator model dynamics. The ENSO cycle is clockwise with the heat content (h) in the vertical direction and sea surface temperature anomalies (T) in the horizontal direction. The three blue arrows in the horizontal planes mark wind anomalies resulting from T and explain the observed phase space characteristics. The ultimate aim of this study is to introduce the ENSO phase space characteristics as an effective way to present and analyse key ENSO dynamics.
The study is organised as follows: The following section introduces the data set used, the ReOsc model equations and the methods of estimating important parameters and statistics. Section 3 presents the results of the observed ENSO phase space, which is followed by a section on the linear ReOsc model and a section on a non-linear ReOsc model. In the final analysis section, we focus on the predictability of ENSO in the context of the ENSO phase space. Then the study is concluded with a summary and discussion.

Data, models and methods
Observed SST data is taken from the HADISST 1.1 data set for the period 1980 to 2019 [Rayner et al., 2003]. The monthly mean SST anomaly index region for T is the NINO3 region (150°W-90°W, 5°S-5°N). The thermocline depth anomalies, h, is estimated on the basis of the 20 o C thermocline depth (Z20) averaged over the equatorial Pacific (130°E-80°W, 5°S-5°N). Given the limitations in subsurface observations of temperature we use a combination of datasets to estimate monthly mean h: the 1980-2019 20 °C isotherm depths from the temperature analyses of the Bureau National Operations Centre (BNOC) at the Australian Bureau of Meteorology [Meinen and McPhaden 2000], the ocean reanalysis from SODA3 1980-2017[Carton and Giese 2008 and the CHOR AS and RL ocean reanalysis 1980-2010 [Yang et al., 2017]. The four datasets are combined to one long time series of T and h, thus repeating each year four times to better capture the variability. We also considered the GECCO2 reanalysis data [Köhl, 2015], but neglected it for this analysis, because it produced significantly different statistics compared to the other four datasets.
The ReOsc model is based on two tendency equations [Burgers et al., 2005]: with the growth rates of T(a 11 ) and h(a 22 ), the coupling parameters (a 12 and a 21 ) and the noise forcing terms (ζ T and ζ h ). The parameters of Eqs. [1][2] are estimated for the combined observations by multivariate linear regression the monthly mean tendencies of T and h against monthly mean T and h, respectively [Burgers et al. 2005;Jansen et al. 2009;Vijayeta and Dommenget 2018]. The residual of the linear regression fit can be interpreted as the random noise forcing, with the standard deviation (stdv) of the residuals being the stdv of the noise forcing for the T and h equations (ζ T and ζ h ). The values are shown in Table 1.
The ENSO phase space is presented by plotting T on the x-axis versus h on the y-axis, see Fig. 2. This Cartesian coordinate system can be transformed into a spherical coordinate system with the phase angle φ = 0 o in the h (y-direction) and 90 o in the T (x-direction). φ follows a clockwise rotation (Fig. 2). For this presentation it is useful to normalise T and h, by their respective standard deviation (Table 2) to get a non-dimensional presentation of the variables (T n and h n ). This normalization can also be applied to the fitted ReOsc model parameters ( Table 2).
In this normalized presentation we can define an ENSO system anomaly, S, as function of the two components T n and h n . The magnitude of S is constant for a constant radius and is not a function of the phase φ . Thus, the ENSO system is now described by the magnitude of S and φ . The tendencies of the ENSO system, as a function of the ENSO phase, are best described by the radial and tangential components. The radial component describes the tendency to move away from the origin (positive values) or towards it (negative values). The tangential component describes the tendency of the system to circle around the origin, with positive values indicating clockwise motion and negative values indicating anti-clockwise motion.
The analysis of observed or simulated statistics is based on monthly mean vales of T and h. The tendencies of T and h are estimated as: from the ReOsc model. This clockwise rotation is present for all phase angles, or more generally, in all four quarters of the diagram. Thus, positive heat content anomalies (h n ) lead to positive SST anomalies (T n ), which subsequently lead to negative h n , which lead to negative T n , and then back to a positive h n to complete the cycle. Therefore, the observed ENSO anomalies and their mean tendency do fit into the ReOsc model idea. However, there are some clear asymmetries present in the observed ENSO phase space diagram that are not expected from the idea of a linear ReOsc model. First of all, we can note that the ENSO system scatters much more towards positive T n values, than towards negative ones, and more towards negative h n values, than towards positive ones. Both asymmetries are expected from the well-known positive skewness in T n and negative skewness in h n [Trenberth 1997;Burgers and Stephenson 1999;Su et al. 2009].
It should be noted that much of the analysis could potentially also be done by analytically analysing Eqs. 1 and 2. However, this is not done here to provide a basis for applying this kind of analysis to any simulated or observed data. Figure 2 shows that the observed monthly mean ENSO phase space values are mostly a chaotic clustering around the origin, but for larger values the transition from one month to the next appears to be circling around the origin. This is indicating a transition in the phase space. To better illustrate how the system is developing in this phase space we compute the mean tendencies of the ENSO anomalies at different sections of the phase space, see vectors in Fig. 2.

Observed Phase Space
The mean tendencies of the ENSO anomalies highlight a clear clockwise rotation in the ENSO system, as expected between T n and h n with a lag zero correlation of zero, is not quite as it is observed.
The mean tendencies, as a function of the phase, are best described by the radial and tangential components, see Fig. 3b. For a stationary system, as ENSO is, the mean radial part over all phases must be zero, as the system is in average around the origin. The radial component is related to the growth rate of the system as it describes the tendency of the system to grow or decay.
The observed mean radial tendency is positive around 0 o and negative around 100 o and 220 o . This can also be noted by the mean tendency vectors in The smaller values indicate that the transition in the ENSO cycle is slowed down on average.
As mentioned above, the radial component of the tendencies are related to the growth rate of S. However, unlike in the ReOsc model (Eqs. [1-2]) where a 11 and a 22 are constant growth rates of T and h, which do not depend on T and h, the mean radial component as presented in Fig. 3b is a function of the mean S for each phase (e.g., the vectors in Fig. 2 depend on the mean S; they increase with distance to the origin).
Analogous to the ReOsc model growth rates, we can estimate a growth rate of S as a function of the phase by dividing the radial component ( Fig. 3b) by the mean S (red line in Fig. 2); see Fig. 5a. The structure of the growth rate of S is very similar to that of the radial component of the tendencies but can now be interpreted in the same way as the growth rate in the ReOsc model. We should note here that this statistical definition of the growth rate by definition is in average over all phases zero, and it represents the combined effect of the dynamics (T and h) and the noise forcing.
Similarly, the tangential component of the tendencies is also a function of the mean S for each phase. We can define a phase transition (angular speed) by dividing the tangential component ( Fig. 3b) by the mean S (red line in Fig. 2), see Fig. 5b. The phase transition is fastest between the El Nino state and the discharged state (~140 o ), and slowest between the discharged state and the La Nina state (~220 o ). The differences in the angular speeds are consistent with the differences in the likelihoods to be at different phases (Fig. 3a). ENSO phases which have large ENSO angular speeds are less likely to occur because ENSO transitions through these phases are relatively fast. In turn, phases which have small In the phase space we can note that we have larger scatter from about 30 o to 240 o and smaller scatter clockwise from about 240 o to 30 o . We can quantify these phase dependent probability distributions by estimating the probability statistics as a function of φ . The mean of S as function of φ is shown in Fig. 2 (red line). The mean is largest around 60 o to 90 o and smallest around 270 o to 360 o . Figure 3a shows the 2-dimensional probability density function. It shows the highest probabilities near the origin in quarter Q4 and larger probabilities for large S values in quarter Q1 to Q3, consistent with the scatter plot in Fig. 2. We can further estimate the probabilities of S values to be at different phase angles φ (black line in Fig. 3a). This shows higher probabilities to be in Q1 or Q3, and lower probabilities to be in Q2 or Q4. The probabilities are somewhat similar for Q2 and Q4 but show some enhanced likelihoods to be in Q1 compared to those in Q3. This shows that ENSO states between a recharge and an El Nino state (40 o ) have the highest probabilities. The lowest probabilities are the state at 90 o , and states before and after this. The probability distribution shifts a bit towards quarter Q2, if we only consider ENSO states with |S| > 1.0 (red line in Fig. 3a). It illustrates that large ENSO anomalies are primarily in phases from about 60 o to 240 o and less so from about 270 o to 360 o . Thus, large ENSO anomalies do not exist from the La Nina to the recharge state phase.
The scatter in Fig. 2 or the probability distribution between T n and h n in Fig. 3 shows an enhance likelihood along the diagonal from lower left (225 o ) to upper right (45 o ), which is reminiscent of a positive correlation between T n and h n . The observed correlation between T n and h n at a time lag of zero is 0.4, see Fig. 4a. Thus, the idealized concept of the ReOsc model which has an out-of-phase relation ReOsc model to illustrate what observed asymmetries are significant. We then discuss the linear ReOsc model with parameters fitted to the observed data to illustrate what kind of structures in the space phase can be explained by the observed linear dynamics. An idealised damped oscillator can be presented by the ReOsc model with all model parameters being symmetrical for T and h. To illustrate the characteristics of an idealised damped oscillator, we create a ReOsc model that is identical for both tendency equations. That is, the growth rates, coupling and strength of the noise forcing are the same magnitudes for both T and h. We chose the following parameters based on the normalised parameters in Table 2: We refer to this model as the idealised linear ReOsc model. The resulting phase space statistics of T n and h n are shown in Fig. 6a-c. Here we can note that all statistics are phase independent. The growth rate is zero at all phases, indicating that the mean tendencies in all phases only have a tangential part. Subsequently, the system is in statistical average moving around the origin in a perfect circle. In contrast, a ReOsc model without coupling (a 12 = a 21 = 0), which reduces to two unrelated red noise processes, does not have any mean ENSO angular speed are more likely to occur because the ENSO system spends more time in these phases.
The time to complete a full cycle (the mean period of ENSO) can be estimated by integrating the angular speed over all angles. This gives us a period for one cycle of about 42 months (3.5 yrs). This is consistent with the observed peak period in the T power spectrum (Fig. 4c). The phase transition time is however strongly variable within the cycle with the slowest transition of about 0.1 per months. This corresponds to a full cycle in about 5yrs. The fastest transition is about 0.3 per months which corresponds to a full cycle in 1.7yrs. These variations in the phase transitions are likely to contribute to the broadening of the power spectrum of ENSO.
The observed data record for the above analysis is only about 30yrs to 40yrs, which opens the question: 'to what extend are the observed characteristics statistically significant?' To address this question, we can use the linear ReOsc model; which will be discussed next.

Linear recharge Oscillator
The ReOsc model can help us understand what underlying dynamics are causing the asymmetries in the ENSO phase space. We start the discussion with an idealised linear idealized linear ReOsc model. Similarly, the distribution of asymmetries in the tangential part for quarter Q2 minus Q4 values (Fig. 7b) is well separated from the observed value, indicating that the observed variations in the tangential part of the tendencies is a clear signal.
We now focus on the linear ReOsc model with parameters as they result from a linear regression to the observed monthly mean T and h data, see Table 2. Using these parameters, we integrate the linear ReOsc model (Eqs. [1-2]) for 10 4 yrs and analyse the resulting normalized anomalies of monthly mean T n and h n . We refer to this model as the observed linear ReOsc model. Figure 8 shows statistics of the ENSO phase space for the observed linear ReOsc model. There are several interesting tendencies, and therefore, shows no mean propagation in the phase space cycle (Fig. 6d-f).
We can use the idealised linear ReOsc model to evaluate the statistical significance of the phase variations we noted for the observed ENSO phase space. For this, we integrate a 30yrs period with the idealized linear ReOsc model and repeat this 10 4 times to estimate distributions of important statistics for a 30 year observational period. In Fig. 7a, we show the distribution of the mean radial component of the tendencies around φ = 0 o ± 45 o from the 10 4 idealized linear ReOsc model integrations in comparison to the observed value. The observed value is clearly outside the modelled distribution, indicating that such large positive radial components of the tendencies cannot happen by chance in an ENSO anomalies, with much higher likelihoods around the Q2 quarter and lower in the Q4 quarter. This is not present in the observed linear ReOsc model (compare Fig. 3a with Fig. 8b). The observed growth rate around the phases of 0 o is much larger than around 180 o , which is not captured by the observed linear ReOsc model (compare Fig. 5a with Fig. 8c).
The phase-depending characteristics of the observed linear ReOsc model result from asymmetries in the ReOsc model parameters. The normalized model parameters (Table 2) allow us to compare the dynamics of the two tendency equations irrespective of the physical units of T and h. Here we can note that the main asymmetries in the two dynamical equations are in the growth rates. The growth rate of T is strongly negative and therefore T is damped. The growth rate of h is slightly positive and therefore unstable. In contrast, the coupling parameters and strength in noise forcing are nearly identical in magnitude for both equations. Thus, it is the asymmetry in the growth rates of T and h that cause the phase-depending characteristics of the observed linear ReOsc model.
The asymmetries in the growth rates have consequences for the growth and decay of the ENSO system at different phases. This is best illustrated if we split the total tendencies of the system (Eqs. [1-2]) into a dynamical part (first two terms on right hand side) and a noise driven part (last term on right hand side). The dynamical part can be calculated based on Eqs. [1-2] for any given T and h, and the noise aspects to point out in these statistics. First, we can note that in this observed linear ReOsc model, all statistics presented are phase dependent (Fig. 8), indicating that the observed linear ReOsc model does create structure in the phase space.
This contrasts with what may be expected from an idealised damped oscillator or the idealised linear ReOsc model (compare with Fig. 6). Secondly, we also note that all statistics are symmetric for opposing phases (e.g., shifts by 180 o ). This is a result of the linear approach in the ReOsc model, which assumes that the sign of T and h are irrelevant, and all feedbacks are symmetrical.
The probability distribution and tendencies in the ENSO phase space of the observed linear ReOsc model have some similarities with the observed statistics (compare Fig. 8a and b with Figs. 2 and 3a). This is also quantified by the correlations of these phase-depending statistics with those observed (see r-values in Fig. 8). The following similarities can be noted: (1) likelihoods are much higher in the Q1 and Q3 quarters relative to the Q2 and Q4 quarters. The phase transition speed is larger in Q2 and Q4 quarters relative to the Q1 and Q3 quarters.
In contrast, the observed linear ReOsc model also has clear deviations from the observed statistics. These are observed statistics that are asymmetric for opposing phases. The following important mismatches can be noted: There is a clear asymmetry in the observed probability of extreme at first seem strange since the noise part is by construction random and therefore should not have a preferred direction.
Here we need to remember that in the phase space diagram, we are considering conditional probabilities. For instance, if we are at S = 1 and φ = 30 o (Fig. 9a), then the ENSO system must have arrived at this point due to its past tendencies. Since the system is overall stationary and damped by the dynamics, it is by statistical average that it would arrive at this point, that is away from the origin, due to the noise. Thus, the noise is overall creating the variability leading to growth in general. This balance between dynamical damping and growth by the noise forcing is also part is estimated as the difference between the total tendencies and the dynamical part. Figure 9 shows the total tendencies and their dynamical and noise driven parts for the observed linear, idealized and uncoupled ReOsc models. Starting with the uncoupled ReOsc model (Fig. 9a) we can see that the mean tendencies are zero. For the radial part, which is related to the growth, we can see that the dynamical and noise parts of the tendency balance each other with the dynamical tendencies damping, therefore, pointing towards the origin. The noise part is pointing away from the origin, indicating that the noise is leading to the growth of the system. This may  Table 2. Therefore, we compute the dynamical tendencies for all phases using S = 1 and Eqs. [1-2] without the noise terms. Figure 10a-c shows the magnitudes, radial and tangential part of the dynamical tendencies. Since S = 1 for all phases, we can interpret the radial part as the dynamical growth rate and the tangential part as the dynamical phase transition (angular speed).
The dynamical growth rate of the system is directly related to the ReOsc model growth rate of T(a 11n ) and h(a 22n ). The strongly negative a 11n leads to strongly negative growth rate of the system when T n is large (at phases 90 o and 270 o ). The weakly positive a 22n leads to near zero growth rates when h n is small and weakly positive growth rates of the system when h n is large (at phases 0 o and 180 o ). The entire phase dependency of the growth rate of the system is directly related to these two extreme cases.
The phase dependencies of the magnitudes of the dynamical tendencies and the phase transition speed have similar structures as those of the dynamical growth rate of the system. However, they have maxima and minima at different phases. The phase dependency of the magnitude of the dynamical tendencies can best be understood if we look at the equation for the magnitude of the dynamical tendencies: Considering that a 21n ≈ −a 12n and |a 11n | |a 22n | we find: present in the observed linear and idealized ReOsc models ( Fig. 9b and c). The uncoupled ReOsc model has also no mean tangential tendencies for transition to another phase (Fig. 9a). Here both the dynamical and the noise part are zero. A mean phase transition in the ReOsc model is caused by the dynamical coupling between T and h [Lu et al. 2018], which is by construction zero in the uncoupled ReOsc model.
For the idealized ReOsc model, the mean dynamical and noise terms add up to have perfectly circular motion with the mean tendencies only having a tangential part and zero radial part. The dynamical part has a negative radial component; which is compensated by a positive radial noise as mentioned above, and a larger tangential component; which leads to the clockwise phase transition of the whole ENSO system. In this idealized ReOsc model all dynamical and noise parts of the tendencies are the same for all phases (Fig. 9b). Hence, it is entirely symmetrical in all parts.
The observed linear ReOsc model is similar to the idealized ReOsc model, but all elements of the mean tendencies are phase dependent, this includes the radial and tangential parts of both the dynamical and noise part (Fig. 9c). Starting with the dynamical part of the tendencies, we can see that the radial part (growth) is pointing towards the origin (is negative) at phases 90 o and 270 o , but is close to zero at phases 0 o and 180 o . Further, we can see that the overall tendencies and the tangential parts are larger at phases 315 o and 135 o , and smaller at phases 45 o and 215 o .
We can best understand these different phase dependencies of the dynamical tendencies by examining the ReOsc model Eqs. [1-2] using the normalized model parameters of T(a 11n ) and h(a 22n ) directly lead to a phase dependency of the dynamical magnitudes of the tendency; with maxima and minima at different phases than the growth rate. Similar computations (not shown) find that the dynamical transition speed has maxima and minima at phases similar to those of the magnitudes but shifted closer to the maxima and minima of the growth rate. ≈ (a 11n T n ) 2 + (a 12n h n ) 2 + (a 12n T n ) 2 − |2a 11n a 12n T n h n | (5) Here we can note that all terms add up if T and h have opposing signs (quarters Q2 and Q4), but if T and h have same signs, then the last term of the equation act against the other terms, reducing the magnitude of the tendencies. Consequently, the asymmetry in the dynamical growth rates In summary, we can say that the dynamical tendencies of the observed linear ReOsc model have phase dependencies resulting from the asymmetries in the dynamical growth rates of T and h. This makes the system anomalies decay when |T| is large and grow when |h| is large. The dynamical phase transition speed is largest at about 45 o before we reach the largest dynamical growth rates, and smallest at about 45 o after we reached the largest dynamical growth rates. This falls in-phase with the minima and maxima of the mean ENSO system anomalies (red line in Fig. 8a or 9c). As a result, the observed linear ReOsc model transitions fast when it is at phases with relatively small mean anomalies, and transitions slow at phases with relatively large mean anomalies.
Interestingly, the noise part of the tendencies of the observed linear ReOsc model is also phase dependent (Fig. 9c), although we have assumed by construction that the noise is purely random and not state dependent. Nevertheless, the phase-dependent dynamical parts of the tendency also lead to noise tendencies that are effectively phase-dependent. The radial part of the noise tendencies is always positive, but smaller at phases 30 o and 210 o , and larger at phases 120 o and 300 o . The tangential part of the noise tendencies is weak, but not zero. It is acting against the clockwise phase transition and is most strongly negative at phases 150 o and 330 o .

Non-linear recharge Oscillator
The above discussion has shown that the observed linear ReOsc model can capture a few characteristics of the observed ENSO phase space but has also illustrated that there are some asymmetries in the phase space that cannot be captured by a linear ReOsc model (e.g., asymmetries for opposing phases).
It is therefore instructive to consider non-linear ReOsc models to study how they would represent the ENSO phase space. Previous studies have suggested several different approaches to incorporate non-linear aspects of ENSO into the ReOsc model [e.g., Frauen and Dommenget 2010;Kim and An 2020;Levine et al. 2016]. These studies focused mostly on non-linear growth rates of T, state dependent noise or considered other non-linear elements in the ReOsc model. It is beyond the scope of this study to explore which non-linear process may explain the observed ENSO nonlinearities. However, we do want to provide an example to illustrate what a non-linear model could do and what it may be missing in the ENSO phase space.
We chose to focus on a non-linear growth rate of T and follow the approach of Frauen and Dommenget [Frauen and Dommenget 2010] by assuming a quadratic function. We non-linear parameters of this model suggest a stronger negative feedback for large negative T values, and a weaker or positive feedback for large positive T values. This is qualitatively similar to models suggested in previous studies [e.g., Frauen and Dommenget 2010;Geng et al. 2019;Kim and An 2020]. We integrate this model with the same noise forcing as for the linear model. We refer to this model as the non-linear ReOsc model.
Several phase space statistics of this model are shown Fig. 11. First, we can note that the non-linear ReOsc model has clear phase-depending statistics. Unlike the linear ReOsc model the statistics are also different for opposing phases. For instance, the phase-depending mean values (red line in Fig. 11a) are different at phases 90 o and 270 o (e.g., positive vs. negative T values).
The non-linear ReOsc model does capture the observed phase-depending characteristic that the observed linear therefore use Eq. [2] of the ReOsc model and change Eq. [1] to include a non-linear growth rate of T: We used a Nelder-Mead optimization scheme [Nelder and Mead 1965] to estimate the non-linear model parameters (a 11−2 , a 11−1 , a 11−0 ). The cost function for this optimization is based on integrating the model for 1000yrs and estimate the monthly mean distribution parameters. These are: the mean, stdv and skewness for T and h, and also the correlation between T and h. The root mean square of the differences in these statistics between the observed and the model values define the cost function of our optimization fit. The values of the non-linear model are shown in Table 1. The and 3a with Fig. 11a, b). The growth rate of the tendencies is larger for phases around 0 o than they are for phases around 180 o . The mean phase transition is slowest in Q3 and fastest in the Q2 quarter. While the non-linear ReOsc model is clearly closer to the observed phase space than the linear model there are also a number of significant mismatches between the non-linear ReOsc model captures, but also captures a few other characteristics. This is also quantified by larger correlation values in the phase-depending statistics (compare r values in Figs. 8 and 11). The following additional similarities to the observed can be noted: The mean and probabilities of the ENSO system are shifted away from the Q4 quarter and towards the other quarters Q1, Q2 and Q3 (compared Figs. 1 would correspond to a full ENSO cycle period of 35yrs to 1.7yrs, respectively. The extremely small phase transition values suggest that the ENSO cycle stalls and is potentially interrupted. This is also reflected in the total phase transition of the system (Fig. 11d). We can further note that small S anomalies transition faster relative to large S anomalies in quarters Q2 and Q3. The opposite is true in quarters Q1 and Q4.
These large variations in the phase transition does affect the power spectrum of T, by reducing the power at the peak oscillation period and increasing the power at all other frequencies (Fig. 4c). In particular, it enhances the decadal variations of T. Thus, the non-linearities in the growth rate of T is broadening the power spectrum, making it more realistic.

Predictability
We would expect that the variations in the ENSO characteristics at different phases of the ENSO cycle would lead to differences in the predictability of ENSO for different phases. We can get some approximation of how the observed ENSO may be predictable at different phases of the ENSO cycle by studying the predictability of the linear and non-linear ReOsc models discussed above.
First, we use the non-linear ReOsc model to start 8 ensembles of 100 members at different phases of the ENSO cycle with an initial S = 2, see Fig. 12. For each ensemble member a different realisation of the noise forcing was used, and the integration of the model was done for 12 months. We can define an ensemble mean S and phase φ for each forecast lead month, defining a mean position in the phase space. This is equivalent to a mean T and h (see solid lines in Fig. 12). The spread can be estimated by the distance of each ensemble member to the mean T and h for each forecast month and is shown as dashed lines in Fig. 12. Note, that in this representation we neglect the fact that the spread is not just in S, but also in the phase φ , as can be seen in the individual ensemble members in Fig. 12a.
The first example, starting at phase φ = 0 o (Fig. 12a), illustrates how the ensemble spreads out in terms of amplitude (S) and phase (φ). Some ensemble members decay in amplitude, while others grow or stay at the same amplitude. There are fairly large variations in the phase propagation, with some ensemble members propagating much further in the ENSO phase than the ensemble mean while others almost do not transition at all in the ENSO phase, but stay close to the initial phase.
The growth of the forecast ensembles is strongly depending on the initial starting phase (see Fig. 12b). Forecast ensembles that start at phases where the growth rate is model and the observations. The following deviations in the phase space can be noted: The probabilities of the nonlinear model are higher in Q3 than Q1, which is the opposite of what is observed. The observed growth rate asymmetry between 0 o and 180 o is much larger than in the non-linear model. In addition to these phase space deviations, we can also note that the cross-correlation between T and h of the non-linear model deviates quite significantly from the observed (Fig. 4a). In particular, when T leads h we find a strong underestimation of the cross-correlation in the nonlinear model.
The phase-dependency of the non-linear ReOsc model dynamical tendencies are overall like those of the linear ReOsc model, but are in detail more complex; see Fig. 10. First, we must note that in the non-linear model the dynamical tendencies do not just scale with S, as in the linear model, but will change its phase-dependency depending on the scale of S. This is illustrated by presenting the dynamical tendencies of the non-linear ReOsc model for three different values of S (S = 0.5, 1, 2; see Fig. 10). Here we can clearly note that the magnitude, radial (growth rate) and tangential part (phase transition) all vary more strongly as function of phases than in the linear model.
The variations in the radial part (growth rate) of the dynamical tendencies are mostly a function of T. This is an expected result given we have only made the growth rate of T non-linear and kept growth rate of h linear (see Eqs. [2 and 6]). The growth rate of T is less negative for large positive T values and much more negative for large negative T values, with a reverse relation for small T values (Fig. 10e). From this difference in the dynamical growth rate, we would have assumed a similar asymmetry in the overall growth rate of the non-linear ReOsc model. However, this is not observed (Fig. 11c). The growth rate of the non-linear ReOsc model is mostly symmetric in respect to T, but is asymmetric in respect to h. This suggests that the interaction with the phase-dependent tangential part of the tendencies and the noise forcing does lead to a significant shift in the asymmetries of the total system growth rate.
The variations in the tangential part (phase transition) of the dynamical tendencies in the non-linear ReOsc model are substantial (Fig. 10f). Before we discuss these large variations, we need to note that the coupling between T and h ( a 12 and a 21 ) in the non-linear ReOsc model are assumed to be still linear (e.g., not phase-depending). So, all non-linear and phase-dependent variations that we can note in the phase transition is a result of the non-linear growth rate of T.
The most extreme variations in the dynamical phase transition can be noted with minima in the Q1 and Q3, and maxima in Q2 and Q4 quarters. The extremes are ranging from 0.015 month − 1 for small values of S in Q1 to 0.315 month − 1 for small values of S in Q2. These extreme values the anomaly correlation between each forecast run and the control run for all forecast whose initial phase falls within ±15 o of the reference phase of T and h at different lead times, see Fig. 13.
Starting with the forecast of the idealised linear ReOsc model we can note a clear structure in the phase space for both T and h anomaly correlation skill ( Fig. 13a and  b). First, we have to recall that the idealised linear ReOsc model has no phase-depending ENSO characteristics, as discussed above, Subsequently, the structure that we see in the anomaly correlation skill scores is a characteristic of the phase space presentation; not a reflection of the characteristic of the idealised linear ReOsc model itself. For instance, at phases 0 o and 180 o , T is zero, and an anomaly correlation skill score for T at these phases must be zero too. Another instance is at a 3-month lead forecast starting at 3 40 o . This will on average end up at phase ~0 o and will therefore have small anomaly correlation skill scores for T. Accordingly, the minima and maxima will shift for different lead times, and the anomaly correlation skill scores for h will be shifted by 90 o .
The anomaly correlation skill scores for the observed linear ReOsc model are very similar to those of the idealized linear ReOsc model ( Fig. 13c and d). However, there are some small differences. Due to the asymmetries in ENSO amplitudes and phase transition speeds for this particular relatively large (see Fig. 11c) will have ensemble means that do not decay as fast, as seen for the ensembles starting at phases 0 o , 135 o , 180 o and 315 o . The opposite holds true for ensembles that start at phases where the growth rate is strongly negative (e.g., 90 o or 270 o ).
The phase transition speed is also strongly depending on the initial starting phase, with the fastest phase transition for the ensembles starting at phases 135 o and the slowest at 225 o (see Fig. 12b). Here it should be noted that all forecast ensembles have the same length in time (6mon), but appear to have different length in the phase space diagram due to their different phase transition speeds. The phase transition speed variations are strongly linked to the mean phase transition speeds (see Fig. 11d). The combination of the growth rate variations and phase transition variations splits the ENSO cycle into phases where the system clearly follows an ENSO cycle (around 315 o to 30 o and 135 o to 210 o ) and phases where the system is more or less collapsing and not propagating much (around 210 o to 300 o and 30 o to 90 o ).
We can evaluate the predictability of T and h in terms of the anomaly correlation skill as a function of phase within the ENSO cycle; based on the linear and non-linear ReOsc models. For this, we integrate a long control simulation, from which we start one additional simulation with different noise forcing every 60 months for a 9-month lead forecast. We do this 3.6 • 10 4 times, which roughly gives us about 100 forecasts for every 1 o of the ENSO cycle. We then estimate than the recharge state (at around 0 o − 30 o ). It is remarkable that we observe stronger non-linearities in the forecast skill of h than in T, considering that the non-linear ReOsc model discussed here is only non-linear in the tendencies of T, but is linear in the tendencies of h.

Summary and discussion
In this study we introduced the ENSO phase space for a detailed analysis of the ENSO dynamics. The observed model, the correlation skill score for h is not 90 o but only about 60 o out-of-phase with those of T.
The non-linear ReOsc model shows some clear asymmetries in anomaly correlation skill scores that are different from those of the linear ReOsc models. Anomaly correlation skill scores for T are in general larger in quarters Q1 and Q2 and lower in Q3 and Q4. Asymmetries are even more pronounced for correlation skill scores in h; with much larger skill scores in Q2 (with values at around 0.6) than in the Q4 quarter (values around 0.2). This suggests that the discharge state of ENSO (at around 180 o ) is much more predictable The positive in-phase (lag zero) correlation between T and h is not ideal in the context of the ReOsc model, suggesting that this is not an accurate presentation of the ENSO phase space as it assumes that T and h should be out-ofphase (zero correlation at lag zero). Other studies assume that the western equatorial thermocline is a good presentation of the ReOsc model [e.g., Jin 1997b; Chen et al. 2021], but western equatorial thermocline has a significant negative correlation with T at lag zero [e.g., Chen et al. 2021].
It is likely that the Z20 estimate of the thermocline depth (h) is causing a problem. Vijayeta [2020] analysed how differences in the estimation of h affects the ReOsc model presentation. The study found that a more accurate estimation of h, by a maximum temperature gradient approach, finds a nearly perfect out-of-phase correlation between T and h. This suggests that a better estimate of h would improve the ENSO phase space presentation.
The ReOsc model with a non-linear growth rate for T can explain most of the asymmetries in the observed phase space that can otherwise not be explained by the linear ReOsc model. A non-linear growth rate for T reproduces the observed shift in the likelihoods for large ENSO anomalies away from quarter Q4 and towards the quarter Q1-Q3. It further reproduces the strongly reduced phase transition speed in quarter Q3 (discharge to La Nina state) and the enhanced phase transition speed in quarter Q2 (El Nino to discharge state).
The variations in ENSO phase transition speed, as captured by the non-linear ReOsc model, lead to a more realistic power spectrum, with a broader interannual peak and enhanced decadal variability. The latter is consistent with earlier studies suggesting that ENSO non-linearity causes decadal ENSO variability [Rodgers et al. 2004;Wittenberg et al. 2014;. Here it is important to note that the mean ENSO period is primarily controlled by the coupling parameters [Lu et al. 2018], which have been kept linear in this model.
However, the ReOsc model with a non-linear growth rate for T cannot explain all aspects of the observed ENSO phase space. In particular, the observed lag-lead cross-correlation between T and h, with enhanced cross-correlation when T leads h are not well captured by the model.
The phase-depending ENSO characteristics should affect the predictability of ENSO. The non-linear ReOsc model suggests that ENSO predictability changes along the phases depending on the lead-time of the prediction. It affects the amplitude and phase transition differently, whilst also being different for T and h, respectively. In particular, h is most predictable in quarter Q2.
The ENSO phase space presentation introduced here provides many opportunities for further studies. A key aspect that needs to be addressed in future studies is the in-phase ENSO phase space showed several interesting asymmetries that reflect important aspects of ENSO dynamics. In agreement with Kessler [2002], we find that the probability distribution of ENSO phases has some clear asymmetries for large ENSO amplitudes, with lower probabilities to be within the Q4 quarter (La Nina to recharge state).
An important aspect of the ENSO phase diagram is that it allows the analyses of ENSO tendencies as a function of the ENSO phase. The spherical coordinate system of the ENSO phase space diagram allows us to define tendencies in the radial and tangential direction. A normalization of the radial tendencies defines the ENSO system growth rate as a function of the phase. While by construction, the mean growth rate in this definition must be zero, the growth rate at different phases shows clear deviations from zero, with positive growth rates around and after the recharge state (330 o to 45 o ) and negative growth rates around the El Nino (70 o to 120 o ) and La Nina states (210 o to 270).
A normalization of the tangential tendencies defines an ENSO system phase transition speed, which, if integrated over the whole cycle, gives an estimate of the ENSO period. The mean observed phase transition speed varies substantially as a function of the ENSO phase; with fast transitions in quarter Q2 (after the El Nino state) and slowest around the La Nina state (220 o to 260). This is somewhat consistent with the argument put forward in Kessler [2002] that ENSO is more event-like rather than a cycle, where it remains in a weak La Nina-like states for longer periods of time. However, the phase transition speed is significantly positive for all phase of the ENSO cycle, also supporting the idea that ENSO is indeed cyclic, though the speed or "clearness" of the phase transitions vary substantially in the ENSO cycle.
The underlying dynamical cause for the observed structures in the ENSO phase space is best analysed with a linear or non-linear ReOsc model. We illustrated that a linear model can explain some of the observed structures in the ENSO phase space and a non-linear ReOsc model can explain most of the remaining asymmetries.
A fit of a linear ReOsc model to the observed data and a normalization of the units reveal an asymmetry in the growth rate of T and h; with a negative growth rate for T and a weakly positive growth rate of h. The coupling parameters and strength of the noise forcing show no asymmetries. The asymmetry in the growth rates reflects a positive in-phase (lag zero) correlation between T and h. This positive correlation explains the observed characteristics in the ENSO phase space that are symmetric for opposing phases (shift by 180 o ); including enhanced growth rates at phases around 0 o and 180 o , and reduced growth rates at phases around 90 o and 270 o . This explains the enhanced phase transition speeds in the quarters Q2 and Q4 and reduced phase transition speeds in the quarters Q1 and Q3. correlation between T and h, which dominates the ENSO phase space characteristics and therefore potentially hides more interesting aspects to phase-depending ENSO dynamics. This is most likely related to how h is estimated by Z20 rather than true vertical profile gradient methods. More generally, other aspects of estimating of h, such as the meridional or zonal range, may affect the ENSO phase space representation.
A further aspect that has not been discussed here is the seasonal changes in ENSO dynamics [e.g., Li 1997;Tziperman et al. 1998;McGregor et al. 2013;Zhu et al. 2015;Dommenget and Yu 2016]. It needs to be considered that each quarter of the ENSO phase space should be transiting through all four seasons of the year. We therefore would expect seasonal variations in the ENSO phase space in all four quarters.
The discussion presented here for the dynamical phase of ENSO can also be applied for other climate modes. The Madden-Julian oscillation (MJO), for instance, has a welldefined dynamical phase space [e.g., Wheeler and Hendon 2004;Oliver and Thompson 2016]. The discussion, presented here for ENSO, could be applied for the MJO or other climate modes in a similar way.
Acknowledgements This study was motived by the Honours Bachelor of Science project of Maryam Al Ansari in 2021. We like to thank Tobias Bayr, Shayne McGregor, Peter van Rensch and the two anonymous reviewers for helpful comments and discussions. This study was supported by the Australian Research Council (ARC), discovery project "Improving projections of regional sea level rise and their credibility" (DP200102329) and Centre of Excellence for Climate Extremes (Grant Number: CE170100023).

Funding Open Access funding enabled and organized by CAUL and its Member Institutions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.