Freddy Bouchet
email: freddy.bouchet@ens-lyon.fr

Committor Functions for Climate Phenomena at the Predictability Margin: The Example of El Niño Southern Oscillation in the Jin and Timmermann Model

Many atmosphere and climate phenomena lie in the gray zone between weather and climate: they are not amenable to deterministic forecast, but they still depend on the initial condition. A natural example is medium-range forecasting, which is inherently probabilistic because it lies beyond the deterministic predictability time of the atmosphere, but for which statistically significant prediction can be made which depend on the current state of the system. Similarly, one may ask the probability of occurrence of an El Niño event several months ahead of time. We introduce a quantity which corresponds precisely to this type of prediction problem: the committor function is the probability that an event takes place within a given time window, as a function of the initial condition. We compute it in the case of a low-dimensional stochastic model for El-Niño, the Jin and Timmermann model. In this context, we show that the ability to predict the probability of occurrence of the event of interest may differ strongly depending on the initial state. The main result is the new distinction between probabilistic predictability (when the committor function is smooth and probability can be computed which does not depend sensitively on the initial condition) and probabilistic unpredictability (when the committor function depends sensitively on the initial condition). We also demonstrate that the Jin and Timmermann model might be the first example of a stochastic differential equation with weak noise for which transition between attractors do not follow the Arrhenius law, which is expected based on large deviation theory and generic hypothesis.

Significance Statement

A key problem for atmospheric and climate phenomena is to predict events beyond the timescale over which deterministic weather forecast is possible. In a simple model of El Nino, we demonstrate the existence of two regimes, depending on initial conditions. For initial conditions in the "probabilistic predictability" regime, the system is unpredictable deterministically because of chaos, but the probability of occurrence of the event can still be predicted because it depends only weakly on the initial condition. In the "probabilistic unpredictability" regime, even predicting probabilities is difficult, because the probability depends strongly on initial conditions. These new concepts of probabilistic predictability and unpredictability should be key in understanding the predictability potential for rare events in climate problems, and in other complex dynamics.

Introduction

It has long become clear that statistics and probability are the natural languages for climate: for given boundary conditions, there is a typical state (or several, in case of bimodality), the climatology, and fluctuations around typical conditions, referred to as climate variability, involving various time and space scales. At first sight, this kind of description may seem orthogonal to the problem of weather forecasting, which consists in predicting the exact state of the atmosphere at a given future time. However, notwithstanding the use of probabilities in numerical weather forecasting for uncertainty quantification, these two approaches meet in several areas of current climate research [START_REF] Jin | An equatorial ocean recharge paradigm for ENSO. Part II: A stripped-down coupled model[END_REF][START_REF] Dijkstra | Nonlinear climate dynamics[END_REF][START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF]. For instance, we are often interested in predicting the occurrence of specific fluctuations of the climate system, be it a given mode of climate variability, such as the El Niño Southern Oscillation (ENSO) [START_REF] Philander | El Niño, La Niña, and the Southern Oscillation[END_REF]), regime changes [START_REF] Tantet | An early warning indicator for atmospheric blocking events using transfer operators[END_REF], or rare events such as heat waves [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF], sudden stratospheric warming, cold spells, extreme precipitations, or any other event of importance. All these events have a probability of occurring any given year, i.e. with respect to climatological conditions, but one may also be interested in their probability of occurrence conditioned on the state of the climate system at the time of the prediction. For instance, given the global impact of events like ENSO, much efforts have focused on developing methods to forecast it several months in advance (e.g. [START_REF] Latif | A review of ENSO prediction studies[END_REF][START_REF] Clarke | An introduction to the dynamics of El Niño & the Southern Oscillation[END_REF]; [START_REF] Chekroun | Predicting stochastic systems by noise sampling, and application to the El Niño-Southern Oscillation[END_REF][START_REF] Ludescher | Very early warning of next El Niño[END_REF][START_REF] Feng | Climate network stability measures of El Niño variability[END_REF][START_REF] Nooteboom | Using Network Theory and Machine Learning to predict El Niño[END_REF]). Similarly, one may want to estimate the probability of occurrence of a summer drought based on soil moisture in the spring, the probability of occurrence of a heat wave a few weeks in advance, based on the observed atmospheric circulation, or the probability of an extreme hurricane season, based on sea surface temperature. Such forecasts are extremely challenging, but would be rewarded with proportionally large benefits, given the socio-economic impact of these events at the local and global scales, especially in a climate change context [START_REF] Aghakouchak | Extremes in a changing climate: detection, analysis and uncertainty[END_REF]2 Coumou and Rahmstorf 2012;[START_REF] Field | Managing the risks of extreme events and disasters to advance climate change adaptation: special report of the intergovernmental panel on climate change[END_REF][START_REF] Herring | Explaining extreme events of 2013 from a climate perspective[END_REF]). While it is not clear that this may be reliably achieved for all the above exemples, due to their different physical nature, conceptually all these events fall in the same class of prediction problems. The goal of this paper is to discuss the mathematical structure of such climate prediction problems.

Indeed, the mathematical structure of weather forecasting is quite clear: it consists in solving an initial value problem (IVP). Given an initial condition 𝑥 0 belonging to an appropriate phase space, we are interested in computing the trajectory 𝑥 𝑡 = 𝜙 𝑡 𝑥 0 , where 𝜙 𝑡 is the flow of the dynamical system, encapsulating the evolution equations. For many dynamical systems, this description only holds for a finite time in practice, due to sensitive dependence on initial conditions. This limitation was already known from mathematicians in the 19th century, such as Poincaré and Hadamard. For low-dimensional chaotic dynamical systems, this time scale, up to which deterministic forecasts are relevant, corresponds the the Lyapunov time [START_REF] Castiglione | Chaos and coarse graining in statistical mechanics[END_REF]). In the atmosphere, the predictability horizon, about two weeks in practice, corresponds to the time it takes for undetectable errors at the smallest scales of the flow to contaminate the large scales [START_REF] Thompson | Uncertainty of initial state as a factor in the predictability of large scale atmospheric flow patterns[END_REF][START_REF] Novikov | On the problem of predictability of synoptic processes[END_REF][START_REF] Lorenz | The predictability of a flow which possesses many scales of motion[END_REF]. Early numerical weather prediction attempts fell short of this predictability horizon, both due to model errors and sparsely constrained initial conditions. As models improved and observational data became much denser, owing in particular to the advent of satellite observations, performance rose and skillfull forecasts are now close to the theoretical barrier [START_REF] Bauer | The quiet revolution of numerical weather prediction[END_REF]. Beyond this limit, the dynamics becomes effectively stochastic. Notwithstanding the fact that the relevant phase space may be different for climate dynamics over geological time scales, climate therefore corresponds to the statistical properties of some stochastic process (𝑥 𝑡 ) 𝑡 >0 . Over very long times, we expect those statistical properties to be independent of the initial condition. In other words, the mathematical concept relevant for climate is the invariant measure of the system. For lack of better techniques, in practice we still compute these properties by averaging over long times and over realizations using ensembles of trajectories obtained by numerical integration of climate models. In any case, the invariant measure only describes the system for times larger than the mixing time, after which the initial condition is forgotten. However, in the applications cited above, the time scale of interest is the intermediate case for which a deterministic forecast is not relevant, but for which some information, more precise than the climate average, might be predicted. We call this range of time scales the predictability margin. We stress that the concept of predictability margin depends on the particular phenomenon under investigation.

Prediction problems at the predictability margin are of a probabilistic nature, because they are concerned with . The associated mathematical problem is an initial value problem (IVP). Long term statistical properties (beyond the mixing time) do not depend on the initial condition, and the corresponding mathematical object is the invariant measure.

In the intermediate range of timescales, which we call the predictability margin here, the appropriate mathematical concept is the committor function, which encodes the probability of a given event to occur, condition on the state of the climate system at the time of the prediction.

time scales beyond the deterministic predictability horizon of the system (e.g. the Lyapunov time). However, we stress that the Lyapunov time scale, a global quantity, is clearly not the relevant dynamical quantity for this predictability problem. By contrast, at the predictability margin, the predictability depends on the current state of the system, by definition. Then, the question is: what is the relevant mathematical concept for prediction problems at the predictability margin? The relevant mathematical concept is called the committor function [START_REF] Ren | Transition pathways in complex systems: Reaction coordinates, isocommittor surfaces, and transition tubes[END_REF][START_REF] Metzner | Illustration of transition path theory on a collection of simple examples[END_REF]. This is a very generic concept: a committor function is the probability for an event to occur in the future, as a function of the current state of the system. Committor functions have first been introduced in climate sciences in [START_REF] Lucente | Machine learning of committor functions for predicting high impact climate events[END_REF], and have been used to study sudden stratospheric warming [START_REF] Finkel | Path properties of atmospheric transitions: illustration with a low-order sudden stratospheric warming model[END_REF][START_REF] Finkel | Learning forecasts of rare stratospheric transitions from short simulations[END_REF] and to understand the flow of ocean debris [START_REF] Miron | Transition paths of marine debris and the stability of the garbage patches[END_REF]. The interest of putting a name, the committor function, to this otherwise very common and generic concept, is two-fold. First it allows to study its mathematical properties and to relate them to key concepts in dynamical systems, for instance the predictability margin, as we do in the present work. Second, it comes with specific theoretical and computational approaches to compute this important quantity. For instance, the committor function is one of the essential components of transition path theory, which is a mathematical framework for studying rare events, see for example [START_REF] Metzner | Illustration of transition path theory on a collection of simple examples[END_REF]; [START_REF] Metzner | Illustration of transition path theory on a collection of simple examples[END_REF][START_REF] Metzner | Transition path theory for markov jump processes[END_REF] and references therein. In atmosphere dynamics, a very interesting use of smart ways of estimating the committor function for a simplified model of sudden stratospheric warming is provided by [START_REF] Finkel | Learning forecasts of rare stratospheric transitions from short simulations[END_REF].

Many problems in medium-range forecasting fall within the predictability margin range; to illustrate the interest of committor functions, we will select only one example of application, the problem of ENSO prediction, using a very simple model. While, as mentioned above, many studies strive to predict the occurrence of El Niño a few months in advance, we shall address here a slightly different problem, focusing on predicting strong El Niño events on longer time scales. This is also a relevant question from the point of view of climate dynamics; while strong El Niño events have been observed almost periodically since the 1950s, with a return time around 15-20 years, historical data and paleoclimatic proxies indicate that ENSO may exhibit high variability over centennial timescales [START_REF] Cobb | El Niño/Southern Oscillation and tropical Pacific climate during the last millennium[END_REF][START_REF] Khider | Assessing el niño southern oscillation variability during the past millennium[END_REF][START_REF] Mcgregor | Inferred changes in El Niño-Southern Oscillation variance over the past six centuries[END_REF] and beyond [START_REF] Rickaby | Cool La Niña during the warmth of the Pliocene?[END_REF][START_REF] Fedorov | The Pliocene paradox (mechanisms for a permanent El Niño)[END_REF][START_REF] Cobb | Highly variable El Niño-Southern Oscillation throughout the Holocene[END_REF]). We study the dynamics of a low-dimensional stochastic model proposed to explain the decadal amplitude changes of ENSO, the Jin and Timmermann model [START_REF] Timmermann | A nonlinear mechanism for decadal El Niño amplitude changes[END_REF][START_REF] Timmermann | A nonlinear theory for El Niño bursting[END_REF]). This model is not aimed at reproducing any precise properties of the real El Niño Southern Oscillation. It is rather used as a paradigm example to introduce the concept of a committor function, and to study its main properties. This will lead us to define probabilistic predictability and unpredictability, some concepts that should be useful for other applications.

We show that probabilistic prediction at the predictability margin depends on the initial state, and that probabilistic predictability is encapsulated in the committor function. This property is analogous to classical, deterministic predictability, which is known to depend on the state of the system: some circulation patterns, such as the positive phase of the North Atlantic Oscillation (NAO), lead to improved predictability. However, we stress that deterministic and probabilistic predictability are different concepts: probabilistic predictability means that the probability of the event does not depend sensitively on the initial conditions. This corresponds to a region of phase space where the committor function has gentle variations with the initial conditions. In these areas, the event occurs with a probability 𝑝 that can be easily determined in practice because of these gentle variations. On the contrary, probabilistic unpredictability corresponds to regions of the phase space with a rough committor function. In these regions, the occurence of the event is also probabilistic. But the probability itself has very rapid variations with the initial conditions, which make the prediction highly dependent on the level of precision with which the initial condition is known. The existence of such features, and especially the new and most interesting probabilistically predictable region, should be generic for most prediction problems in climate dynamics.

This paper also discusses relations between qualitative properties of the committor function, finite time Lyapunov exponents, and the stability properties of trajectories with respect to noise perturbations. It also discuss methodological aspects for a data-based approach for the computation of committor functions.

The dynamics of the Jin and Timmermann model, when perturbed by a weak noise, is characterized by rare transitions between a limit cycle and a strange attractor [START_REF] Roberts | Mixed-mode oscillations of El Niño-Southern Oscillation[END_REF][START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF]. Based on large deviation theory, and with generic hypothesis, the average transition time E[𝜏 𝑐 ] to see such transitions is expected to change following an Arrhenius law:

E[𝜏 𝑐 ] 𝜎→0 𝐴 exp Δ𝑉/𝜎 2 ,
where 𝜎 is the noise amplitude. Using numerical simulations, we demonstrate that the Jin and Timmermann transition times do not follow the expected Arrhenius law for a very large range of small noise amplitudes. We conjecture that this very interesting phenomenon might be the first observed counterexample to the expected generic result, for deterministic dynamics perturbed by weak noises. We argue that this is related to the intricated entanglement between the basins of attraction of the limit cycle and the strange attractor.

The paper is organized as follows: in Sec. 2 we define the Jin and Timmermann model [START_REF] Timmermann | A nonlinear mechanism for decadal El Niño amplitude changes[END_REF][START_REF] Timmermann | A nonlinear theory for El Niño bursting[END_REF]. In this model, the occurrences of strong ENSO events correspond to noiseinduced transitions between a strange attractor and a limit cycle [START_REF] Roberts | Mixed-mode oscillations of El Niño-Southern Oscillation[END_REF][START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF]). We study in Sec. 3 the statistics of such transitions, and we show that they do not obey an Arrhenius law. Finally, in Sec. 4 we introduce the committor function, we compute it for the Jin-Timmermann model, and characterize the regions of the phase space with qualitatively different predictability properties. In the regime of intermediate noise amplitude, at the predictability margin, we delineate four regions (see Fig. 7): two regions of deterministic predictability (where the event occurs with probability 0 or 1), one probabilistically predictable region (where a value of the probability 0 < 𝑞 < 1 can clearly be predicted with very mild dependence with respect to initial condition), and finally a region which is unpredictable in practice, because the strong dependence with respect to the initial condition prevents any precise prediction, either deterministic or probabilistic.

The Jin and Timmermann model

El Niño Southern Oscillation (ENSO) is one of the most important mode of climate variability at the interannual time scales [START_REF] Philander | El Niño, La Niña, and the Southern Oscillation[END_REF]. El Niño events consist in an increase of the Sea Surface Temperature in the eastern equatorial Pacific Ocean, leading at the local scale to reduced thermocline depth, reduced upwelling and reduced nutrient supply, thereby affecting marine life. Such events are correlated with a reorganization of the Walker circulation in the atmosphere, known as the Southern Oscillation.

The global phenomenon, referred to as ENSO, has major impacts all over the world. However, the nonlinear coupled atmosphere-ocean dynamics of ENSO makes it very difficult to predict (e.g. [START_REF] Mcphaden | The curious case of the el niño that never happened: a perspective from 40 years of progress in climate research and forecasting[END_REF]. Models of various complexities have been constructed to capture the dynamics of El Niño at different levels of realism [START_REF] Clarke | An introduction to the dynamics of El Niño & the Southern Oscillation[END_REF][START_REF] Sarachik | The El Nino-Southern Oscillation phenomenon[END_REF]. In order to introduce and illustrate the concept of committor function in the simplest possible framework, we shall consider here one of the most idealized models, consisting of a low-dimensional stochastic process. This simple dynamical model, introduced by Jin (1997a,b), accounts for the recharge-discharge mechanism which is at the basis of ENSO. This model was later extended by [START_REF] Timmermann | A nonlinear theory for El Niño bursting[END_REF] and was related to the decadal amplitude changes of ENSO [START_REF] Timmermann | A nonlinear mechanism for decadal El Niño amplitude changes[END_REF]. The model describes the evolution of three variables:

1. 𝑇 1 , the Sea Surface Temperature in the western equatorial Pacific Ocean, 2. 𝑇 2 , the Sea Surface Temperature in the eastern equatorial Pacific Ocean, 3. ℎ 1 , the thermocline depth anomaly in the western Pacific.

Assuming a thermal relaxation towards a radiativeconvective temperature 𝑇 𝑟 , the equations of motion can be written as [START_REF] Timmermann | A nonlinear mechanism for decadal El Niño amplitude changes[END_REF][START_REF] Timmermann | A nonlinear theory for El Niño bursting[END_REF]:

d𝑇 1 d𝑡 = -𝛼(𝑇 1 -𝑇 𝑟 ) -𝜖 𝛽𝜏(1 -𝜎𝜂 𝑡 ) (𝑇 2 -𝑇 1 ), d𝑇 2 d𝑡 = -𝛼(𝑇 2 -𝑇 𝑟 ) + 𝜁 𝛽𝜏(1 -𝜎𝜂 𝑡 ) (𝑇 2 -𝑇 𝑠𝑢𝑏 ), dℎ 1 d𝑡 = 𝑟 -ℎ 1 - 1 2 𝑏𝐿𝜏 , (1)
where 𝜖 and 𝜁 represent the strength of the zonal and vertical advection, 𝑇 𝑠𝑢𝑏 denotes the temperature being upwelled into the mixed layer, 𝜏 represents the zonal wind stress, 𝐿 denotes the basin width, 𝑏 captures the efficiency of wind stress in driving thermocline tilt, 1/𝛼 measures a typical thermal damping timescale and 1/𝑟 is the dynamical adjustment timescale of the thermocline depth. The term 𝜂 𝑡 in the equations for temperatures is a Gaussian white noise with unit variance and the level of stochasticity is controlled by the noise amplitude 𝜎. This term takes into account the fluctuating component of wind stress. In the last equation the noise has not been considered because wave processes are filtered out in the thermocline equations of the model [START_REF] Timmermann | A nonlinear theory for El Niño bursting[END_REF].

The expressions of 𝑇 𝑠𝑢𝑏 and 𝜏 are

𝑇 𝑠𝑢𝑏 = 𝑇 𝑟 + 𝑇 𝑟 0 2 + 𝑇 𝑟 -𝑇 𝑟 0 2 tanh 𝐻 + ℎ 2 -𝑧 0 ℎ * , 𝜏 = 𝜇(𝑇 2 -𝑇 1 ) 𝛽 , (2)
where 𝑇 𝑟 0 is a reference temperature, ℎ 2 is the thermocline departure from its reference value 𝐻, 𝑧 0 represents the depth at which 𝜁 takes its characteristic value, ℎ * measures the sharpness of the thermocline. The relation between the eastern and western thermocline depth anomalies is

ℎ 2 = ℎ 1 + 𝑏𝐿𝜏. (3)
In order to study the dynamical behavior of the system it is useful to perform a change of variables from physical to dimensionless ones [START_REF] Roberts | Mixed-mode oscillations of El Niño-Southern Oscillation[END_REF]. So, we define

𝑥 = 𝑇 2 -𝑇 1 𝑇 0 , 𝑦 = 𝑇 1 -𝑇 𝑟 𝑇 0 , 𝑧 = ℎ 1 + 𝐻 -𝑧 0 ℎ * , t = 𝑡 𝑡 * ,
where 𝑇 0 = ℎ * 𝛽 𝑏𝐿 𝜇 and 𝑡 * = 𝑏𝐿 𝛽𝜁 ℎ * . After the change of variables, the equations ( 1) read

𝑥 = 𝜌𝛿(𝑥 2 -𝑎𝑥) + 𝑥 [𝑥 + 𝑦 + 𝑐 -𝑐 tanh (𝑥 + 𝑧)] -𝐷 𝑥 (𝑥, 𝑦, 𝑧)𝜂 𝑡 , 𝑦 = -𝜌𝛿(𝑥 2 + 𝑎𝑦) + 𝐷 𝑦 (𝑥, 𝑦, 𝑧)𝜂 𝑡 , 𝑧 = 𝛿 𝑘 -𝑧 - 𝑥 2 , ( 4
)
where

𝐷 𝑥 (𝑥, 𝑦, 𝑧) = [(1 + 𝜌𝛿)𝑥 2 + 𝑥𝑦 + 𝑐𝑥(1 -tanh (𝑥 + 𝑧))]𝜎, 𝐷 𝑦 (𝑥, 𝑦, 𝑧) = 𝜌𝛿𝑥 2 𝜎,
and the new control parameters 𝛿, 𝜌, 𝑐, 𝑘, and 𝑎 are defined as follows:

𝛿 = 𝑟𝑏𝐿 𝜁 𝛽ℎ * , 𝜌 = 𝜖 ℎ * 𝛽 𝑟𝑏𝐿 , 𝑎 = 𝛼𝑏𝐿 𝜖 𝛽ℎ * , 𝑐 = 𝑇 𝑟 -𝑇 𝑟 0 2𝑇 0 , 𝑘 = 𝐻 -𝑧 0 ℎ * .
The deterministic version (𝜎 = 0) of equations ( 4) was widely studied in literature. For some parameter values, the system has only one attractor, a periodic orbit. 2016) also analyzed the mechanism through which this limit cycle arises. [START_REF] Roberts | Mixed-mode oscillations of El Niño-Southern Oscillation[END_REF] defined strong El-Niño events for this model as periods in this limit cycle for which the temperature is large. Figure 2 shows a qualitative comparison of the eastern Pacific sea surface temperature anomaly for this limit cycle with the El-Niño3 index. Both the measurements and the model display positive temperature anomaly excursions with a return time of approximately 15 years. Varying the parameter 𝛿, a strange attractor emerges through a period doubling cascade, as shown by [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF].

Moreover, [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF] show that for some parameter values this strange attractor coexists with a limit cycle similar to the one mentioned in the previous paragraph (with a different period).

Following [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF], we use the parameter values [𝛿, 𝜌, 𝑐, 𝑘, 𝑎] = [0. 225423, 0.3224, 2.3952, 0.4032, 7.3939] all along this paper. While [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF] considered only the deterministic model (𝜎 = 0), we also consider later on the stochastic model (𝜎 ≠ 0). For 𝜎 = 0, the strange at-tractor and the limit cycles are shown in Fig. 3. They are intertwined in a complex way.

For this dynamics, we define a strong El-Niño event as any situation when 𝑥 becomes larger than the threshold 𝑥 𝑐 = -1. As can be seen from Fig. 3, this only happens in the limit cycle, and not in the strange attractor. Note that with this choice of parameters, the return period of strong El Niño events on the limit cycle is around 50 years (the time unit is 𝑡 * = 105 days and the period is 186 nondimensional time units), rather than 15 for the parameters studied by [START_REF] Roberts | Mixed-mode oscillations of El Niño-Southern Oscillation[END_REF].

For 𝜎 = 0, when the parameter 𝑎 is time periodic rather than constant, mimicking a seasonal forcing, [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF] observed transitions between the strange attractor and the limit cycle. In the following we consider a constant 𝑎, and rather study noise induced transitions between these two attractors, for 𝜎 ≠ 0.

We now discuss qualitatively the effect of the noise level 𝜎. For small 𝜎 > 0, the dynamics can switch from one regime close to the strange attractor to another regime close to the limit cycle. This is evident by looking at Fig. 4, where the time evolution of the variable 𝑥 for 𝜎 = 5 • 10 -5 is shown. As can be seen, the system switches between epochs where the variable 𝑥 is bounded around the value 𝑥 = -2.5 (strange attractor regime) and epochs where there are large oscillations with a period around 50 years (limit cycle regime). Such shifts can be interpreted as a toy model for hypothesized regime changes in the long-term ENSO dynamics. Strictly speaking, for 𝜎 ≠ 0, there is not anymore two attractors. However, as illustrated in Fig. 4, for small enough 𝜎, the trace of the two deterministic attractors is clearly visible. For small 𝜎, we will thus continue to discuss the strange attractor and the limit cycle, for simplicity. Strong El-Niño events occur only during The dynamics exhibits a switching between a regime close to the strange attractor (bounded regions around 𝑥 = -2.5) and a limit cycle regime.

the limit cycle regime. In the next section, we will study transitions from the strange attractor regime to the limit cycle regime (or equivalently to the strong El Niño regime).

The frequency of transitions from the strange attractor to the strong El-Niño regime increases as the amplitude of the noise increases. For large values of 𝜎, the dynamics is completely dominated by the noise and the distinction between the two attractors becomes meaningless.

Statistics of the first exit times for transitions to strong El Niño regimes

As discussed in the previous section, we define strong El Niño events as periods of time when 𝑥 𝑐 > -1, which occur along the limit cycle. In this section, we study transitions from the strange attractor regime to the strong El Niño regime, and how their statistics change when the noise amplitude 𝜎 is varied.

We consider X(t) = (𝑥(𝑡), 𝑦(𝑡), 𝑧(𝑡)) solutions to the stochastic Jin and Timmermann model (4). We define first exit times from a point x to the strong El Niño regime as

𝜏 𝑐 (x) = inf{𝑡 > 0 : 𝑥(𝑡) > 𝑥 𝑐 | X(0) = x}. (5)
The random variable 𝜏 𝑐 (x) depends both on the realization of the noise and on the initial condition x. The statistics are understood as averages over both the noise realization and the invariant measure of x over the strange attractor of the deterministic system (𝜎 = 0), the so called SRB measure.

For instance the mean first exit time E[𝜏 𝑐 ] is defined as

E[𝜏 𝑐 ] = ∫ dx 𝜌 SRB (x)E noise [𝜏 c (x)]. (6)
where E 𝑛𝑜𝑖𝑠𝑒 [•] is the expectation with respect to the noise realization and dx 𝜌 SRB (x) is the SRB measure. The SRB measure is defined through time averages of the deterministic dynamics (𝜎 = 0). In practice, we thus compute a very long trajectory of the deterministic dynamics. We then choose a set of 1000 initial conditions x taken randomly among all the points of this deterministic trajectory. Then, for any fixed value of 𝜎 > 0, for any initial condition x, we compute the first-passage time 𝜏 𝑐 for several noise realizations.

In function is close to an exponential: 𝑝(𝜏 𝑐 ) = 𝜆𝑒 -𝜆𝜏 𝑐 . The parameter 𝜆 is then equal to the inverse of the mean first exit time:

𝜆 -1 = E[𝜏 𝑐 ].
Because typically 𝜏 𝑐 (x) is much larger than the relaxation time to the strange attractor, one might expect that for most of the points of the strange attractor the dependence of 𝜏 𝑐 (x) on x is practically irrelevant. Indeed, we have verified numerically that except for a small region around the transition paths, the statistics are independent from the initial condition, up to numerical accuracy. Hence we have

E[𝜏 𝑐 ] E 𝑆𝑅𝐵 [𝜏 𝑐 (x)] E 𝑛𝑜𝑖𝑠𝑒 [𝜏 𝑐 (x)]
for generic points x close to the strange attractors.

As illustrated in Fig. 5, the mean first exit time E[𝜏 𝑐 ] is of the order of 1,000 in non-dimensional time units. The measured value is 𝜏 𝑐 = 1, 039. As 𝜏 𝑐 is much larger than the mixing time of the SRB measure, of order 1, then it is natural to expect that the first exit times should be random and distributed, with a very good approximation, according to a Poisson statistics. The observed exponential distribution is consistent with such a Poisson statistics. Similar exponential distributions for first exit times were observed for the deterministic dynamics with periodic modulation of the 𝑎 coefficient [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF].

We now study how the mean first exit time E[𝜏 𝑐 ] varies when the noise amplitude 𝜎 is changed. One generally expects an Arrhenius law:

E[𝜏 𝑐 ] ∼ 𝜎→0 𝐴𝑒 Δ𝑉 𝜎 2 . ( 7
)
Arrhenius laws where first derived by Kramers for gradient dynamics forced by white noise x = -d𝑉/dx+ √ 2𝜎𝜂(𝑡), where Δ𝑉 (the potential barrier) is the difference of potential between the original attractor and the saddle-point separating the basins of attraction of the two attractors (see for example the textbook by [START_REF] Gardiner | Handbook of stochastic methods[END_REF]). The Jin and Timmermann model is however not a gradient dynamics, and the function 𝑉 is not explicit. For such non gradient systems, the exponential factors of the Arrhenius law can be justified through a Laplace principle for a path integral representation of the transition probabilities, or asymptotic studies of Fokker-Planck operators [START_REF] Graham | Macroscopic potentials, bifurcations and noise in dissipative systems[END_REF], or through large deviation theory [START_REF] Freidlin | Random Perturbations of Dynamical Systems[END_REF]. The function 𝑉 is then called the quasipotential, which can be computed through a variational problem, or computing viscosity solutions of a Hamilton-Jacobi equation. The sub-exponential prefactor 𝐴 in Eq. ( 7) can be computed through Eyring-Kramers formulas, derived either for gradient [START_REF] Gayrard | Metastability in reversible diffusion processes i: Sharp asymptotics for capacities and exit times[END_REF] or non-gradient dynamics [START_REF] Bouchet | Generalisation of the Eyring-Kramers Transition Rate Formula to Irreversible Diffusion Processes[END_REF], for transitions from a point attractor and through a point saddle. Many generalizations exist, for instance for periodically modulated systems [START_REF] Dykman | Activated escape of periodically modulated systems[END_REF] or systems approaching a bifurcation [START_REF] Herbert | Predictability of escape for a stochastic saddle-node bifurcation: when rare events are typical[END_REF]. In large dimensional systems related to climate dynamics, effective Arrhenius laws have been observed numerically, for instance in transitions in beta-plane turbulence [START_REF] Bouchet | A rare event algorithm links transitions in turbulent flows with activated nucleations[END_REF] or in a simplified climate model with ice-albedo feedback [START_REF] Lucarini | Transitions across melancholia states in a climate model: reconciling the deterministic and stochastic points of view[END_REF].

When one studies a transition from a strange attractor rather than a point attractor, or when the saddle set between the two basins of attraction is a not a simple saddle point, then no theorem exists to put the Arrhenius law (7) on firm mathematical ground. However it has been argued for a long time [START_REF] Graham | Macroscopic potentials, bifurcations and noise in dissipative systems[END_REF], that if a finite distance 𝑑 > 0 exists between the strange attractor and the saddle set, then there is a non-zero quasipotential difference Δ𝑉 > 0, and an Arrhenius law should generically be expected.

In Figure 6: Mean first exit time E[𝜏 𝑐 ] for the transition from the strange attractor regime to the strong Niño event regime, as a function of the noise amplitude 𝜎, in loglog coordinates. In the limit of small noise amplitude 𝜎 < 10 -4 , E[𝜏 𝑐 ] seems to be closer to a power-law 𝜎 -1 (green line) than to the standard Arrhenius law. The inset shows the mean first exit time E[𝜏 𝑐 ] as a function of 𝜎 -2 for small values of 𝜎. This clearly shows that E[𝜏 𝑐 ] does not follow an Arrhenius law (represented by the green line). dimensional units for weak noise). Fig. 6 clearly shows that the mean exit time from the strange attractor to the regime of strong Niño events does not follow an Arrhenius law of the form (7). The mean exit time seems much closer to a power law E[𝜏 𝑐 ] ∝ 𝜎 -1 . An even clearer visualization of the violation of the Arrhenius law is provided in the inset of Fig. 6 where E[𝜏 𝑐 ] is shown as a function of 𝜎 -2 , for 𝜎 1.

In order to check the numerical robustness of our result, we computed E[𝜏 𝑐 ] using two different schemes of integration, and we checked numerical convergence with respect to the time step Δ𝑡 in the integration schemes. The first integration scheme is the fourth order Runge Kutta method to which a zero mean gaussian white noise is added. The variance of the noise is proportional to the integration time step Δ𝑡. In this way we have a precision of Δ𝑡 4 for the deterministic part while we make an error of order

√

Δ𝑡 for the statistics. The second integration scheme is the stochastic Runge Kutta method which has an error of order Δ𝑡 for a stochastic dynamics [START_REF] Roberts | Modify the improved euler scheme to integrate stochastic differential equations[END_REF].

Let us note that many other behaviors than exponential ones have been observed for mean exit times. For instance a power-law has been observed for flow reversals in numerical simulations of inviscid turbulent flows [START_REF] Shukla | Statistical theory of reversals in two-dimensional confined turbulent flows[END_REF]. However for this last example, as the dynamics is not a deterministic system with attractors perturbed by weak noise, it was not clear why one should have expected an Arrhenius law in the first place.

We observe a breakdown of the Arrhenius law for the Jin and Timmermann model which is a deterministic system with attractors perturbed by weak noise. This is striking. Indeed we stress again that if a finite distance 𝑑 > 0 exists between the strange attractor and the saddle set, then there is a non-zero quasipotential difference Δ𝑉 > 0, and an Arrhenius law should be expected. The distance 𝑑 might be expected to be generically strictly larger than 0.

We see two possible heuristic explanations for this interesting breakdown. The first explanation might be that a finite distance 𝑑 > 0 and a quasipotential barrier Δ𝑉 > 0 between the strange attractors and the basin of attraction limit cycle do actually exist, but they are extremely small. Then the explanation of the observed breakdown in Fig. 6 would be that we have not studied small enough values of 𝜎. We note however that we computed first exit times of order E[𝜏 𝑐 ] = 5.10 5 for values of 𝜎 as small as 10 -6 . If this first explanation is valid, this means that the Arrhenius law is practically irrelevant even if it might be mathematically correct.

The second possible explanation might be that there exists no finite distance between the strange attractor and a possible fractal boundary between the basins of attractions. Then for any small values 𝑑 and 𝑣, there always exist points in the strange attractor and in the boundary of the basin of attraction at a distance smaller than 𝑑 and a quasipotential differences Δ𝑉 smaller than 𝑣. Many phenomenologies could then be imagined, for instance with a distribution of a large number of transition paths, possibly infinite, leading to a power law or an effective behavior of the first exit times described by any function. Those conjectures are not based on any mathematical results yet. However the possibility of a breakdown of the Arrhenius law is a very interesting problem, that should be studied further either through theory and mathematics, or through numerical simulations.

Committor function of the Jin and Timmermann model

In Sec. 3, we have shown that, in the stochastic Jin and Timmermann model, transitions between the strange attractor regime and the strong El-Niño regime occur at random times, in the limit of small noise 𝜎 → 0. In this section, we focus on the associated prediction problem: What is the probability that a strong El Niño event occurs within a given timeframe, given the state of the system at the time of prediction? We will address this question for any finite value of the noise amplitude 𝜎.

We consider solutions X(𝑡) = (𝑥(𝑡), 𝑦(𝑡), 𝑧(𝑡)) of the stochastic Jin an Timmermann model (4). We remind the reader that we identify a strong El Niño with an event when 𝑥 > -1. For a solution X(𝑡) that starts from x, that is X(0) = x, we want to predict the probability 𝑞(x) that a strong El Niño event occurs within a fixed time 𝑇. This is

𝑞(x) = P max 0≤𝑡 ≤𝑇 𝑥(𝑡) > -1 | X(0) = x . (8)
Recalling the definition of the first passage time to a strong El Niño regime, Eq. ( 5),

𝜏 𝑐 (x) = inf{𝑡 > 0 : 𝑥(𝑡) > -1 | X(0) = x}, (9)
we note that 𝑞(x) = P[𝜏 𝑐 (x) < 𝑇] is the cumulative distribution function (CDF) of the first-passage time. We now define committor functions and explain that 𝑞 is a committor function.

Committor functions. For a Markov stochastic process {Y(𝑡)} which takes values in Γ, we define the first hitting time of the set 𝐶 as 𝜏 𝐶 (y) = inf{𝑡 : Y(𝑡) ∈ 𝐶 | Y(0) = y}. For two disjoint subsets 𝐴, 𝐵 ⊂ Γ, the committor function q(y) is defined as the probability to hit the set 𝐵 before hitting the set 𝐴: q(y) = P(𝜏 𝐵 (y) < 𝜏 𝐴 (y)).

(10)

Considering the auxiliary process {Y(𝑡)}, with Y(𝑡) = (X(𝑡), 𝑡), and the two sets

𝐴 = {y = (x, 𝑡) | 𝑥 > -1 and 𝑡 ∈ [0,𝑇]} and 𝐵 = {y = (x,𝑇) ; 𝑥 ≤ -1}, (11)
we see that 𝑞(x) = q(x, 0). Hence 𝑞, in Eq. ( 8) is a committor function.

For an ergodic process, replacing statistical averages by temporal averages in (10), and using y = (x, 𝑡), we have

𝜌(x)𝑞(x) = lim 𝑡→∞ 1 𝑡 ∫ 𝑡 0 d𝑡 𝛿 (X 𝑡 -x) 1 {𝜏 𝐵 ≤𝜏 𝐴 } , and 𝜌(x) = lim t→∞ 1 t ∫ t 0 dt 𝛿 (X t -x) , (12)
where 𝜌(x) is the stationary distribution function of 𝑋, 𝛿 is a Dirac delta function, and 1 {𝜏 𝐵 ≤𝜏 𝐴 } takes value 1 if 𝜏 𝐵 ≤ 𝜏 𝐴 and 0 otherwise. The formulas ( 12) can be used to estimate 𝑞(x) from an observed trajectory {X(𝑡)} of the dynamical system. For the sake of completeness, it should be said that when the dynamics is a stochastic differential equation, the committor function 𝑞(x) is the solution of the Dirichlet problem [START_REF] Ren | Transition pathways in complex systems: Reaction coordinates, isocommittor surfaces, and transition tubes[END_REF][START_REF] Thiede | Galerkin approximation of dynamical quantities using trajectory data[END_REF].

To illustrate the concept of predictability margin introduced in Sec. 1, we choose the value 𝑇=200 in nondimensional time units, which is slightly larger than the period of the limit cycle (the "natural" periodicity of strong El-Niño events which is 186), and of the order of the Lyapunov time. This choice guarantees that for the deterministic dynamics, 𝜎 = 0, each trajectory starting in one point of the limit cycle almost certainly will reach the threshold 𝑥 𝑐 = -1. Note that for 𝜎 = 0 the same behavior would occur for any 𝑇 larger than the limit cycle period. For 𝑇 smaller than this period, there would be a portion of the cycle that would not reach the threshold in the given time. Instead, for 𝜎 ≠ 0, changing 𝑇 has the same effect on the results as changing 𝜎 as it allows the noise more or less time to cause a transition.

a. Description of the committor function: deterministic and probabilistic predictability

Figure 7 shows the committor function 𝑞(x), for different values of 𝜎. As 𝑞 is a function of 3 variables (𝑥, 𝑦, 𝑧), we have chosen to represent cuts of 𝑞 in different planes. We will discuss in detail the cut of 𝑞 along the plane 𝑥 = -2.831 (Fig. 7) and also cuts along the planes 𝑦 = -1.1580 and 𝑧 = 1.3409 (Fig. 8 and Fig. 9, respectively). To compute the committor function 𝑞(x) on the different planes we adopted the following strategy:

1. Discretize the plane into 𝐾 = 𝐿 × 𝐿 cells C 𝑘 .

2. For each cell C 𝑘 , generate 𝑁 = 1000 trajectories starting from a point x ∈ C 𝑘 .

3. Count the number of trajectories 𝑁 1 that reach the threshold 𝑥 𝑐 before time 𝑇.

4. Estimate 𝑞(x) for x ∈ C 𝑘 as q(x) = 𝑁 1 𝑁 . This method is less efficient, from a computational point of view, than the one based on Eq. ( 12). In fact, in the former we only use the information carried by the initial condition while in the latter we use the information carried by a much more significant part of the trajectory. The committor function computed from long trajectories using Eq. ( 12) will be discussed in Sec. 4 c.

1) C

F -(𝜎 = 0)

In the deterministic case (𝜎 = 0, Fig. 7a), as the future is completely determined by the initial condition, 𝑞 can only take values 0 or 1. On Fig. 7a, we can distinguish three regions corresponding to two very different situations. First, two regions correspond to uniform values of the committor function: in the yellow area 𝑞 = 1, when trajectories reach the threshold within a time 𝑇, corresponding to large values of 𝑧; in a thick purple band 𝑞 = 0, when no trajectory reaches the threshold. In those two regions, the occurrence of strong El-Niño events is easily predicted. Everywhere else, we see very fine filaments of alternating yellow and purple values. In this third region, because of the sensitive dependence on the initial conditions, a small but finite initial perturbation, of the order of 1% of the values of 𝑥 or 𝑦, leads to a different outcome. In this region, the occurrence of strong El-Niño events is very difficult to predict.

We stress that a precise definition of this third area is not intrinsic, it depends on the actual precision with which the values of 𝑥 and 𝑦 can be measured. However the distinction between areas with easy predictability and areas with difficult predictability, might be crucial at a practical level. This important practical point will be discussed more in section 4.d.

One might ask what relationship exists between the regions outlined above and the basins of attraction of the system. However, this relationship is less trivial than one might expect. Although some regions reflect the structure of the basins of attraction, this is not true in general. In fact, there are points in the basin of attraction of the strange attractor which pass the threshold before reaching the strange attractor, as well as points in the basin of attraction of the limit cycle which do not reach the threshold within the time T.

2) C

F 𝜎 ≠ 0 Figures 7b and7c show the committor function in the case where a finite noise amplitude 𝜎 ≠ 0 is considered. As can be seen by comparing Figs. 7a, 7b and 7c, adding a small noise blurs the visible structures of the deterministic case. For larger noise values (𝜎 = 10 -3 ), Fig. 7c shows that the committor function looks smooth nearly everywhere (mathematically it is smooth everywhere, smooth here is used qualitatively and means with mild variations). This means that the deterministic predictability is lost for most initial conditions as (0 < 𝑞 < 1). Then one cannot expect to predict the outcome in the way of a deterministic forecast. However, the occurrence of strong El-Niño events is probabilistically predictable: the value of the probability can be determined in practice with an excellent precision as it changes very slowly when one changes the initial conditions. It can also be seen on the figure that the occurrence of strong El-nino events is frequent (𝑞 > 0.6 almost everywhere). This is an indication that for such a value of 𝜎 we are in the noise-dominated regime.

The most interesting case is probably the one with the intermediate noise amplitude 𝜎 = 0.00005. In Fig. 7b, we delineate 4 regions: first, two regions of perfect (deterministic) predictability, where the event will occur with probability very close to 1 (region 1) or to 0 (region 2). Second, there exists a probabilistically predictable region (region 3) with good (probabilistic) predictability properties, where a value 0 < 𝑞 < 1 can clearly be predicted with very mild dependence with respect to the initial conditions. Finally, the region 4, which is unpredictable in practice. In this region, the strong dependence with respect to the initial condition prevents any practical prediction, either deterministic or probabilistic, of the precise value of 𝑞. While regions 1, 2 and 4 are reminiscent of their deterministic counterparts (Fig. 7a), region 3 is not. Instead, the behavior in this region is similar to the strong noise case shown in Fig. 7c. It is a region where the stochasticity is large enough to smooth out the deterministic values of 𝑞. The fact that it occurs even at very low noise amplitude is probably related to extremely unstable parts of the phase space, for instance for trajectories passing close to unstable fixed points or orbits. We study the dynamical mechanisms which lead to this area with good probabilistic predictability in section 4.b.

We stress that the existence of these 4 types of predictability regions, and especially the new and most interesting probabilistically predictable region (region 3), should be generic for most prediction problems in climate dynamics. In the present section, we have assumed a perfect model and implicitly assumed an observation scale for determining the initial conditions. In section 4.d, we will discuss the effect of a changing observation scale, and of possible model errors, on the meaning of these 4 types of predictability regions.

b. Dynamical characterization of the probabilistically predictable region

In order to understand the reason for which the probabilistically predictable region arises, it can be useful to introduce a quantity which characterizes sensitivity to small perturbations. This quantity is the largest finite time Lyapunov exponent and it is a measure of the sensitivity to initial conditions for the deterministic system. Let us consider two trajectories of the deterministic system (𝜎 = 0) with initial conditions x and x + 𝛿x and let Δ(𝑡) be the value of the euclidean distance between the two trajectories at time 𝑡. The largest Lyapunov exponent Λ 𝐿 is defined as

Λ 𝐿 = lim 𝑡→+∞ lim Δ(0)→0 1 𝑡 log Δ(𝑡) Δ(0) . (13)
Since we are dealing with predictions with a time horizon 𝑇, we believe it is more appropriate to define a finitetime version of Λ 𝐿 . Hence, we compute the largest finite time Lyapunov exponent 𝜆 𝐿 as 𝜆 𝐿 = 1 𝑇 log ( Δ(𝑇 ) Δ( 0) ). Note that the initial perturbations have to be considered small but finite as we have taken x and x + 𝛿x into the same cell C 𝑘 . Positive values of 𝜆 𝐿 mean that the distance between the trajectories grows exponentially. This quantity is shown for the Jin-Timmermann model in the 𝑥 = -2.831 plane in Fig. 10. Comparing this figure with Fig. 7a is enlightening: the regions of perfect predictability (1 and 2) are associated to values of 𝜆 𝐿 which are either negative or close to 0. Moreover, 𝜆 𝐿 < 0 for almost all regions for which the deterministic committor is equal to 1 (see the two thick yellow bands at the bottom left of Fig. 7a). From this description, it seems that 𝜆 𝐿 is not the correct quantity to explain the emergence of the probabilistically predictable region. However, from a careful analysis of Fig. 10 it can be noted that there is a region, close to the left boundary of region 2, for which the values of 𝜆 𝐿 are positive also for points x such that 𝑞(x) = 1. It means that this region is quite unstable with respect to small initial perturbations. Since this region belongs to the probabilistically predictable region, it is reasonable to say that the appearance of region 3 is related with this instability. In fact, although the instability region is a subset of the probabilistically predictable region, it should be noted that for 𝜎 ≠ 0 the system is perturbed at any time. The ensemble of these small perturbations gives rise to a finite perturbation which could explain the growth of the region of instability.

To reinforce this conjecture, we compute the averaged value of the euclidean distance between two different trajectories x 1 (𝑡) and x 2 (𝑡), with the same initial conditions but different realisations of the noise: 𝑑 max = max 𝑡 ∈ [0,𝑇 ] x 1 (𝑡)x 2 (𝑡) 2 . Figure 11 shows 𝑑 max as a function of the initial condition in the plane 𝑥 = -2.831 for 𝜎 = 0.00005: it can be seen that the typical values of 𝑑 max are larger in the probabilistically predictable region than in other regions. This means that trajectories starting space. Therefore, the emergence of the probabilistically predictable region is qualitatively associated to an instabil-ity present in the deterministic system which is accentuated as the amplitude of the perturbations increases.

c. Committor function computed from long trajectories

In this section we discuss the committor function computed from an ensemble of long trajectories by means of Eq. ( 12). The motivation is that the very precise strategy adopted in Sec. 4a is unlikely to be adapted for real-world problems as well, because it requires a dynamical model and has a very large computational cost. Indeed, it requires to sample an ensemble of trajectories for every point in the phase space. For a data-based approach, it is usually possible to observe only the evolution of a trajectory (or of an ensemble of trajectories) over a very long time. However, individual trajectories, regardless of their length, do not fill the whole phase space. Indeed, they usually concentrate on the region where the invariant measure of the stochastic system is concentrated. Hence, the strategy adopted in Sec. 4a is appropriate for computing the committor function in an arbitrary plane while the use of Eq. ( 12) allows the computation of the committor function in the region where the invariant measure is concentrated.

The committor function computed using Eq. ( 12) is shown in Fig. 12, for the same three values of the noise amplitude as in Sec. 4a. For this figure, we used 10 000 trajectories of length T = 10 4 (nondimensional time units) initialized in the strange attractor. Note that T is 10 times bigger than the first exit time 𝐸 [𝜏 𝑐 ] for 𝜎 = 0.00005 and 100 times greater than 𝐸 [𝜏 𝑐 ] for 𝜎 = 0.001 (see Fig. 6). This guarantees that the trajectories will be distributed according to the invariant measure of the system.

As already mentioned, Fig. 12 shows that the trajectories do not cover all the phase space but they are concentrated in a certain region. Furthermore, we can see that it is more appropriate to call it a manifold rather than a region. In fact, if it were an object of dimension 3, its intersection with a plane should define an area on that plane. Instead, it appears that the intersections between the object and planes are lines rather than areas. This is illustrated in Fig. 13 which shows an intersection between the object in Fig. 12b and the plane 𝑥 = -2.831 (the same plane as in Fig. 7). This leads us to conclude that the trajectories are distributed over a manifold of dimension smaller than 3.

The comparison between Fig. 13 and Fig. 7b allows us to make another important remark on the committor function on the manifold. As can be seen in Fig. 13, the committor function takes values between 0 and 1 at the two ends of the line. For the rest, the line is made up of segments on which 𝑞 takes the values 0 or 1. It is straightforward to recognize that the two ends of the line belong to the probabilistically predictable region while the central part of the line belongs to region 4. This highlights that the committor function computed from long trajectories provides useful information for many of the states of the system, that is, for all typical states. However, for atypical conditions where we have little or no information, using 𝑞 to make predictions can lead to erroneous results.

Having made these necessary considerations, we can continue the description of the committor computed from long trajectories. By comparing Fig. 12 with Fig. 3 we can immediately identify the strange attractor in the manifold on which the committor function is represented. The limit cycle is more difficult to visualize but its presence can be deduced from the spiral behavior present in the region inside the strange attractor and from the shape of the manifold boundaries that follow the shape of the limit cycle in Fig. 3.

The qualitative structure of the committor function on the manifold is similar to the one observed on plane cuts: in the deterministic case, we observe regions of perfect predictability, and regions where the sensitivity to initial conditions make it unpredictable in practice. When the noise is sufficiently strong, regions of probabilistic predictability appear, where a finite value of the probability of a strong El Niño event 0 < 𝑞 < 1 depends only mildly on the initial conditions (Fig. 12c). It can be noted, however, that the intermediate case, analogous to Fig. 7b, with coexistence of a region of deterministic predictability, a region of probabilistic predictability, and an unpredictable region, is more difficult to observe in this visualization.

Finally, we underline that the region of unpredictability, made by thin filaments where the committor is a highly fluctuating function, emphasizes again how the two attractors are intertwined in a complex way.

d. The effect of the observation scale and of the model resolution on the predictability areas

In section section 4.a.1 and in section 4.a.2, we have described 4 regions depending on their predictability properties. In Fig. 7b, the distinction between the probabilistically predictable region (region 3) with good (probabilistic) predictability properties, and region 4 which is unpredictable in practice, is very clear visually. However we stress again that this distinction, of high practical relevance, is not intrinsic and that the probabilistic predictability properties can be characteristic of a system or emerge due to ignorance of the initial conditions or the lack of a sufficiently precise model.

Indeed, for practical applications, the initial condition is known only up to some errors. Let us for instance assume that the initial condition is known up to a given error that is characterized by a scale 𝑙 𝑜 . For this case, the practitioner would see an effective committor which would be the real committor function (for instance as represented in Fig. 7b)) which would be spatially averaged at the observation scale 𝑙 𝑜 . If 𝑙 𝑜 is larger than the typical size of the filaments in region 4, then the unpredictable region will appear practically as having a smooth committor function at the scale 𝑙 𝑜 , and the practitioner will not see qualitative differences between the regions 3 and 4. However, if the practitioner is able to improve her knowledge of the initial condition, and decrease 𝑙 𝑜 below the size of the filaments, she will learn that the value of the committor function varies abruptly with the initial condition in the unpredictable region 4, while it will maintain smooth and more or less constant values when changing 𝑙 𝑜 in the probabilistically predictable region (region 3). The prediction in the probabilistically predictable region remains robust, although probabilistic.

In the case of an imperfect model, for instance due to a finite resolution 𝑙 𝑟 , the amplitude of the effective noise should decrease when 𝑙 𝑟 decrease. This idea of somehow too simplistic, as the deterministic term itself should be slightly modified. However it helps understand qualitatively the effect of the model resolution on the committor function. If the practitioner is able to improve his model resolution by decreasing 𝑙 𝑟 , he will be able to progressively decrease the noise amplitude and uncover the fine structure of the committor function in the unpredictable region 4, and he will observe abrupt variations of the committor function when the initial condition is changed. By contrast, in the region of probabilistically predictable region, the change of the resolution will not lead to qualitative or important quantitative changes for the committor function, the prediction will be robust even with coarser resolution.

Conclusion

In this paper, we have introduced a mathematical concept, the committor function, encoding the probability that an event occurs within a given time, conditioned on the state of the system at the time of prediction. We believe it is an appropriate concept for many prediction problems in climate science in a range of time scales which we call the predictability margin. It corresponds to timescales for which a deterministic description of the system is no longer relevant, because of the sensitive dependence to initial conditions, but for which more precise probabilistic predictions than the climatological one can be made, because the system has not yet forgotten completely the initial condition.

In the context of a simple, low-dimensional stochastic model, the Jin-Timmermann model, in a regime of coexistence of a limit cycle and a strange attractor found by [START_REF] Guckenheimer | Un) predictability of strong El Niño events[END_REF], we have shown that noise could induce transitions between the two attractors. These transitions correspond to regime shifts regarding the occurrence of strong El Niño events, which are periodic in the limit cycle, with a return time close to 15 years, and which do not occur at all in the strange attractor (in the deterministic case). In the stochastic case, the occurrence of strong El Niño events therefore becomes random, and the waiting times follow a Poisson statistics.

In this example, we have shown that the probability of occurrence of strong El Niño events had different predictability properties depending on the state of the system at the time of prediction. The most important result is that there exist regions of probabilistic predictability, where the event has a finite probability of occurring 0 < 𝑞 < 1, and this probability does not depend sensitively on the initial state, and regions of probabilistic unpredictability where the probability changes a lot if one changes by a small and finite amount the initial condition. We expect the existence of this dichotomy between probabilistically predictable and probabilistically unpredictable regions to be a generic feature for climate prediction problems at the predictability margin. We stress that this notion depends on the precision with which the initial condition can be assessed and on the model resolution. We also discussed that in the regions of probabilistic predictability, the prediction is expected to be more robust that in the regions of probabilistic unpredictability, when either the initial conditions, the observation scale, or the model resolution scale, are changed.

We have also discussed the methodological aspects for computing the committor functions. For our example, a small stochastic perturbation of a chaotic deterministic system, we have computed the committor function using two approaches. First, by direct sampling of ensembles of initial conditions close to any point in phase space, and second, through a data based approach using observed trajectories. As soon as the number of degrees of freedom increases, the first method will become impossible to use in practice, because of the numerical cost. The second method may sometimes be associated with sampling issues, as one can get meaningful results only for the parts of the phase space that have been visited many times. Another method mentioned in Sec. 4, would be to solve a backward Kolmogorov equation by employing classical finite difference methods. This method is impractical for systems with more than a few degrees of freedom. To be able to sample efficiently committor functions in large dimensions, more efficient data-based methods will be necessary, relying either on classical statistical methods or machine learning methods [START_REF] Lucente | Machine learning of committor functions for predicting high impact climate events[END_REF]. Furthermore, such data-based methods are also effective for solving backward Kolmogorov equations [START_REF] Thiede | Galerkin approximation of dynamical quantities using trajectory data[END_REF][START_REF] Finkel | Learning forecasts of rare stratospheric transitions from short simulations[END_REF][START_REF] Lucente | Coupling rare event algorithms with data-based learned committor functions using the analogue markov chain[END_REF].

A natural question is to understand the difficulty to compute committor functions using the best up to date atmosphere or climate models. In a series of forthcoming papers, we develop a method based on machine learning for computing committor functions for the prediction of midlatitude extreme heat waves using data from one of the best current climate model, CESM2.1, the model used by NCAR for the CMIP experiments.

Figure 1 :

 1 Figure1: Schematic illustrating the concept of predictability margin: deterministic predictability is only possible until a finite time (e.g. the Lyapunov time). The associated mathematical problem is an initial value problem (IVP). Long term statistical properties (beyond the mixing time) do not depend on the initial condition, and the corresponding mathematical object is the invariant measure.In the intermediate range of timescales, which we call the predictability margin here, the appropriate mathematical concept is the committor function, which encodes the probability of a given event to occur, condition on the state of the climate system at the time of the prediction.

  Figure 2-b illustrates such a periodic orbit, with the parameter values [𝛿, 𝜌, 𝑐, 𝑘, 𝑎] = [0.2625, 0.3224, 2.3952, 0.4032, 6.8927] and dimensional normalization constants [𝑇 0 , 𝑡 * , ℎ * ] = [2.8182 • C, 104.9819 days, 62 m]. Roberts et al. (

Figure 2 :

 2 Figure 2: Top plot: observed sea surface temperature anomalies over the last decades, spatially averaged over the Niño-3 region. Bottom: eastern Pacific sea surface temperature anomalies simulated with the deterministic Jin and Timmermann model. The red segments in the upper panel highlight the growing sequence of oscillations leading to a strong El-Niño event, which are qualitatively reflected in the behavior of the model shown in the lower panel.

Figure 3 :

 3 Figure 3: The two intertwined attractors of the Jin-Timmermann model (limit cycle in blue and strange attractor in red).

Figure 4 :

 4 Figure 4: Time series of the variable 𝑥 for 𝜎 = 0.00005.The dynamics exhibits a switching between a regime close to the strange attractor (bounded regions around 𝑥 = -2.5) and a limit cycle regime.

Figure 5 :

 5 Figure 5: Logarithm of the Probability Density Function of the first exit time between the strange attractor and the limit cycle for 𝜎 = 5 × 10 -5 sampled by direct integration (red) of the stochastic differential equation (4), and the exponential distribution 𝑝(𝜏 𝑐 ) = 𝜆𝑒 -𝜆𝜏 𝑐 with 𝜆 -1 = E[𝜏 𝑐 ] (green).

  Fig. 6, we show the mean first exit time E[𝜏 𝑐 ] as a function of the noise amplitude 𝜎. It ranges from transition times of about 25 years (𝑇 = 100 in non-dimensional time units, for the strongest values of the noise amplitude) to several millenia (about 3000 years for 𝑇 = 10000 in non-

Figure 7 :

 7 Figure 7: Color plot of the committor function 𝑞(𝑥, 𝑦, 𝑧) in the plane 𝑥 = -2.8310, for three values of the noise amplitude, 𝜎 = 0 (left, deterministic), 0.00005 (middle) and 0.001 (right). Regions with uniform 𝑞 = 0 or 1 values correspond to deterministic predictability, smooth regions with 0 < 𝑞 < 1 to probabilistic predictability, and regions with sensitive dependence on initial conditions to unpredictable parts of phase space.

Figure 8 :

 8 Figure 8: Color plot of the committor function 𝑞(𝑥, 𝑦, 𝑧) in the plane 𝑦 = -1.1580, for 𝜎 = 0.00005 (left) and 0.001 (right).

Figure 9 :

 9 Figure 9: Color plot of the committor function 𝑞(𝑥, 𝑦, 𝑧) in the plane 𝑧 = 1.3409, for 𝜎 = 0.00005 (left) 0.001 (right).

Figure 10 :

 10 Figure 10: Color plot of the maximum finite time Lyapunov exponent 𝜆 𝐿 as a function of the initial condition (𝑥, 𝑦, 𝑧) in the plane 𝑥 = -2.831.

Figure 12 :

 12 Figure 12: Color map of the committor function 𝑞(𝑥, 𝑦, 𝑧) in correspondence to the most visited region of phase space, for 𝜎 = 0 (left, deterministic), 0.00005 (middle) and 0.001 (right).

Figure 13 :

 13 Figure 13: Color plot of the committor function for 𝜎 = 0.00005 at the intersection of the manifold with the plane 𝑥 = -2.831.

Acknowledgments. The work of D. Lucente was funded through the ACADEMICS grant of IDEXLYON, project of the Université de Lyon, PIA operated by ANR-16-IDEX-0005. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 753021 and the ANR grant SAMPRACE, project ANR-20-CE01-0008-01. Computer time was provided by the "Pôle Scientifique de Modélisation Numérique" in Lyon. We acknowledge useful scientific discussions with P. Abry, P. Borgnat, C.-E. Bréhier, J. Finkel, T. Lelièvre, J. Rolland and J. Weare.

Data availability statement. This theoretical work does not use external datasets. The datasets on which this paper is based as well as the numerical code are not archived on a public data repository. All the information for reproducing the results are provided in the main body. Additionally, the corresponding author may provide the data upon reasonable request.