Detection of limit cycle signatures of El Niño in models and observations using reservoir computing

The Jin–Timmermann (JT) model [38] was developed to study variations in the amplitude of ENSO variability. For convenience, the original version of the model is shown in appendix A. Here, we use the non-dimensional version described in [31] of which the governing equations are

$$\begin{align} \frac{\mathrm dx}{\mathrm dt} & = \rho\delta\left(x^2-ax\right)+x\left[x+y+c-c\text{tanh}{\left(x+z\right)}\right] -D_x\left(x,y,z\right)\eta_t ~, \end{align} \tag{ 1a }$$

$$\begin{align} \frac{\mathrm dy}{\mathrm dt} & = -\rho\delta\left(ay+x^2\right) +D_y\left(x,y,z\right)\eta_t ~, \end{align} \tag{ 1b }$$

$$\begin{align} \frac{\mathrm dz}{\mathrm dt} & = \delta\left(k-z-\frac{x}{2}\right) ~.\end{align} \tag{ 1c }$$

Here x is the dimensionless zonal temperature difference, y is the dimensionless temperature deviation between western Pacific and radiation equilibrium temperature and z is the dimensionless western thermocline depth anomaly. The physical meaning of the parameters ρ, k, c, a is explained in appendix A. The parameter δ is a measure of the strength of the ocean-atmospheric coupled processes and is taken as the main parameter here.

To include noise in the model, the approach described in [21] has been followed, where the noise components D_x and D_y are defined as:

$$\begin{align} D_x(x,y,z)& = \left[ (1+\rho\delta)x^2+xy+cx(1-\text{tanh}{(x+z)}\right]\sigma ~ , \end{align} \tag{ 2a }$$

$$\begin{align} D_y\left(x,y,z\right)& = \rho\delta x^2\sigma ~, \end{align} \tag{ 2b }$$

and η_t in (1) is a white-noise component with unit variance. For σ = 0, we obtain the deterministic version of the model. The equations (1) have been integrated using the Euler–Maruyama scheme, with a dimensional time step of one month (dimensionless time step is 0.29).

The ZC model is a so-called intermediate complexity ocean-atmosphere model and is one of the first models that was able to reasonably simulate ENSO [44]. The original version of the model was deterministic and the seasonal cycle was prescribed from observations. It consists of a set of seven two-dimensional partial differential equations for the surface temperature T, the zonal and meridional velocities in the surface ocean (u, v), those of the lower atmosphere (U, V), the ocean thermocline depth h and the atmospheric geopotential φ. In a subsequent variant, the mean state was also modeled and noise was introduced [9]; this is the model we use here. For convenience, the equations of the model are shown in appendix B. The main dimensionless parameter is the coupling strength µ. To numerically solve the model's equations the technique described in [40] has been applied.

The community earth system model (CESM) version 1 is a coupled climate model consisting of different model components, including ocean, sea ice, land surface and land ice models [17]. These different components exchange fluxes through a central coupler. The atmosphere is simulated by the community atmosphere model (CAM, [20]), while the oceans are simulated by the Parallel Ocean Program model (POP, [34]). We have used results from simulations with the CESM version 1.1.2 that is composed of 30 vertical layers in the atmosphere with a 0.9^∘ horizontal resolution, and 60 vertical layers in the ocean with a horizontal resolution of 1^∘. We use monthly mean data from a 100 years CESM simulation experiment under constant year 2000 forcing [41].

The observations analyzed are derived from the NOAA ERSST v5 [14], the PEODAS D20 datasets [43] and the ICOADS dataset [10]. From the NOAA ERSST v5 dataset, the NINO3.4 index has been extracted. From the PEODAS dataset the monthly mean thermocline anomalies have been extracted. Finally, from the ICOADS dataset, the monthly mean zonal wind-stress anomalies have been extracted. We focus on the period between January 1960 and September 2020.

The idea behind our approach is that the ML technique will be able to extract properties of the underlying dynamical system and hence can indicate whether a certain time series (characterized by parameters in a model, or by the present-day/historical climate in observations) contains limit-cycle characteristics. The RC framework is ideally suited for this, as it generates a high-dimensional dynamical system from input data, whose relevant patterns and dynamics are learned by the model during the training procedure.

Although the procedure to generate a RC has been well described elsewhere [22, 27], we briefly summarize the approach here, and introduce our notation. Given an input signal $u(n)\in R^{N_u}, n = 1,..,N_t$, where N_t is the total number of time steps, and a desired output signal $y^{\mathrm{target}}(n)\in R^{N_y}$, the model learns how to estimate an output signal $y(n)\in R^{N_y}$ as close as possible to $y^{\mathrm{target}}(n)$. In order to do that during the training procedure, an error measure $E(y,y^{target})$ is minimized. Different loss functions can be used and we have decided to apply the mean squared error (MSE) loss, i.e.

$$\begin{align} {\mathrm{MSE}}\left(y,y^{\mathrm{target}}\right) = \frac{1}{N_y}\sum_{i = 1}^{N_y}\left(\frac{1}{N_t}\sum_{n = 1}^{N_t}\left(y_i\left(n\right)-y^{\mathrm{target}}_i\left(n\right)\right)^2\right), \end{align} \tag{ 3 }$$

which is a common choice for regression problems.

Before the training procedure the input data u(n) are non-linearly expanded into a higher dimensional so-called reservoir space, generating a new signal $x(n)\in R^{N_x}$. This new representation of the data also contains temporal information and is based on the following update equations:

$$\begin{align} \tilde{x}\left(n\right)& = \text{tanh}\left( W^{in}\left[1;u\left(n\right)\right]+Wx\left(n-1\right)\right), \end{align} \tag{ 4a }$$

$$\begin{align} x\left(n\right)& = \left(1-\alpha\right)x\left(n-1\right)+\alpha \tilde{x}\left(n\right), \end{align} \tag{ 4b }$$

where $[1;u(n)]$ is short for stacking (here adding a top element 1 to the vector u(n)), which is done to introduce a bias term, and $\alpha \in (0,1]$ is the leaking rate. The two matrices $W_{in} \in R^{N_x\times(N_u+1)}$ and $W \in R^{N_x\times N_x}$ are generated randomly according to chosen hyperparameters. The non-zero elements of W and W_in are sampled from a symmetrical uniform distribution over the range $[-a,a]$. W is usually a sparse matrix, while W_in is dense, and the tuple $(W_{in},W,\alpha)$ defines the reservoir. The output layer is given by

$$\begin{align} y\left(n\right) = W^{out}\left[1;u\left(n\right);x\left(n\right)\right] , \end{align} \tag{ 5 }$$

and, during the training procedure, only the weights of the output matrix $W^{out}\in R^{N_y \times(N_u+N_x+1)}$ are estimated, minimizing $E(y, y^{target})$ through a linear regression procedure. We use a ridge regression to avoid overfitting, leading to the loss function ${\cal L}$:

$$\begin{align} {\cal L}\left(W^{out}\right) = E\left(y, y^{\mathrm{target}}\right)+\epsilon \sum_{i = 1}^{N_y} \sum_{j = 1}^{\left(N_u+N_x+1\right)}\left(W^{out}_{i,j}\right)^2 , \end{align} \tag{ 6 }$$

where ε is a regularization parameter. The hyperparameters are given by the dimension of the reservoir (N_x ), the sparsity of W's connections (as measured by an average degree of the network), the input scaling a, the leaking rate α, the regularization parameter ε and the spectral radius ρ_s of W. To determine the best set of hyperparameters in each problem, a Bayesian search has been applied based on criteria specified below. We will not list the values of the best hyperparameters in each of the results below, and the reader is referred to the codes and data which are publicly available (see Data Availability statement at the end of the paper).

After having determined optimal hyperparameters the RC can be transformed into an autonomous evolving system that can be used for prediction [18, 27]. Given an input sequence u(n), we train the RC to predict the target sequence $y(n) = u(n+1)$, i.e. perform single-step forecasts. After training, we can introduce feedback connections between the outputs at time step n and the inputs at the subsequent time step n + 1. This way, we generate a model that autonomously evolves in time without being driven by external inputs. This results in the following map:

$$\begin{align} u\left(n+1\right)& = y\left(n\right)\: = \:W^{out}\left[1;u\left(n\right);x\left(n\right)\right] , \end{align} \tag{ 7a }$$

$$\begin{align} x\left(n+1\right)& = \left(1-\alpha\right)x\left(n\right)+\alpha\text{tanh}{\left(Wx\left(n\right)+W^{in}\left[1;u\left(n+1\right)\right]\right)}.\end{align} \tag{ 7b }$$

Here, x(n) and $x(n+1)$ are the reservoir states at time step n and n + 1, while y(n) is the output at time step n and $u(n+1)$ is the input at the subsequent time n + 1.

Consider a simple system of two stochastic differential equations (SDEs):

$$\begin{align} \frac{\mathrm dx}{\mathrm dt}& = \mu x - \omega y - x^3 - xy^2 + \sigma \eta_{t} \end{align} \tag{ 8a }$$

$$\begin{align} \frac{\mathrm dy}{\mathrm dt}& = \omega x +\mu y - y^3 - x^2y + \sigma \eta_{t} \end{align} \tag{ 8b }$$

where µ and ω are parameters and η_t is a white-noise component with unit variance. The quantity σ is the noise amplitude and for σ = 0, a supercritical Hopf bifurcation occurs at µ = 0.

The values of µ we consider are a value far before the Hopf bifurcation ($\mu = -0.3$) and one (µ = 0.3) far after. Using the Euler–Maruyama method, we perform simulations of 2000 time steps $\Delta t = 0.1$ for σ = 0.08 and delete the first 1000 steps before training to guarantee that a statistical equilibrium has been reached. The best set of hyperparameters of the RC has been determined using a Bayesian search algorithm. For a particular set of hyperparameters, we evaluate the RC performances for both $\mu = -0.3$ and µ = 0.3 as follows. For each µ, the RC is trained using a noisy solution of (8). Next, the RC performances for a lead time of six time steps are evaluated using a different solution and using the root MSE (RMSE) as metric. Finally, the mean of the results obtained for the two µ values is computed. This approach allows us to estimate a single best hyperparameter set for both the µ values, increasing the generalization capabilities of the model and reducing the risk of overfitting. Table 1 shows that the prediction skills of the RC for both $\mu = -0.3$ and µ = 0.3 using the best set of hyperparameters, are very good. To account for the random nature of the RC, 100 runs have been performed (for each µ value separately), training and evaluating every time the RC with two different solutions of (8).

To better analyze the behavior of the self-evolving RC, we eliminate the bias term (see section 2.2), and define the output y(n) at time step n as $y(n) = W^{out}x(n)$. In this way, the self-evolving RC is given by:

$$\begin{align} u\left(n+1\right)& = y\left(n\right)\: = \:W^{out}x\left(n\right) , \end{align} \tag{ 9a }$$

$$\begin{align} x\left(n+1\right)& = \left(1-\alpha\right)x\left(n\right)+\alpha\text{tanh}{\left(Wx\left(n\right)+W^{in}u\left(n+1\right)\right)}, \end{align} \tag{ 9b }$$

where $W\in R^{N_x\times N_x}$, while $W^{in} \in R^{N_x \times N_u}$ and $W^{out} \in R^{N_y \times N_x}$ (the inputs and the outputs of the self-evolving RC have the same dimension so $N_u = N_y$). Using this formulation, we can compute the Jacobian matrix J of this discrete dynamical system (after the training procedure, W^out will be fixed so that the self-evolving RC will be deterministic) and compute the corresponding eigenvalues at the equilibrium point (where $x(n+1) = x(n)$) before and after the bifurcation. The number of eigenvalues depends on the reservoir dimension N_x , here $N_x = 8$.

Figures 1(a) and (b) show the output of the self-evolving RC for $\mu = -0.3$ (before the bifurcation) and µ = 0.3 (after the bifurcation), respectively. Figure 1(c) and (d) also show the eigenvalues of the self-evolving RC's Jacobian matrix of the mapping (9) at the equilibrium point. Before the Hopf bifurcation ($\mu = -0.3$), all the eigenvalues are inside the unit circle in the complex plane, and the self-evolving RC output is indeed a damped oscillation. After the bifurcation (µ = 0.3), the dominant eigenvalues have crossed the unit circle, and the self-evolving RC approaches a periodic solution.

Zoom In Zoom Out Reset image size

The first Hopf bifurcation in the deterministic JT model (σ = 0) for $[a,\rho,c,k]$ = [8,0.3224,2.3952,0.4032] occurs at $\delta = \delta_c = 0.21$. Hence periodic solutions are found slightly above $\delta \gt \delta_c$ and a damped oscillatory solution is found for $\delta \lt \delta_c$. If we add noise to the model with noise amplitude σ = 0.0095, the JT model shows an irregularly oscillatory solution for δ = 0.194 (figure 2(a)), while for δ = 0.214 bursting behavior is found (figure 2(b)).

Zoom In Zoom Out Reset image size

In the spectra for δ = 0.194, a peak at a frequency around 0.16 (corresponding to about 1.8 years) appears. For δ = 0.214, new frequency components appear but the frequency components related to the corresponding deterministic periodic solution are still clearly visible.

We now use the noisy solutions of the JT model to train RCs. For the two different values of δ (0.194 and 0.214), an RC was trained using 100 years of training data and all the three features $[x,y,z]$. Thereto we perform simulations of 300 years and delete the first 200 years before training to guarantee that a statistical equilibrium has been reached. The best set of hyperparameters for this scenario has been determined using a Bayesian search algorithm. For a particular set of hyperparameters, we evaluate the RCs performances for both δ = 0.194 and δ = 0.214 as follows. For each δ, a simulation is used to train the RC. Next, a new simulation (with a new noise realization) is used to evaluate the performance of the RC using a lead time of 6 months. Finally, the mean over the loss function for the two values of δ is computed. In this way, we determine a single best set of hyperparameters for both δ values. This was done again to avoid overfitting and increase the RC's generalization capabilities. The indicator we used to evaluate the RCs performances throughout the Bayesian search was the NRMSE (Normalized Root MSE), calculated by taking into account all features. The RMSE has been normalized separately for each variable, dividing it by the variance of the target values. Table 2 shows the prediction skills of the RC for both δ = 0.194 and δ = 0.214 (using the best set of hyperparameters). To account for the random nature of the RC, 100 runs have been performed (for each δ value separately), and the mean NRMSE values is shown (table 2) for each of the variables $[x,y,z]$ separately. The RCs prediction skills increase after the bifurcation (δ = 0.214), and NRMSE values are the smallest for the variable z (i.e. representing the thermocline depth).

The time series produced by the self-evolving RC (for best hyperparameters) are shown in figures 3(a) and (b). As expected from the results in section 3.1, the self-evolving RC displays behavior that corresponds well to that of the deterministic model, although it has stochastic features in the training procedure and in its random construction. Spectra of the time series of the self-evolving RC are shown in figures 3(c) and (d) using 150 different simulations. For each simulation, 1100 years of data has been generated and the spectrum has been computed only on the last 100 years. For δ = 0.214 (figure 3(d)), the RC has been able to identify the main frequency components of the corresponding deterministic JT solution (which is periodic for δ = 0.214). For δ = 0.194, the self-evolving RC most of the time converges to a constant value, but for a small number of runs, the self-evolving RC's spectrum is characterized by a small amplitude frequency component around 0.16 (figure 3(c)).

Zoom In Zoom Out Reset image size

Using the JT model, we can efficiently investigate the effects of the number and type of features on the RCs behavior. For each case, the same procedures described above have been used to both optimize the hyperparameters and evaluate the prediction skills of the RC. Using all the three variables $x, y$ and z as features, we first decreased the number of years used to train the model and found that the RC can distinguish subcritical and supercritical conditions down to a training data set of 50 years. When trained with 40 years of data instead of 100, the prediction skills of the RC decrease for both δ = 0.214 and δ = 0.194 (table 2). In this case, the self-evolving RC is not able anymore to reconstruct the underlying deterministic behavior of the JT model. For example, for δ = 0.214, the spectrum (not shown) displays variability on a wide range of frequencies. For δ = 0.194, the RC most of the time converges to a periodic solution, instead of a damped oscillation.

The prediction skills of the RCs when only the variable x is used to train it decrease for both δ = 0.214 and δ = 0.194 compared to those of the full features set (table 2). Also when we combine the variables $[x,y]$ the performances for both δ values decrease (table 2); the same holds when the RC is trained using the variables $[x,z]$ (table 2). When only an input feature x is used and 50 years of training data, the self-evolving RC produces very different results compared to using the full feature vector. The same result was found for a feature vector $[x,y]$, However, for a feature vector $[x,z]$, the self-evolving RC produces similar results as using the full feature vector. So both the zonal temperature difference (x) and the thermocline depth (z) are important for the RC to generate the correct behavior. This is in correspondence with the recharge-discharge oscillator view of ENSO, where both features are crucial [15, 38].

In the ZC model, we use the dimensionless coupling strength µ as the control parameter. For the standard values of the other parameters, a supercritical Hopf bifurcation occurs at µ = 3.0 in the deterministic model. Hence, for µ = 2.8, the Pacific climate state is stable (figure 4(a)), while for µ = 3.1, a periodic solution with a time scale of about 4 years appears (figure 4(b)). Under the wind-stress noise forcing with σ = 0.7 (see [9] for the dimensional values), the case µ = 2.8 shows substantial variability with a mean period of about 4 years (figure 4(c)). For µ = 3.1 the variability becomes more irregular due to the noise (figure 4(d)). These results are basically similar to those in [9], and are in correspondence with a stochastic Hopf bifurcation interpretation.

Zoom In Zoom Out Reset image size

The ZC-model results for all the combinations of five µ ($\mu = [2.5,2.6,2.7,2.8,3.1]$) and three σ ($\sigma = [0.5,0.7,1.0]$) values are used to train RCs, optimize the hyperparameters using a Bayesian search, assess the RCs performances and subsequently study the behavior of the self-evolving RC. This procedure is the same as used for the JT model described in the previous section. For the ZC model, we use a feature set consisting of the sea surface temperature in the western (T_W ) and eastern (T_E ) equatorial Pacific and the thermocline depth in the eastern (h_E ) and western (h_W ) equatorial Pacific, where the western and eastern regions are given by 120^∘ E–180^∘ E × 5^∘ S–5^∘ N and 150^∘ W–90^∘ W × 5.0^∘ S–5.0^∘ N, respectively.

We first consider the case where 60 years of training data are used. In this case, for each combination of µ and σ, data from 10 different stochastic ZC simulations of 100 years were generated. To avoid capturing the transient behavior, only the last 60 years of each simulation were considered. The best set of hyperparameters for this scenario has been determined again using a Bayesian search algorithm. In this case, we optimize the hyperparameters for each value of σ separately as follows. For a specific set of hyperparameters, the RC performances are evaluated for each σ, using all values of $\mu = [2.5, 2.6, 2.7, 2.8, 3.1]$. For each µ, the RC has been trained and evaluated using two different stochastic ZC simulations (sampled from the ten that were available). Finally, the mean of the different µ results was used for each σ to determine a best set of hyperparameters. This has been done again to avoid overfitting and generate a model able to generalize better. It was too difficult for the RC to generalize well also for all values of σ and for this reason different best hyperparameters sets have been determined. To evaluate the RC performances during the Bayesian search, a lead time of 6 months was considered, and the NRMSE was again used as the performance measure.

Table 3 shows the performances of the RCs (using the best sets of hyperparameters) for a lead time of 6 months. To account for the random nature of the reservoir 100 runs have been performed (for each combination of σ and µ separately) and the mean of the results have been computed. Similarly to the JT model results, the RCs show better prediction skills for supercritical values of µ, where the behavior of the ZC model is more regular and more robust with respect to the noise. Moreover, for the same value of µ, the performances remain similar when increasing the noise amplitude σ. One of the reasons behind this apparent resilience to the noise is connected to the strategies applied by the RC to avoid overfitting the noisy patterns in the training data. During training, the connections W_out , which are the only ones that can be influenced by the presence of the noise in the input signal (see section 2.2), are not estimated using a simple linear regression, but using a ridge regression, which limits the weights amplitude reducing the overfitting. Moreover, W_out does not directly map the inputs to the desired outputs as the input signal is first mapped into a non-linear high-dimensional space where it is easier to separate the noise components and identify the most relevant underlying patterns. The weights W_in and W (responsible for mapping the input data into the Reservoir space) are fixed and not influenced by the training data. This is another strategy that enhances the generalization capabilities of the RC, reducing the overfitting and the influence of noise.

As in the JT model results, the self-evolving RC solution converges either to a fixed point or a periodic solution. To show the results in a concise way, we define the criticality index C as follows:

$$\begin{align} C = \frac{\sqrt{\frac{1}{n_T}\sum_{n = 1}^{n_T} r^2\left(n\right)}}{\sqrt{\frac{1}{n_T}\sum_{n = 1}^{n_T} z^2\left(n\right)}} ~, \end{align} \tag{ 10 }$$

where n_T is the number of data points considered. Here, r(n) is the self-evolving RC solution for T_E and z(n) the same variable for the random ZC model realization, that has been used to train the RC. Before computing C, both r(n) and z(n) are centered around zero by subtracting the mean. Under subcritical conditions, we expect the RC solution to converge to a fixed point, so the nominator in (10) is much smaller than the denominator and C ≈ 0. For supercritical conditions, the RC solution is expected to converge to a periodic solution so C ≈ 1.

To compute the index C, we used 100 years simulations of the self-evolving RC and made sure that the amount of self-evolving RC data considered was the same as the amount of data used to train the RC. For example, for the 60 years case, the first 40 years of both the self-evolving RC and the training ZC simulation have been eliminated. The criticality index C is plotted as a function of µ for different σ values in figure 5(a). To produce this boxplot, for each combination of µ and σ, 100 simulations were performed. The index C is near zero for values of µ = 2.5 and µ = 2.6, for all σ values. in agreement with the subcritical conditions (see figure 5). For µ = 3.1, the value of C is near unity, in agreement with supercritical conditions. For the subcritical values µ = 2.7 and µ = 2.8 the value of C shows a large variance.

Zoom In Zoom Out Reset image size

The robustness of the results to the number of training years was investigated by using 40 and 20 years of training data, instead of 60 years. To estimate the best set of hyperparameters and evaluate the RC performances after the Bayesian search, the same methodology as the one applied for the 60 years case has been followed. The RCs, when trained with either 40 or 20 years of training data show prediction skills comparable to the ones obtained when it is trained using 60 years of data (not shown). The self-evolving reservoir results for 20 years of training data (as well as those for 40 years, not shown) are also fairly similar (figure 5(b)) to the ones obtained using 60 years of data, and hence the self-evolving RC produces adequate results also if we train it using a smaller amount of data.

For the deterministic ZC and JT models considered in the previous subsections, we precisely know the Hopf bifurcation and hence can use the criticality index C to determine whether the ENSO variability is sub- or supercritical. For the 100 years of de-seasonalized CESM data (model years 1200–1300) [41], we do not know whether a limit cycle is central for the ENSO variability. For the features, we use the NINO3.4 index (figure 6(a)), the anomaly of the depth of the 20 ^∘C isotherm over the area 120^∘ –280^∘ E × 5^∘ S–5^∘ N (figure 6(b)) and the zonal wind-stress anomaly over the area 120^∘ E–200^∘ E × 2.5^∘ S–2.5^∘ N (figure 6(c)). The depth of the 20 ^∘C isotherm for each grid point was determined using linear interpolation from the temperature field provided by the CESM. The surface zonal wind stress and NINO3.4 are direct outputs from the CESM model.

Zoom In Zoom Out Reset image size

Again, a Bayesian search has been applied to determine the best set of hyperparameters. For each set of hyperparameters, the RC performances for a lead time of 3, 6, and 9 months have been evaluated. To assess these RC performances, a cross-validation technique has been applied over the first 80 of the 100 years of the CESM data time series. The mean of the results obtained for each lead time has been taken as a measure of the RC performances using as metric the NRMSE computed in the same way described in section 3.2. After the validation phase, 5 different overlapping intervals of 40 years have been defined from the 80 years of validation data, using a sliding window of size 40 years and step 10 years. To evaluate the prediction skills of the RC, it has been trained using each of these periods to then assess its prediction performances on the last 20 years of the CESM time series, which have not been considered during the validation phase. This has been done to determine if the RC was able to perform well on completely unseen data. From these results, it was found that the RC has difficulties in predicting large amplitude ENSO events, especially for larger lead times (not shown).

For this reason, we decided to weight the contribution of different samples to the loss function. In particular, we define the NINO3.4 classes I: $\left[-\rho \sigma_w \mu_w,\rho \sigma_w \mu_w \right]$, II: $\left[\rho \sigma_w \mu_w, 2\rho \sigma_w \mu_w \right]$ and $\left[-2 \rho \sigma_w \mu_w, -\rho \sigma_w \mu_w \right]$, III: $\left[2 \rho \sigma_w \mu_w, 3 \rho \sigma_w \mu_w \right]$ and $\left[-3 \rho \sigma_w \mu_w, -2\rho \sigma_w \mu_w \right]$ and IV: $\left[3\rho \sigma_w \mu_w, \infty \right]$ and $\left[-\infty,-3\rho \sigma_w \mu_w \right]$. Here, µ_w and σ_w are the mean and standard deviation of the distribution of the training data and ρ is an additional hyperparameter optimized during the Bayesian search. With the number of samples belonging to each class being s_j , the weight w_j of each sample belonging to class j is determined from

$$\begin{align} w_j = \frac{\sum_{j = 1}^{n_c} s_j}{n_c s_j} \end{align} \tag{ 11 }$$

where here $n_c = 4$. Using this approach, the RC shows good prediction skills even for a lead time of 9 months, and is able to well capture the strong El Niño events, independently from the 40 years period used to train it (figure 7). To quantify the performance, we use the MAE (Mean Absolute Error) for the 9-month lead time predictions of the two largest El Niño (NINO3.4 $\unicode{x2A7E}2.5\,\,^\circ$C) and the 3 largest La Niña (NINO3.4 $\unicode{x2A7D}-2\,\,^\circ$C) events of the evaluation period. This has been done for all the training periods, and table 4 summarizes the results obtained. Applying the weights method, the RC is clearly able to predict strong El Niño events better. The RC still has difficulties in correctly predicting strong La Niña events. This is because our weighting scheme is based on the training data distribution, so the performance asymmetry reflects the asymmetry in the distribution of strong El Niño and La Niña events. The strongest Niño events have, in general, a larger amplitude than the strongest La Niña events (this is the case in CESM data as well as in real observations). In general, the number of strong amplitude El Niño events is larger compared to the number of strong amplitude La Niña events. The weighting method intensifies this asymmetry, leading to better performances in predicting positive events and to worse performances in predicting negative events. Table 4 also shows the RC prediction skills over the whole test interval (model years 1280–1300). Applying the weights method gives a small performance decrease. This is because, with the weights, we are giving more importance to strong events, so the predictions for the samples related to weaker temperatures will be less accurate.

Zoom In Zoom Out Reset image size

For these reasons, the weighting method was also applied to generate the self-evolving RC solution. The criticality index C has been computed taking into consideration the five overlapping periods of 40 years used to evaluate the RC prediction performances. Figure 8 shows the distribution of the results we have obtained in terms of the criticality index C for the different training periods. To produce this boxplot, for each training period, 300 different simulations were performed. A self-evolving RC solution of 100 years has been generated every time to then eliminate the first 60 years to avoid capturing the transient behavior. The RC has been trained every time using the best set of hyperparameters.

Zoom In Zoom Out Reset image size

For some of the intervals analyzed (figure 8(a)), the self-evolving RC produces a periodic solution, while for other intervals, its output converges most of the time to a constant value. More specific, for the intervals 1200–1240, 1210–1250 and 1230–1270, C is on average larger than 0.5. This means that the self-evolving RC converges most of the time to a periodic solution. The intervals 1220–1260 and 1240–1280 are characterized by an average value of $C \sim 0.1$.

To characterize the behavior of the self-evolving RC in more detail, we consider the change in amplitudes of the peaks of the NINO3.4 values to estimate the damping rate ζ of its solution r(n), with peak values p(n), as

$$\begin{align} \zeta = \frac{\delta}{\sqrt{4\pi^2+\delta^2}} ~ ; ~ \delta = \frac{1}{m}\text{log}{\frac{p\left(n\right)}{p\left(n+m\right)}} ~. \end{align} \tag{ 12 }$$

Here $p(n+m)$ is the amplitude of the peak of r after m cycles of the dominant period, with n fixed. When the time series is a sustained periodic signal, $\delta = \zeta = 0$. Under strong damping of the oscillation, δ will become larger than zero and, consequently, also ζ. Of course, the larger the value of m, the better the estimation of the damping ratio. Despite this, for small values of m (m = 5 was used from the starting solution with n = 0) it is already possible to clearly see the difference between a sustained and a damped oscillation.

In comparison with the index C, the damping ratio is less dependent on the self-evolving RCs output length. This advantage is significant, especially when the RC needs a lot of steps to converge. For the intervals 1200–1240 and 1230–1270 characterized by a large C in figure 8(a), the values of ζ are close to zero (figure 8(b)). This confirms that the self-evolving RC's output for these periods, on average, converges to a periodic solution. For the intervals 1220–1260 and 1240–1280, the damping rate is larger in correspondence with the lower value of C (figure 8(b)). Values of ζ for the interval 1210–1250 are larger than those of the intervals 1200–1240 and 1230–1270, but smaller than those of 1220–1260 and 1240–1280.

Without applying the weights method, the self-evolving RC produces damped oscillations for all the different intervals (not shown). This is due to the lower RC performances in predicting strong El Niño events for larger lead times (figure 7), so the RC seems less able to capture the relevant dynamics of the data properly. The logarithmic decrement technique to estimate the self-evolving solutions' damping ratio can be applied only to underdamped oscillations. When trained with the CESM data, the self-evolving RC always produces an underdamped oscillation or a sustained periodic solution (not shown), so the method can be applied. This is not the case when the RC is trained with the ZC data for small values of µ (e.g. µ = 2.5), where the self-evolving RC solutions are usually overdamped. For this reason, the behavior ζ for the self-evolving RC trained with the ZC data has not been analyzed in section 3.3.

The observational data set was linearly detrended, the seasonal cycle was removed, and then the same features as in the CESM were used, i.e. NINO3.4, the 20 ^∘C thermocline depth anomalies over the area 120^∘ E–280^∘ E × 5^∘ S–5^∘ N and the zonal surface wind-stress anomalies over the area 120^∘ E–200^∘ E × 2.5^∘ S–2.5^∘ N. A 3-month moving average was also applied to the observational time series (figure 9).

Zoom In Zoom Out Reset image size

The observations cover a period of approximately 60 years, from January 1960 to September 2020. The first 40 years have been used for validation and training, while the last 20 years have been used only for testing. The same hyperparameter search as described in the previous section has been applied.

In this case, we have determined a different set of best hyperparameters for each lead time by applying a separate Bayesian search. Figure 10 shows that the prediction skills of the RCs for a lead time of 3 and 6 months are relatively good, while they become low for a lead time of 9 months. However, due to the weighting method, the large amplitude El Niño events are better captured, especially for larger lead times. To quantify this performance, the MAE for the nine-month lead time RC's predictions of the three strongest El Niño (NINO3.4 $\unicode{x2A7E}1^{\,\circ}$C) and the strongest La Niña (NINO3.4 $\unicode{x2A7D}-1{\,^\circ}$C) events, over the test period January 2000- September 2020 has been computed (table 5). With the weights method, the RC can better predict strong El Niño events for a larger lead time. The model still has difficulties capturing the behavior of strong La Niña events, which are also poorly captured without applying the weights method. Table 5 also shows the RC prediction skills over the whole test period for different lead times. Applying the weights method, the RC performances over the whole test period decrease. As in the CESM case, we are giving larger importance to the samples related to strong events, so the performances in capturing the behavior of weaker events decrease.

Zoom In Zoom Out Reset image size

Because a separate validation procedure was performed for each lead time, also the criticality index C was computed for the best sets of hyperparameters identified for each lead time (figure 11). To estimate the index C and the damping ratio ζ, the same procedure described in the previous section has been followed, but only the period from 1960 to 2000 has been considered. Again to avoid capturing transient behavior, before computing the criticality index C, 100 years of self-evolving RC data have been generated, and only the last 40 years were retained. The self-evolving RC solutions for the real observations usually converge to a fixed point in less than 40 years (not shown). This is not a problem when we compute C, but for ζ this can be problematic. For this reason we have also considered the transient behavior to compute ζ, as otherwise it was not possible to adequately estimate the damping rate. The value of C (figure 11(a)) is close to zero for each set of hyperparameters, suggesting that the observed ENSO variability is subcritical.

Zoom In Zoom Out Reset image size

This is also confirmed by the damping ratio values (figure 11(b)) which are in the same order of magnitude as the ones obtained for the CESM intervals characterized by damped self-evolving RC solutions (figure 8(b)).