Application of the Conditional Nonlinear Local Lyapunov Exponent to Second-Kind Predictability

The atmosphere is a chaotic system (Lorenz, 1963), which inherently possesses a limit to predictability due to its solutions’ sensitivity to initial conditions. According to Lorenz (1975), predictability studies mainly fall into two categories based on error sources. Those studies of first-kind predictability analyze the effects of the atmospheric system’s initial uncertainties on future state predictions, while studies of second-kind predictability examine the system's predictability response to boundary conditions (Lorenz, 1975; Schneider and Griffies, 1999). Since the predictable information derived from initial conditions tends to dissipate quickly, first-kind predictability limit of daily weather forecast is typically around two weeks. However, long-term predictability of the atmosphere can be achieved through interactions with external forcings, such as the ocean, land, and sea ice, which exhibit low-frequency variability (Shukla, 1981; Meehl et al., 2009). For example, numerous studies illustrate that El Niño-Southern Oscillation (ENSO) is a major global source of seasonal to interannual predictability (Tang et al., 2018).

Various methods have been proposed to estimate second-kind predictability (Lorenz, 1969; Kumar, 2009; Ren et al., 2009; Meehl et al., 2014; Zhang et al., 2022). Among them, the signal-to-noise ratio (SNR), which is defined as the fraction of variance of the “climate signal” due to the external forcings to that of the “climate noise” due to intrinsic random variations, has been widely used to quantify the degree of predictability caused by external forcings (Jones, 1975; Madden, 1976; Kumar, 2009; Strommen and Palmer, 2019). However, this method has limitations. It fails to account for the nonlinear interaction between the climate signal and noise, which may lead to underestimation of predictability (Scaife and Smith, 2018). Also, it only provides a qualitative (not quantitative) measure of the predictability due to external forcings, such as the information on regions with either higher or lower potential predictability (Zhang et al., 2022).

Numerous studies have attempted to quantify the influence of external forcings directly on the predictability limit through experiments on numerical models (Reichler and Roads, 2003; Guo et al., 2012; Zhao et al., 2021; Yong et al., 2023). However, these models are imperfect representations of the real world due to the errors associated with numerical and physical approximations, which can lead to uncertainty in the results (Bauer et al., 2015; Zhang and Kirtman, 2019; Duan et al., 2022). Given the vast amounts of atmospheric and oceanic observational datasets that contain almost all aspects of their movement and evolution (Guo et al., 2015; Runge et al., 2019), it becomes feasible to investigate the limit of predictability using these datasets.

Based on the theoretical framework of the Lyapunov exponent (Oseledec, 1968; Krishnamurthy, 2019), Ding and Li (2007) proposed a new method known as the nonlinear local Lyapunov exponent (NLLE) to quantify the predictability limit of the first kind. Along with the proposed new method, an algorithm called local dynamical analogs (LDAs) was devised to estimate the NLLE based on observational or experimental time series (Li and Ding, 2011). The LDAs algorithm was designed to search for LDAs from time series, thereby eliminating the possibility of generating a spurious estimate of the NLLE. Based on the new algorithm, the NLLE has been applied to quantify the predictability limit arising from initial weather and climate conditions across various timescales in recent years (Zhong et al., 2018a, b, 2021; Hou et al., 2022; Mengist and Seo, 2022; Duan et al., 2023).

To separately investigate the effect of external forcings on predictability, Li and Ding (2015) extended the NLLE method and proposed the conditional nonlinear local Lyapunov exponent (CNLLE). As an extension of the NLLE, this method has shown advantages in quantifying the impact of external forcings on the limit of predictability when applied in simple models. In those cases where the system’s equations of motion were explicitly known, such as the Lorenz-63 model coupled with different external forcings (Lorenz, 1963; Peña and Kalnay, 2004), we can directly obtain the mean CNLLE by numerical integration of these coupled models (Ding and Li, 2012; Zhao et al., 2021). However, a practical problem is that existing numerical models inevitably have model errors because the equations governing the evolution of the atmosphere and the influences of external forcings are either incomplete or unknown. Therefore, the estimation of the CNLLE from observational datasets can serve as a supplementary method for existing predictability studies. Unfortunately, an algorithm for estimating the mean CNLLE from observational time series has not been established.

The present study aims to fill this gap, by developing an algorithm to estimate the CNLLE based on the observational data. As an example, the CNLLE approach is applied to quantify the contribution of ENSO to the predictability limit of global oceanic and atmospheric fields [e.g., sea surface temperature (SST), sea level pressure (SLP), and geopotential height (GHT)]. As mentioned above, ENSO, mainly driven by abnormal SST in the central-eastern Pacific Ocean, functions as the primary predictable signal for the atmosphere. Knowledge of its contribution to the global atmospheric and oceanic prediction horizon can provide useful guidance to improve the forecasting skills for weather and climate.

According to the NLLE, the predictability limit can be quantified by the rate of error growth. However, this predictability may vary under different external conditions. Consequently, the CNLLE aims to further investigate the error growth rate influenced by external forcings. Let x and y represent the internal and external systems, respectively. The evolution of initial perturbations within the internal system, constrained by external forcings, is given by

where $ {{\boldsymbol{\delta}} _{{\mathrm{x}}|{\mathrm{y}}}}\left( t \right) = {\left[ {{\delta _1}\left( t \right),{\delta _2}\left( t \right), \cdots ,{\delta _n}\left( t \right)} \right]^{\mathrm{T}}} $ represents a perturbation vector of x constrained by external forcings, in which the superscript ${}^{\mathrm{T}} $ represents the transpose and the subscript x|y indicates it is constrained by external forcings (the same below). $ {\boldsymbol{\eta}} \left[ {{{\bf{x}}_{{\mathrm{x}}|{\mathrm{y}}}}\left( {{t_0}} \right),{{\boldsymbol{\delta}} _{{\mathrm{x|y}}}}\left( {{t_0}} \right),\tau } \right] $ is the nonlinear error propagator, which is responsible for propagating the initial perturbation $ {{\boldsymbol{\delta}} _{{\mathrm{x}}|{\mathrm{y}}}}\left( {{t_0}} \right) $ at the initial time ${t_0}$ toward ${{\boldsymbol{\delta}} _{{\mathrm{x}}|{\mathrm{y}}}}\left( {t + \tau } \right)$ after the evolution time $\tau $ (Ding and Li, 2007). $ {{\bf{x}}_{{\mathrm{x}}|{\mathrm{y}}}}\left( {{t_0}} \right) $ represents the initial value vector of x when under the influence of external forcings. The CNLLE ${\lambda _{{\text{x|y}}}}$ is defined as:

where $ {\lambda _{{\text{x|y}}}} $ depends on the initial state $ {{\bf{x}}_{{\text{x|y}}}}\left( {{t_0}} \right) $, the initial perturbation ${{\boldsymbol{\delta}} _{{\text{x|y}}}}\left( {{t_0}} \right)$ and the integration time $ \tau $, and $\left\| {\;} \right\|$ indicates the value of the norm for a vector. Equation (2) represents the average logarithmic growth of initial error from ${t_0}$ to ${t_0} + \tau $ driven by external forcings. The mean CNLLE over the global attractor of a dynamical system is given by:

where ${{\Omega }}$ represents the domain of the global attractor of the system and ${\left\langle {\;} \right\rangle _N}$ indicates an ensemble average of sufficiently large size samples. The mean relative growth of the initial error (RGIE) can be obtained by:

The RGIE is a function of time. Experiments in a nonlinear model coupled with external forcings show that the RGIE will first grow rapidly and then increase slowly due to nonlinearity before reaching saturation (Ding and Li, 2012; Sun and Zhang, 2020). The saturation point can be considered as the moment when all information from predictable signals is lost and prediction becomes meaningless, which is the limit of predictability (Li and Ding, 2011).

As stated above, the CNLLE can be determined if the governing equations of the dynamical system are explicitly known. However, knowledge of the equations for atmospheric motion and its interaction with external forcings remains limited, making it difficult to calculate the CNLLE theoretically. Therefore, we developed an algorithm based on LDAs to estimate the CNLLE from experimental or observational data (Fig. 1). The idea of the CNLLE method is to examine the growth rate of the distance (or error) between two states where the trajectories of corresponding external forcings are close in phase space. Let $ {\bf{Y}} = \left( {{\bf{y}}\left( {{t_i}} \right),i = 1,2, \cdots ,N} \right) $ be the observed data of the external forcings, represented by thick arrows in Fig. 1, with a time length of N. ${\bf{y}}\left( {{t_i}} \right) = {\left[ {{y_1}\left( {{t_i}} \right),{y_2}\left( {{t_i}} \right), \cdots ,{y_L}\left( {{t_i}} \right)} \right]^{\mathrm{T}}}$ is the state of external forcings at time ${t_i}$, where $L$ represents the number of observed external factors. In the schematic (Fig. 1), the state of external forcings is represented by the variation of the thick arrow around the black dots. Let ${\bf{X}} = \left( {{\bf{x}}\left( {{t_i}} \right),i = 1,2, \cdots ,N} \right)$ be the observed data of affected variables (in contrast to the external forcings), represented by thin arrows in Fig. 1. ${\bf{x}}\left( {{t_i}} \right) = {[{x_1}\left( {{t_i}} \right),{x_2}\left( {{t_i}} \right), \cdots ,{x_M}\left( {{t_i}} \right)]^{\mathrm{T}}}$ denotes the state of the affected variable at time ${t_i}$, where M represents the number of observed affected variables. Similarly, the state of affected variables is represented by the variation of the thin arrow around the dots in the schematic (Fig. 1).

Figure 1.  Schematic of the evolutionary procedure used to estimate the CNLLE from observed (experimental) data. The black arrow represents the reference trajectories, while the lighter gray arrows correspond to the analogous trajectories. The external forcing is indicated by thick arrows, while the affected variable is illustrated by thin arrows.

Moreover, we call the observed time series the reference trajectories, illustrated by the black arrow in Fig. 1. In this representation, the external forcing is denoted by a thick arrow, while the affected variable is depicted by a thin arrow. As expounded earlier, the purpose of the CNLLE is to examine the growth rate of the distance between two states—specifically when the trajectories of corresponding external forcings are close in phase space. Consequently, it is necessary to find trajectories that are closely aligned with each reference state of the external forcings; these are referred to as analogous trajectories, illustrated by the lighter gray arrows in Fig. 1. Therefore, for each state of the reference trajectories its paired analogous trajectories share analogous external forcings initially, but with the progression of time, the initially close states will gradually spread in different directions.

There are four steps to calculate the CNLLE from experimental or observational data, as described below:

As previously mentioned, the fundamental concept of the CNLLE is to examine the distance between two states that initially share analogous external forcings. Thus, the first step is to employ the LDAs algorithm to find the analogous initial states based on the data of external forcings. Thus, let ${\bf{y}}\left( {{t_i}} \right)$ be a reference state at time ${t_i}$, representing the variation around the black dots in Fig. 1, and ${\bf{y}}\left( {{t_j}} \right)$ be a state at another time in the observed data. To ensure analogous initial states are not merely due to persistence of time series, ${\bf{y}}\left( {{t_j}} \right)$ must meet the requirement $\left| {{t_j} - {t_i}} \right| > {t_{\mathrm{D}}}$, where ${t_{\mathrm{D}}}$ is the time taken for autocorrelations of all the variables in Y to drop to around 0 (Li and Ding, 2011). The distance between ${\bf{y}}\left( {{t_i}} \right)$ and ${\bf{y}}\left( {{t_j}} \right)$ is given by

where ${d_{{{\mathrm{in}}} ,{\text{x|y}}}}$ is the initial distance. To ensure the similarity of the external forcings not only at the initial stage but also in a short subsequent period, the initial evolutionary distance ${d_{{{\mathrm{ie}}} ,{\text{x|y}}}}$ between the trajectories of the external forcings within a short time is considered, which is given by

where K is number of samples over the evolutionary interval $\tau $[$\tau = K\Delta $, in which $\Delta $ is the sampling interval of the time series (i.e., $\Delta = {t_i} - {t_{i - 1}}$)] and $\tau $ is dependent on the persistence of the variables of interest (Li and Ding, 2011; Zhong et al., 2018a, b). The total distance ${d_{{\text{t,x|y}}}}$ is defined by the average of $ {d_{{{\mathrm{i}}} {\text{n,x|y}}}} $ and ${d_{{{\mathrm{ie}}} ,{\text{x|y}}}}$:

The ${d_{{\text{t,x|y}}}}$ for each state ${\bf{y}}\left( {{t_j}} \right)$ that satisfies the mentioned requirement can be calculated following the above procedures. If the ${d_{{\text{t,x|y}}}}$ minimizes at ${t_{k,{\text{x|y}}}}$, it is likely that the corresponding state ${\bf{y}}({t_{k,{\text{x|y}}}})$ is the analogous initial state of ${\bf{y}}\left( {{t_i}} \right)$. Thus, there is minimum uncertainty in the external forcings at these two states. Accordingly, the initial error of X between the reference trajectory and analogous trajectory is given by

where ${L_{{\text{x|y}}}}\left( {{t_0}} \right)$ is the initial error of X under the influence of external forcings (the dashed line in Fig. 1). In dynamical systems theory, if time-series variables (such as X and Y) belong to the same dynamic system, they share a common attractor manifold and the state of the one variable can be detected from the other (Sugihara et al., 2012), which means that the variation of X will follow that of Y to some extent. Therefore, analogous external forcings between the two states will lead to a similar distribution of the affected variables that corresponds to both states. In the CNLLE method, the LDAs algorithm is applied to identify analogous states of the external forcings so that the magnitude and growth of the initial error ${L_{{\text{x|y}}}}\left( {{t_0}} \right)$ between the two states will be constrained. This is different from the NLLE, which employs the LDAs algorithm to directly identify the most analogous state of the affected variables themselves (Li and Ding, 2011).

Given that the external forcings of the state at ${t_{k,{\text{x|y}}}}$ are analogous to those of the reference state ${\bf{x}}\left( {{t_i}} \right)$, the evolution of error ${L_{{\text{x|y}}}}\left( {{t_n}} \right)$ between the two states over a time interval ${t_n}$ is given by

In Fig. 1, the ${L_{{\text{x|y}}}}\left( {{t_n}} \right)$ is represented by dashed lines, where ${t_n}$(${t_n} = n \times {{\Delta }},n = 0,1,2, \cdots ,T$, where T is the number of time steps) is the evolution time, ${t_{i + n}} = {t_i} + {t_n}$, and ${t_{\left( {k + n} \right){\text{,x|y}}}} = {t_{k,{\text{x|y}}}} + {t_n}$. The growth rate of the initial error at reference state ${t_i}$ over a time interval ${t_n}$ is given by

The above steps are repeated to calculate the error growth rate for each reference state ${\bf{y}}\left( {{t_i}} \right) $ $(i = 1, 2, \cdots , M - T)$. Then, the error growth rates of all reference states are averaged to obtain the mean CNLLE:

where $\overline {{\xi _{{\text{x|y}}}}} \left( {{t_n}} \right)$ corresponds to the theoretical ${\lambda _{{\text{x|y}}}}$ in Eq. (3). This term indicates the mean divergence rate of X trajectories when it is constrained by analogous external forcings.

To visually show the error growth, the mean relative RGIE $\ln \left( {\overline {{{{\Phi }}_{{\text{x|y}}}}} } \right)$ is given by

By investigating the evolution of $\ln \left( {\overline {{{{\Phi }}_{{\text{x|y}}}}} } \right)$ with increasing ${t_n}$, the impact of external forcings on the predictability limit of X can be determined by the time at which it surpasses the large error threshold. In this study, the threshold of the error growth is set at 99% of the saturation level in the theoretical model. For the actual observational data, which exhibit large error fluctuations, a threshold of 97% of the saturation level is chosen.

As previously mentioned, the NLLE can also be estimated by following the above steps if Y is replaced by X to search analogous states in step 1, and let the mean NLLE be given by

where ${\xi _{k,{\text{x}}}}\left( {{t_n}} \right)$ is the growth rate of small initial errors at each reference state ${\bf{x}}\left( {{t_k}} \right)$. The mean RGIE of the NLLE over a period ${t_n}$ is given by

Similarly, the time when the RGIE ${\text{ln}}\left( {\overline {{{{\Phi }}_{\text{x}}}} } \right)$ grows close to its saturation is considered as the predictability limit (Li and Ding, 2011). However, comparing the RGIE of the CNLLE [Eq. (12)] and that of the NLLE [Eq. (14)], it is evident that the evolution of the REIG for both the CNLLE and NLLE methed is determined by two components: the initial error $ {L_{\text{x}}}\left( {{t_0}} \right) $ and the evolution error ${L_{\text{x}}}\left( {{t_n}} \right)$ for the NLLE. Analogously, for the CNLLE, these components are respectively represented as ${L_{{\text{x|y}}}}\left( {{t_0}} \right)$ and ${L_{{\text{x|y}}}}\left( {{t_n}} \right)$. Obviously, the error evolution of the ${L_{\text{x}}}\left( {{t_n}} \right)$ is different from that of the $ {L_{{\text{x|y}}}}\left( {{t_n}} \right) $, as the latter is constrained by external forcings, while the former is not. However, they will converge the same saturation value as they belong to the same system (Ding et al., 2008). That is,

On the other hand, the initial errors ${L_{{\text{x|y}}}}\left( {{t_0}} \right)$ and $ {L_{\text{x}}}\left( {{t_0}} \right) $ are different due to the variables X and Y being used in the LDAs method, respectively. This difference makes a direct comparison of error growth between the CNLLE and NLLE inappropriate. To overcome this issue, the initial error is replaced by that of the initial error of the NLLE in the calculation of the RGIE of the CNLLE. That is, ${L_{{\text{x|y}}}}\left( {{t_0}} \right)$ in Eq. (12) is replaced by the value of ${L_{\text{x}}}\left( {{t_0}} \right)$ from Eq. (14), and the amended RGIE of the CNLLE is given by

Therefore, the RGIE equations of the amended CNLLE and NLLE share the same initial error, and the corresponding evolution errors are projected to converge to the same values as ${t_n}$approaches infinity. As a result, a consistent saturation level of the RGIE for both methods can be attained—$\ln \left( {\overline {{{{\Phi }}_{{\text{x|y}}}^*}} } \right) = \ln \left( {\overline {{{{\Phi }}_{\text{x}}}} } \right)\left( {{t_n} \to \infty } \right)$—thus enabling a direct comparison of RGIE as estimated by the CNLLE and NLLE methods. Consequently, it allows for the quantification of the impact of external forcings on the predictability limit.

For comparative analysis, three prominent monthly SST datasets are used for comparison: (1) version 5 of the Extended Reconstructed SST (ERSST; Huang et al., 2017) on a 2º×2º spatial grid covering the period January 1854 to December 2020; (2) version 1 of the Hadley Center Sea Ice and SST (Hadley; Rayner et al., 2003) on a 1º×1º spatial grid, covering the period January 1870 to December 2020; and (3) version 2 of the Kaplan extended SST (Kaplan; Kaplan et al., 1998) on a 5º×5º spatial grid, covering the period January 1856 to December 2020.

The atmospheric global fields of monthly mean GHT at 17 pressure levels and SLP were acquired from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis dataset (Kalnay et al., 1996) over the period 1948–2020 on a 2.5º×2.5º spatial grid. To validate the effectiveness of our method, we also use SLP data from the fifth major global reanalysis produced by ECMWF (ERA5) for the period 1979-2020, with the spatial resolution interpolated to 2.5º×2.5º (Hersbach et al., 2020). For all datasets used, the seasonal cycle of the chosen period is removed before calculating the predictability limit. Moreover, the quadratic trend is removed from the SST data to obtain the SST anomaly at each grid point (Cai et al., 2020).

To demonstrate the effectiveness of the CNLLE method in quantifying the impact of the external forcings on the predictability limit, we employ a conceptual model composed of two distinct Lorenz subsystems as a testbed (Boffetta et al., 1998). The equations of the model are given by

where the superscripts (s) and (f) represent the slow and fast dynamics, respectively. The coupling between the two systems is achieved through $ - {\varepsilon _{\mathrm{s}}}x_1^{\left( {\mathrm{f}} \right)}x_2^{\left( {\mathrm{f}} \right)} $ and $ {\varepsilon _{\mathrm{f}}}x_1^{\left( {\mathrm{f}} \right)}x_1^{\left({\mathrm{ s}} \right)} $, where $ {\varepsilon _{\mathrm{s}}} $ and ${\varepsilon _{\mathrm{f}}}$ represent the coupling parameters. The choice of the form of coupling is constrained by the physical request that the solution remains in a bounded region of the phase space (Boffetta et al., 1998). When $ {\varepsilon _{\mathrm{s}}} = {\varepsilon _{\text{f}}} = 0 $, the slow dynamics and the fast dynamics are actually two Lorenz systems (Lorenz, 1963). The Lorenz system was initially derived from thermal convection, which is related to atmospheric convection in the real climate system (Saltzman, 1962). Due to its chaotic regime, this system is widely used for atmospheric predictability. For each Lorenz model, $ x_1^{\left( {\mathrm{f}} \right)} $ and $ x_1^{\left( {\text{s}} \right)} $ are proportional to the intensity of the convective motion, while $ x_2^{\left( {\mathrm{f}} \right)} $ and $ x_2^{\left( {\mathrm{s}} \right)} $ are proportional to the temperature difference between the ascending and descending currents. The variables $ x_3^{\left( {\mathrm{f }}\right)} $ and $ x_{\text{3}}^{\left( {\text{s}} \right)} $are proportional to the distortion of the vertical temperature profile from linearity (Lorenz, 1963). The parameters $ {r_{\mathrm{s}}} $ and $ {r_{\text{f}}} $ are the Rayleigh numbers, representing the intensity of the heating (Rayleigh, 1916). The Rayleigh number also serves as a key parameter determining whether chaotic behavior occurs in the Lorenz system. With the present choice of couplings, the fast dynamics is driven by means of the effective Rayleigh number, ${r_{{\text{eff}}}} = {r_{\text{f}}} + {\varepsilon _{\mathrm{f}}}x_1^{\left( {\mathrm{s}}\right)} (t) / c$, and one can recognize in the temporal evolution the slow-varying component of the driver (see Fig. 4a). Here, our focus centers on the coupled system's two distinct time scales. Specifically, the fast system tends to fluctuate 10 times as quickly as the slow dynamics, as indicated by the relative time scale factor c = 10. By considering the fast system as the atmosphere and the slow system as the oceanic system, this model can be approximated as simulating the interactions between the fast-varying atmosphere and slow-varying (external forcings) components of the Earth system. Therefore, by comparing the result of the NLLE calculated by uncoupled (${\varepsilon _{\mathrm{s}}}$=${\varepsilon _{\mathrm{f}}}$= 0) and coupled (${\varepsilon _{\mathrm{s}}} = {10^{ - 2}}$ and ${\varepsilon _{{\mathrm{f}}} } = 10$) cases, we can theoretically assess the influence of external forcing on predictability (Peña and Kalnay, 2004; Ding and Li, 2012).

Figure 4.  (a) Evolution of variables $x_1^{\left( {\mathrm{s}} \right)}$, $x_3^{\left({\mathrm{ f}} \right)}$, and 11-point moving average $\bar x_3^{\left( {\mathrm{f}} \right)}$of the fast subsystem in the coupled model over time. (b) The noise variance, signal variance, and variance of the slow system used in the calculation of the SNR.

In the following case study, the remaining parameters are chosen as $ \sigma $= 10 and b = 8/3. The Rayleigh numbers for the slow dynamics and the fast dynamics are $ {r}_{{\mathrm{s}}}=28 $ and ${r_{{\mathrm{f}}} } = 45$, respectively (Boffetta et al., 1998). The time series of the fast and slow dynamics of the coupled Lorenz system can be obtained by using the fourth-order Runge–Kutta method with a time step of 0.005. In this paper, the results obtained from $6 \times {10^5}$ steps of numerical integration are treated as the observational series. Furthermore, we report time in terms of model time units (MTUs) in the analysis, where 1 MTU = 200∆t. Therefore, the duration of the time series is 3000 MTU.

In weather forecasting, a smaller initial error implies more predictable information, indicating a higher likelihood of accurate and reliable forecasts (DelSole, 2004). The probability distributions of the logarithmic total distances between reference states and the corresponding analogous states are shown in Fig. 2, which is calculated using the NLLE, CNLLE, and random selection methods (the error between two randomly selected states points) based on the experimental data from integration of the coupled Lorenz models. Here, we take the fast subsystem to represent the affected variables, while the slow subsystem symbolizes the external forcings. On average, the initial distance of affected variables calculated by the NLLE and CNLLE is smaller than that of random selection (Fig. 2a), suggesting the initial states calculated by the CNLLE method contain more predictability information than when randomly selected. That is because the initial states in the CNLLE are the moments of analogous external forcings, and the reduced initial error can be attributed to these external forcings constraining the variability of the affected variables. Essentially, this illustrates the controlling influence of the external forcings. Moreover, when the NLLE method is applied to the fast subsystems, it is observed that the total distance of the slow subsystem is smaller than that obtained by the randomly selected method (Fig. 2b). This is because of the two-way feedback between the fast and slow subsystems, which allows the fast system to influence the variation of the slow subsystem.

Figure 2.  Probability distributions of the total distance ${d_t}$ for the (a) fast subsystem and (b) slow subsystem based on the NLLE (light blue line), CNLLE (orange line), and random selection (pink line) methods. The dashed lines indicate the positions of the average location, with their values given in parentheses.

In addition, the initial error of the affected variables, as calculated by the NLLE, is smaller than that calculated by the CNLLE (Fig. 2a), whereas the initial error of external forcings in the NLLE is larger than that in the CNLLE (Fig. 2b). This suggests that the initial states of analogous trajectories in the NLLE method contain more internally predictable information but less externally predictable information. Conversely, similar trajectories in the CNLLE method indicate more externally predictable information but relatively less internally predictable information. As such, the NLLE quantifies the predictability limit arising from the initial conditions, while the CNLLE serves to quantify the external-forcings-induced predictability limit.

Besides the initial error, the error growth rate is another crucial factor in determining the limit of predictability. Figure 3a shows that the initial error, as calculated by the NLLE, increases rapidly with time and saturates at the time of 1.67 MTU, which is close to the predictability limit (2.19 MTU) of the uncoupled fast subsystem. Such a large error growth rate is primarily due to the persistence of the initial information, unconstrained by external forcings. However, the error growth rate calculated by the CNLLE is smaller and its error growth curve resembles the slower error growth of the coupled fast subsystem calculated by the theory (Fig. 3b). The predictability limit obtained by the CNLLE is 5.84 MTU, which is longer than that obtained by the NLLE and close to the theoretical predictability limit (6.96 MTU) of the fast subsystem in the coupled model. This implies that the CNLLE provides a method to quantify second-kind predictability by capturing the error growth rate modulated by external forcings.

Figure 3.  (a) RGIE of the fast subsystem as calculated theoretically for the coupled (blue) and uncoupled (cyan) models, and estimated by the NLLE (orange) and CNLLE (red) from the experimental data. (b) A partial enlargement of the area between the dotted lines in (a), with the predictability limit (unit: MTUs) given in parentheses and marked as dots. The black lines represent the saturation value for coupled and uncoupled models, respectively.

To deepen our understanding of the CNLLE method, a comparison was made with the SNR method, another widely used method to evaluate the influence of external forcings on predictability from experimental data (Jones, 1975; Madden, 1976). Here, given that the slow dynamics exhibit a timescale that is 10 times longer than that of the fast dynamics and to compare variations on the same time scale for both fast and slow subsystems, the 11-point moving average is used to extract the predictable signal of the fast subsystem (see Appendix). Figure 4a shows that the low-frequency component (the signal of the SNR) of the fast subsystem is similar to the variation of a slow subsystem component. The variances of the signal and the slow subsystem are consistent and exceed the variance of the noise (Fig. 4b). In fact, the calculated SNR is 1.12 (larger than 1), indicating that the fast subsystem has high potential predictability that can be attributed to the contribution of the slow subsystem. However, due to the linear assumption, the SNR method may underestimate the contribution of the external forcings to the predictability (Zhang and Kirtman, 2019). Therefore, the SNR is smaller than the theoretical ratio (3.18), which is defined as the RPL of the coupled system to that of the uncoupled system. In contrast, the RPL estimated from the time series by the CNLLE to that of the NLLE is around 3.50, which is close to the theoretical ratio. Therefore, the CNLLE method can not only quantify the predictability limit due to the external forcings, but in combination with the NLLE, it is possible to calculate the RPL—a metric that can then be utilized to assess the potential for a long-term predictability limit. Specifically, the RPL larger than 1 indicates that the external forcings have the potential to increase the predictability limit. For an RPL around or less than 1, the suggestion is that the external forcings may not significantly benefit the long-term predictability limit.

ENSO, arising from interactions between the atmosphere and ocean, represents the most robust seasonal-to-interannual predictable signal in the tropical Pacific. It has significant impacts on the global ocean and atmosphere through the atmospheric bridge effect (Cai et al., 2011; Yang et al., 2018; Kug et al., 2020). In this section, the CNLLE method is applied to quantify the impact of ENSO, considered as an external forcing, on oceanic and atmospheric predictability from observational datasets. We use the Niño-3.4 index as the observed external forcing signal to identify the LDAs of ENSO (Bamston et al., 1997).

A single atmospheric variable at a given grid point is denoted as ${x_{i,j}}$, where $i$ and $j$ denote the latitudinal and meridional grid point positions, respectively. Every single grid point in the real-world atmosphere and ocean fields is closely related to the neighboring grid points (Li and Ding, 2011). Thus, the information of the eight nearest grid points around the grid point $\left( {i,j} \right)$ is considered when searching for the local analogs of ${x_{i,j}}$. For example, in Eq. (9) the distance is calculated by nine variables (M = 9).

Figure 5 shows the spatial distribution of the predictability limit of the SST field across different datasets calculated by the NLLE and CNLLE methods. The ENSO-induced predictability limit (calculated by the CNLLE) shows a distribution similar to the ENSO correlation pattern (the area of high correlation coefficients between the Niño-3.4 index and SST fields; Figs. 5a, c and e) for the three datasets. In comparison, the spatial distributions of the SST predictability limit calculated by the NLLE method is relatively inhomogeneous, but resembles that calculated by the CNLLE (Figs. 5b, d and f). This suggests that ENSO is the main predictability source of SST. Specifically, the regions of the ENSO-induced high predictability limit for SST include the central and eastern tropical Pacific, the tropical Indian Ocean, and the tropical western Atlantic, with the maximum exceeding one year. Moreover, the generally consistent distribution of the ENSO-driven predictability limit among the three different observational datasets highlights the reliability of the results.

Figure 5.  Spatial distribution of the predictability limit (PL; units: months) calculated by the CNLLE (right column) and NLLE (left column) methods based on the (a, b) ERSST, (c, d) Hadley and (e, f) Kaplan data. The contour lines in (a, c, e) represent the correlation coefficients between the SST and the ENSO index. Positive and negative correlations are shown by red solid and blue dashed contours, respectively, with an interval of 0.3, while a correlation coefficient of 0 is represented in black.

Figures 6a, c, and e show the differences in the predictability limit calculated by the two methods. It is observed that ENSO can extend the prediction horizon of the Indian Ocean’s SST by four months. This is consistent with previous research, demonstrating that the correlation is highest when ENSO leads the Indian Ocean basin mode by about four to five months (Zhang et al., 2021). ENSO can also extend the prediction horizon in the tropical western Atlantic by up to two months compared to that calculated by the NLLE. However, in the central and eastern Pacific Ocean, the predictability limit calculated by the CNLLE is at most one month longer than that calculated by the NLLE, with the RPL being around one (Figs. 6b, d and f). This suggests that the predictability limit of the region arising from ENSO is equivalent to that arising from its initial conditions. This equality is because the mean SST in the area is used to define the Niño-3.4 index, which represents the ENSO signal. Outside the tropics, ENSO makes little contribution to predictability, other than along the two great circle routes: the Pacific–North American (PNA) and Pacific–South American teleconnection (PSA) patterns (Wang, 2019). Figures 6b, d and f show that the RPL of the Southern Hemisphere is larger than that of the Northern Hemisphere, suggesting there is an asymmetric contribution of ENSO to the limit of predictability between the hemispheres. This can be attributed to the inconsistency in the PNA and PSA response to ENSO. Previous studies have shown that the ENSO-induced PNA exists only in winter and spring, while the PSA can be present year-round (Liu et al., 2012). Furthermore, as ENSO typically peaks in the austral summer, during the season, the warmest SST and the associated ascending branch of the Hadley cell shift to the Southern Hemisphere (Guo and Tan, 2018). These enhanced ocean-air interactions in the Southern Hemisphere make the ENSO signal clearer.

Figure 6.  Spatial distribution of the difference in the predictability limit (PL; left column, unit: months) and RPL (right column) distribution of SST using the CNLLE and NLLE methods based on the (a, b) ERSST, (c, d) Hadley and (e, f) Kaplan data. The contour lines represent the correlation coefficients between the SST fields and the ENSO index. Positive and negative correlations are shown by red solid and blue dashed contours respectively, with an interval of 0.3, while a correlation coefficient of 0 is represented in black.

To investigate the impact of ENSO on the atmospheric predictability limit, the spatial distributions of the predictability limit of the SLP fields are calculated by the CNLLE and NLLE methods (Figs. 7a–d). Figures 7a and b show that the predictable signal from ENSO for the SLP fields is mainly concentrated in the tropics, providing seven to ten months predictability limit for areas such as the Maritime Continent, the central and eastern tropical Pacific, and the tropical Atlantic. These regions exhibit high correlation coefficients between the SLP fields and the ENSO index, while the ENSO-induced predictability limit of SLP fields over the central and eastern tropical Pacific is smaller than that in the western tropical Pacific Ocean. This discrepancy can be attributed to the warmer background SST in the tropical western Pacific Ocean compared to the central-east Pacific Ocean (Yan et al., 1992; Yang et al., 2023). The higher SST promotes stronger convection, and indicates more intense ocean-air interaction in the tropical western Pacific Ocean, making ENSO signals easier to detect in that region (Fu et al., 1994). The spatial distribution of the predictability limit calculated by the NLLE (Figs. 7c and d) resembles that calculated by the CNLLE. However, it exhibits a shorter predictability limit resulting in the RPL being greater than 1 in these areas (Figs. 7e and f). This suggests that in those regions where there are relatively strong correlations between the ENSO index and SLP fields, ENSO can be a significant predictable signal. However, there is no clear ENSO-driven predictability limit over the PNA and PSA regions in the midlatitude SLP fields, especially when comparing the distribution of the ENSO-driven predictability limit of SST fields (Figs. 5a, c and e). The lack of a clear impact over the PNA and PSA regions in the midlatitude SLP fields may be partly attributable to the short duration of observations of SLP fields. Additionally, the correlation between the ENSO index and HGT fields is lower than with SST fields, suggesting that midlatitude SLP is influenced by other factors like the subtropical jet and Aleutian low system (Overland et al., 1999; Liu et al., 2012). Moreover, the movement of the atmosphere is relatively stochastic, making ENSO impacts less apparent in the SLP fields (Wang et al., 2021).

Figure 7.  Spatial distribution of the predictability limit (PL; units: months) of SLP calculated by the (a, b) CNLLE and (c, d) NLLE methods, along with (e, f) the RPL of SLP based on the NCEP (left column) and ERA5 (right column) data. The contours lines represent the correlations between the SLP and the ENSO index. Positive and negative correlations are shown by red solid and blue dashed contours, respectively, with an interval of 0.2.

Figure 8 displays the spatial distributions of the predictability limit of the GHT fields at different levels, calculated by both the CNLLE and NLLE methods. The results indicate that the spatial patterns of the ENSO-induced predictability limit for the GHT at 850 hPa are very similar to those for the SLP (Fig. 8a). Specifically, the ENSO-induced predictability limit is six to eight months over the Indian Ocean, six to twelve months over the Maritime Continent and adjacent areas, four to eight months in the central Pacific, and six to ten months in the tropical Atlantic. For the upper atmospheric layers, the ENSO-induced predictability limit of the GHT appears to have a zonal distribution over the equatorial regions between 20°N and 20°S, with the peak predictability limit reaching about one year (Figs. 8c and e). Additionally, it is evident that ENSO can enhance the predictability of subtropical GHT in regions where the correlation coefficients between the ENSO index and the GHT fields show significant values. Compared to the predictability limit calculated by the NLLE (Figs. 8b, d and f), the spatial distribution influenced by ENSO is analogous but exhibits greater uniformity. In certain areas where ENSO is acting as an external driver, it has the potential to add up to four months to the prediction horizon derived from the initial GHT conditions (Figs. 9a, c and e). Over most of the tropical regions, the RPL values are larger than 1 (Figs. 9b, d and f), indicating that ENSO can contribute to long-term atmospheric predictability in these areas. Although we can detect ENSO’s predictable signal over the PNA and PSA regions in the midlatitudes, with predictability limits of driven by ENSO of two to four months, it is shorter than that driven by the initial condition. This may reflect that ENSO’s impacts may not be the major predictable signal; other weather and climate systems, such as the subtropical jet and Aleutian low, also have significant impacts on the movement in those regions (Overland et al., 1999; Liu et al., 2012).

Figure 8.  Spatial distributions of the GHT predictability limit (PL; units: months) calculated by the CNLLE (left column) and NLLE (right column) methods for the (a, b) 850 hPa, (c, d) 500 hPa, and (e, f) 300 hPa. The contour lines represent the correlation coefficients between the GHT fields and ENSO index, where positive and negative correlations are denoted by red solid and blue dashed contours, respectively, with an interval of 0.2.

Figure 9.  The difference in the predictability limit (left column) and RPL (right column) distribution of GHT using the CNLLE and NLLE methods at (a, b) 850 hPa, (c, d) 500 hPa, and (e, f) 300 hPa. The contour lines represent the correlation coefficients between the GHT fields and ENSO index, where positive and negative correlations are denoted by red solid and blue dashed contours, respectively, with an interval of 0.2.

Figure 10 shows the vertical distributions of the GHT predictability limit. The high predictability limit obtained by the NLLE is mainly concentrated between 100 and 50 hPa (up to 16 months), with the predictability limit in the northern and southern subtropical regions being lower than in the tropical regions (Figs. 10c and d). The high predictability limit induced by ENSO is uniformly concentrated in the latitudinal direction at altitudes between 400 and 150 hPa, with a predictability limit of up to 10 months. However, near the surface, the predictability limit decreases to four to eight months. At the lower atmospheric levels of the eastern region, there is a lower predictability limit compared to the western region (Figs. 10a and e). Additionally, above the 50 hPa level, our results show there is an absence of predictability induced by ENSO (Figs 10a and b). In the longitudinal direction, the ENSO-induced predictability limit is predominantly in the tropical region and extends to upper levels where a positive correlation with ENSO is observed. Furthermore, the contribution of ENSO to GHT predictability in the southern region is greater than that in the northern region (Figs. 10b and f).

Figure 10.  Zonal-height (left column) and meridional-height (right colunm) distribution of the predictability limit (PL; unit: months) calculated by the (a, b) CNLLE, (c, d) NLLE methods and (e, f) the vertical distribution of GHT’s RPL. Note that the zonal-height cross-section plots are averaged over 20°S–20°N and meridional-height cross-section plots are averaged over 90°W–180°. Contour lines represent the correlation coefficient between the GHT fields and the ENSO index, where positive and negative correlations are denoted by red solid and blue dashed contours, respectively, with an interval of 0.2.

In this work, we develop a specific algorithm of the CNLLE to investigate second-kind predictability based on the NLLE method. Different from the NLLE, which mainly focuses on the predictable information arising from the initial conditions, the CNLLE algorithm captures the slow error growth constrained by external forcings, and can therefore quantify the predictability limit induced by external forcings. To assess whether external forcings can improve long-term predictability, we introduce the RPL, which is defined as the ratio of the predictability limit calculated by the CNLLE to that calculated by the NLLE. Similar to the SNR method, if the RPL is significantly greater than 1, the external forcings are potential predictors for long-term forecasting. However, if the RPL is significantly less than 1, the external forcings may not contribute to the long-term predictability. If the RPL is around 1, the external forcings and affected variables are inextricably intertwined and difficult to separate. However, the effectiveness is only examined in a coupled Lorenz model. There exists another form of the coupled model (Peña and Kalnay, 2004), and further work is required to assess its sensitivity to the structure of the coupled equations.

Then, the spatial distributions of the ENSO-induced predictability limit for monthly SST, SLP, and GHT are estimated using the CNLLE method, and compared to the predictability limit arising from the initial-conditions calculated by the NLLE method. We find that ENSO contributes to the predictability limit of atmosphere and ocean fields in regions with relatively high correlation coefficients between those fields and the ENSO index. For SST, ENSO can extend the predictability limit arising from the initial conditions by up to four months in the tropical Indian Ocean and by one to two months in the mid- to high-latitude teleconnection areas of the subtropics. For the atmospheric fields, the spatial distribution of the ENSO-induced predictability limit of the SLP appears similar to that of the low-level GHT, mainly concentrated in the eastern and western tropical Pacific, with a predictability limit of over eight months, which is about four months longer than that arising from initial conditions. For upper-level GHT, the patterns of the ENSO-induced predictability limit exhibit a zonal distribution over the tropical region within the troposphere below 50 hPa.

ENSO can impact the predictability limit of both the ocean and atmosphere in the midlatitudes, but the magnitude of the impact differs. ENSO provides significant contributions to the predictability limit of SST fields in the midlatitudes along the PNA and PSA ENSO-teleconnection routes. However, the impact of ENSO’s predictable signal on atmospheric fields is weakly observed in the midlatitudes of both the Northern and Southern Hemispheres, especially in the SLP fields. Several factors could contribute to this. Firstly, the quality and length of observational data are crucial factors in determining whether the CNLLE can identify true analogous trajectories for the reference trajectories. In this work, the length of atmospheric observation data is relatively shorter than those of the ocean, which may result in the ENSO signal not being properly detected. Secondly, the correlation between the ENSO index and HGT/SLP fields is lower than with SST fields, indicating that the weather in the midlatitudes may be influenced by factors like the subtropical jet and Aleutian low system (Overland et al., 1999; Liu et al., 2012; Zhang et al., 2023), making ENSO impacts less apparent. Additionally, compared to the ocean's large thermal inertia, the movement of the atmosphere is relatively stochastic (Constantinou and Hogg, 2021; Huang et al., 2023). For instance, it has been found that there is an asymmetric response of the PNA to La Niña and El Niño, representing the cold and warm phases of ENSO (Wang et al., 2021), respectively, suggesting inconsistent ENSO influence in the midlatitudes. Hence, consistently detecting the low-frequency variation of ENSO signals in the atmosphere proves to be challenging.

This study is a preliminary application of the CNLLE in investigating the impact of ENSO to atmospheric and oceanic predictability based on observational data. Our findings suggest that ENSO could be a crucial long-term predictor for regions where the value of RPL exceeds 1, which might provide a guidance for long-term statistical weather and climate prediction. However, knowledge about the ENSO-induced predictability limit of oceanic and atmospheric fields based on observational data is not sufficient to improve the numerical model. Further research is needed to comprehend the practical predictability of numerical weather forecasting models, revealing and reducing the gap between the predictable horizon of state-of-the-art models and their intrinsic predictability limit. Hence, the next focus is to fully explore how predictable signals propagate from the ocean to the atmosphere, especially understanding the physical processes involved. For instance, to investigate why ENSO's predictable signal is clearly detected in SST fields but not in the atmosphere over midlatitudes through the PNA and PSA routes. Such knowledge will inform the development of coupled models and reduce model bias in long-term weather forecasting. Moreover, it is important to recognize that, in addition to ENSO, there are other types of external forcings (Bellucci et al., 2015) within the Earth system, such as greenhouse gas emissions and volcanic eruptions, exerting significant impacts on the climate system. Unlike ENSO, which is a phenomenon involving a two-way feedback between the ocean and atmosphere, these external forcings are less impacted by atmospheric motion. Further exploration is necessary to assess the generalization ability of the CNLLE method in estimating the impacts of those one-way external forcings on the predictability limit, both in conceptual models and real climate systems. A more generalizable approach can be used to identify the most important sources of predictable signals from these interwoven external forcings, which will provide more valuable insights for weather forecasting.

Acknowledgements. This research was jointly supported by the National Natural Science Foundation of China (Grant Nos. 42225501 and 42105059), and the National Key Scientific and Technological Infrastructure project “Earth System Numerical Simulation Facility” (EarthLab).

The SNR method has been extensively used in climate science to measure the strength of a signal compared to the level of noise present in the data. For example, let ${x_t}$ represents the time series of a variable, which can be decomposed as follows:

where $\mu $ is the average value, and${x_{{\mathrm{st}}}}$ and ${x_{{\mathrm{ft}}}}$ donate the slow-frequency signal and high-frequency noise, respectively. The assumption of small nonlinear interaction between the slow-frequency signal and high-frequency noise can be drawn because of their distinct timescales, and therefore the covariance between them is ${\mathrm{Cov}}\left( {{x_{{\mathrm{ft}}}},{x_{{\mathrm{st}}}}} \right) = 0$. By taking the variance of both sides of Eq. (A1) and applying the formula for variance calculation, we can derive:

where $\mu $ is constant, so ${\mathrm{Var}}\left( \mu \right) = 0$, ${\mathrm{Cov}}\left( {\mu ,{x_{{\mathrm{st}}}}} \right) = 0$, and ${\mathrm{Cov}}\left( {\mu ,{x_{{\mathrm{ft}}}}} \right) = 0$. Then we can obtain:

where ${\mathrm{Var}}\left( {{x_{{\mathrm{st}}}}} \right)$ is the variance of the predictable signal, while $ {\mathrm{Var}}\left( {{x_{{\mathrm{ft}}}}} \right) $ is the variance of the unpredictable noise. Due to the relatively weak nonlinear interaction between the signal and noise, represented by the term ${\mathrm{Cov}}\left( {{x_{{\mathrm{st}}}},{x_{{\mathrm{ft}}}}} \right)$, it can be considered negligible. Consequently, Only the variance of the signal is deemed potentially predictable, so potential predictability can be measured by the ratio of the variance of the signal to the variance of the noise. Thus, the SNR can be defined as:

Typically, when the SNR is greater than 1, it suggests that the contribution of the signal to the variability of the system is larger than that of the noise, and the system can be considered highly predictable. Conversely, when the SNR is less than 1, the noise is greater than that of the signal, leading to lower predictability of the system. In this study, the 11-point moving average is used to extract the signal component (${x_{{\mathrm{st}}}}$) of the coupled fast system, while the high-frequency noise component (${x_{{\mathrm{ft}}}}$) can be obtained using ${x_{{\mathrm{ft}}}} = {x_t} + \mu - {x_{{\mathrm{st}}}}$.