The evolution of a non-autonomous chaotic system under non-periodic forcing: a climate change example

Complex Earth System Models are widely utilised to make conditional statements about the future climate under some assumptions about changes in future atmospheric greenhouse gas concentrations; these statements are often referred to as climate projections. The models themselves are high-dimensional nonlinear systems and it is common to discuss their behaviour in terms of attractors and low-dimensional nonlinear systems such as the canonical Lorenz `63 system. In a non-autonomous situation, for instance due to anthropogenic climate change, the relevant object is sometimes considered to be the pullback or snapshot attractor. The pullback attractor, however, is a collection of {\em all} plausible states of the system at a given time and therefore does not take into consideration our knowledge of the current state of the Earth System when making climate projections, and are therefore not very informative regarding annual to multi-decadal climate projections. In this article, we approach the problem of measuring and interpreting the mid-term climate of a model by using a low-dimensional, climate-like, nonlinear system with three timescales of variability, and non-periodic forcing. We introduce the concept of an {\em evolution set} which is dependent on the starting state of the system, and explore its links to different types of initial condition uncertainty and the rate of external forcing. We define the {\em convergence time} as the time that it takes for the distribution of one of the dependent variables to lose memory of its initial conditions. We suspect a connection between convergence times and the classical concept of mixing times but the precise nature of this connection needs to be explored. These results have implications for the design of influential climate and Earth System Model ensembles, and raise a number of issues of mathematical interest.


Introduction
The theory of non-autonomous nonlinear dynamical systems has enjoyed great popularity over the past few decades, particularly within the climate modelling community [1].This is because complex global climate models, or rather Earth System Models (ESMs), which are widely used to make projections of the 21st century and to support the IPCC's climate assessment reports, are subject to non-periodic, climate-change-like forcing, which inevitably breaks their autonomy.These models are also high-dimensional, multi-component, multi-scale, chaotic nonlinear systems and as a consequence, any forward computation -that is to say, projection of the future within the model -is highly sensitive to the finest details † of the initial state, making climate prediction a non-trivial task.
Uncertainty in the state from which to initialise ESMs is known as initial condition uncertainty (ICU).The sensitivity of such climate system models to ICU is well known since the early 60s' [2] and has led to the development of ensemble weather forecasting [3].Its relevance for climate forecasting is also increasingly being recognised [4,5,6,7], as it is the necessity of using large climate initial condition ensembles (ICEs) to characterise ICU [8].Nevertheless, it is often assumed that the uncertainty arising from ICU can be addressed by taking statistics from a single, long trajectory, which it is assumed, over time, would explore all possible states in phase space.In a stationary system ‡ this is essentially an ergodic [9], or "kairodic" [8], assumption: that averages and distributions of states over long periods (e.g. 30 years for IPCC) are representative of any particular instant -with the caveat that it would require infinite time for convergence.However, under non-periodic, climate-change-like forcing -such as increasing atmospheric greenhouse gas concentrations -the system is not ergodic, and hence cannot be studied in this way [9].
The non-autonomous nature of ESMs under anthropogenic, non-periodic climate change forcing means that, in general, such a system do not possess an attractor in the classical sense, because we cannot take the asymptotic limit as time tends to infinity.Recent years have seen the emergence of a number of approaches from the mathematical community to address this issue [1,10].Central to these approaches is the idea that a model's climate can be formally seen as an evolving probability distribution constructed from an ensemble of simulations which have been initialised from different ICs, initialised in the very remote past.This can be thought of as multiple "evolutions" of the same Earth System (that is to say, they all obey the same physical laws) but with each one starting from different initial points [10].
For a wide class of nonautonomous systems, it has been shown that, in this "parallel climate realisations" approach, the correct concept to describe a time-dependent set in the phase space as the "limit" of a set of ICs is the pullback attractor.[11,12,13,14,15].Many climate models (including the one discussed here § ) satisfy some form of energy balance which typically implies the core structural hypotheses required to establish the existence of pullback attractors.At any instant in time, the system's 'climate' can therefore be taken as an instantaneous slice of the pullback attractor -this slice is the so-called snapshot attractor.Furthermore, in the same way that the (pullback) attractors are some form of "limit" for a set of IC's, the initial distribution of IC's might converge to a time dependent "pullback" probability measure supported on the pullback attractor.Invariant and pullback measures are typically not unique but here we are specifically interested in so-called natural or physical pullback measures, which emerge as the limit of smooth IC distributions ¶ [16].However, while mathematically appealing, these concepts are of limited use in supporting the construction of climate change ensembles of ESMs, and therefore in making climate projections and ultimately supporting society.By definition, the pullback attractor depends on initialisation infinitely far in the past ‖ .Generally, this problem can be overcome by noting that in most cases we can assume that mixing happens on finite time scales, which, however long, can be taken as providing a convergence time: the time taken for the ensemble dynamics to forget its initial state.We do not therefore require infinitely long simulations, only sufficiently long, where "sufficient" is defined by this convergence time.Nevertheless this means that the pullback attractor is only applicable for long term climate analyses -longer than the convergence time.This convergence time can be small (around 5 years) for a simple conceptual low-dimensional atmospheric model system [17] but rather long (over 150 years) even for fast-mixing atmosphere variables in an intermediate-complexity ESM [18].In other words, the pullback attractor approach might give us a good description of our idealised model system's climate by the end of the next century (i.e., in about 150 years time), but it can not tell us how we will get there.
This means that, while the pullback attractor represents the internal variability of the mathematical system on timescales beyond the convergence time, it is not the relevant object to represent climate on shorter timescales because it does not reflect knowledge regarding the current state of the climate system.On shorter timescales, the representative distribution is more constrained.The set of trajectories that make up this constrained distribution is a subset of those making up the pullback attractor, but it is not clear how the two distributions relate to each other.
Here we consider how to quantify this initial response and how such forward distributions can depend on both our knowledge of the current state and the characteristics of the non-autonomous forcing.These issues are critical to understanding what is required to make climate projectionseven in the perfect model scenario [19] -and in characterising the behaviour of non-autonomous, non-periodic, nonlinear systems more broadly.To do so, we use a low-dimensional system with characteristics of an ESM [20].The concept of an evolution set is introduced to describe the set on which a more constrained distribution would be supported.We also introduce the concept of an evolution distribution to describe the more constrained distribution and we consider the convergence time for this evolution distribution to become indistinguishable from the pullback invariant distribution.
The paper is divided as follows.In Section 2, we describe the model used in this study, as well as the experiments performed.In Section 3, we elaborate on the concept of the pullback attractor, demonstrate it with examples from our model, and define and illustrate the convergence time for different variables in a stationary situation.In Section 4, we approach the transient climate change problem in combination with some hypothetical, highly-constrained knowledge of the initial state -so-called micro ICUs [4,7].In Section 5 we consider situations where the initial state is not well constrained -so called macro ICUs [4,7], while revisiting the concept of convergence time in the non-autonomous situation.In Section 6, we explore the influence of the forcing on the evolving distributions.We then conclude the paper with Section 7, where we discuss further questions and future directions for the this study.
the importance of this concept is demonstrated, showing that such a proof would be worthwhile.

Model
We use a low-dimensional coupled ocean-atmosphere model, which is taken as a conceptual representation of a climate model.In this model, the ocean domain is presented as two connected but distinct basins, say, one representing the ocean at high latitudes and another representing it at low latitudes in the same hemisphere, with its dynamics given by the Stommel '61 (S61) model [21].The S61 model is based on the free convection controlled by density differences maintained by heat and salt exchange between the reservoirs.The atmosphere is represented by a simplified description of its large scale circulation in one hemisphere, given by the Lorenz '84 (L84) model [22,23].The L84 model is based on the interaction of the westerly, mid-latitude wind current and large scale, pole-ward eddies.The L84 model and the S61 model form the coupled ocean-atmosphere model used in this study, which we shall refer to as Lorenz 84-Stommel 61 (L84-S61) model.
Mathematically, the L84-S61 model consists on the following five coupled ODEs 2) where and (2.8) In the above, the variables X, Y, Z represent the high-frequency, atmospheric variables from the L84 model: X represents the intensity of the symmetric westerly wind, Y and Z are the Fourier amplitudes characterising a chain of large-scale eddies, which transport heat towards the pole at a rate proportional to their amplitude.The variables T, S are the slow ocean variables as in the S61 model: T and S denote the pole-equator temperature and salinity differences, respectively.The function f (T, S) represents the strength if the thermohaline circulation (THC), while F 0 (t) is the forcing due to seasonal variation in the heating contrast between the pole and equator.The latter corresponds to an average forcing equals F m which varies seasonally according to a cosine function with amplitude M , and can be forced towards another value at a rate H.All the variables in the model are non-dimensional.The model parameters and their reference values are described in Table S.1, except the forcing function F 0 (t) which are presented separately in Table 1.
While t denotes the non-dimensional time, we note that the characteristic time for the this model is 5 days, and hence, one time unit in this model corresponds to 5 days, as originally assumed by Lorenz (1984) [22].We refer to this as 1 [8].
The L84-S61 model is a nonlinear, non-autonomous system of ODEs [24].Using vector notation, this system can be written as where X = (X, Y, Z, T, S), and F(X, t) is a time-dependent, nonlinear vector function of X given by the right-hand side of Equations (2.1) to (2.5).Its solutions are bounded, i.e. ||X|| < C, with C being a positive constant.The system is conditionally dissipative, i.e. ∇ • F(X, t) < 0 under certain conditions, meaning that finite-volume attractors might exist.Despite being a simplified representation of the ocean-atmosphere system, the L84-S61 model retain some of main characteristics of a state-of-the-art ESM: it is nonlinear, multiscale, multicomponent, complex and chaotic.Hence, conceptual results obtained from this model can be insightful, if not informative, of general properties of ESMs.However, contrary to complex ESMs, which are high-dimensional (normally with billions or even trillions of degrees of freedom), the L84-S61 model consists of only 5 ODEs, making it an affordable model to be (extensively) studied computationally -in particular, allowing for very large ensembles to be run.

Numerical solver, parameter values and ensemble design
The L84-S61 model is solved using the 4th-order Runge-Kutta method, with time step 0.01 LTUs (1.2 hours).The output frequency is 0.2 LTUs (1 day).All results, whether single trajectories or ensembles, are presented as 1-year averages * * .
All simulations use the parameter values as shown in Tables 1 and S.1, except for some simulations in Sections 3 and 6, in which H = 0 and 0.0025, respectively.Regarding the forcing, note that the values presented in Table 1 means that the forcing oscillates seasonally around an * * Strictu sensu, averages are not a solution to the L84-S61 system's IVP.However, for the concepts and computational results presented in this paper, this difference is of little importance.In fact, for temperature and salinity, the difference between annual averages and actual values is small, and hence the latter can be used instead as proof of concept.For the atmosphere, it only matters if we were to look at observables where the annual average of the observable is very different from the observable of the annual average.Such a function would have to be nonlinear to begin with.The only point where the difference potentially matters is in the convergence time for the atmosphere.
average value F m = 7 with seasonal amplitude M = 1, while being driving to another value at a rate of H = 0.01 units per year, or 1 unit per 100 years.
The ensembles run in this work are designed as follows.Given an initial condition X 0 = (X 0,1 , X 0,2 , X 0,3 , X 0,4 , X 0,5 ) in the phase space, we randomly sample another 1,000 initial conditions such that, for each dependent variable, the sample is normally distributed around X 0,j with variance given by σ X0,j -with σ X0,j being two orders of magnitude lower than X 0,j .Hence, each ensemble has 1,001 members.The details of each individual experiment, including duration and parameter values, can be found in the Supplementary Materials.

The pullback attractor and convergence time
The pullback attractor [11,12] is a mathematical object that generalises the concept of attractor to non-autonomous dynamical systems.This approach consists on the idea that, for most non-autonomous systems, there exists a time-dependent object in the phase space, to which trajectories that started in the infinite past will converge.Such object presents therefore a natural distribution for the internal variability of the system.
A formal definition can be presented as follows.Let us denote the solution to the initial value problem (IVP) given by Equation (2.9) and X(t 0 ) = X 0 by X(t; X 0 , t 0 ), and the corresponding phase space by X.A set A = A(t) in the phase space is said to "pullback" attract a set, or ensemble of points for all t, where dist X (•, •) denotes the Hausdorff semi-distance between sets in the phase space.The time-dependent set A(t), if also invariant with respect to the dynamics, is called pullback attractor.When pullback attractors exist, there might also exists an invariant probability distribution supported on this set, so-called the pullback invariant measure (or distribution), which we will generically denote by µ A [14].An explicit, rigorous computation of both A(t) and µ A is only viable for very simple dynamical systems, and usually not possible for most nonlinear ones, including L84-S61.However, for nonconservative systems, a more practical approach is possible.This relies on the fact that, in general, a solution (or ensemble) starting near or on the attractor takes only a finite time to lose most of its dependency on the initial condition and run through (span) most of the attractor.The time for this convergence is dependent on the system and its relevant time scales, and can also be estimated numerically, as we shall see below.
Figures 1(a-c) illustrate this convergence to the pullback attractor for some of the variables of L84-S61.There, the pullback attractor and its natural distribution are computed from a micro ICE normally distributed around a central IC point X 0 in the attractor, with variance σ X0 being O(10 −2 ) for atmosphere variables, O(10 −3 ) for the ocean temperature and O(10 −4 ) for ocean salinity (as per Daron and Stainforth, 2013 [8]; see also Supplementary Materials).Note that, soon after the simulation starts, the initial micro cluster of trajectories disperses quickly and cover most of the attractor within a few years.The exact number of years depends on the variable of consideration though.For example, the time taken is visibly long for the ocean temperature (Figure 1(a)), and even longer for the salinity (Figure 1(b)), but very short for the fast, atmospheric variable X (Figure 1(c)).The latter is in line with what has been reported by Drotós et al. (2015) [17] and Tél et al. (2020) [10] for the L84 atmospheric model.

Convergence time
The convergence time, which we shall denote as t conv , can be loosely defined as the time taken by a localised ensemble to become indistinguishable from the pullback attractor.A statistically formal way to compute t conv is by comparing, at each instant of time, the distribution of interest with a snapshot of the numerically estimated pullback invariant distribution, via a hypothesis test using some suitable statistics, where the null hypothesis H 0 is that both distributions come from the same population.If we define a function of time h such that h(t) = 1 if the null hypothesis is rejected at time t and h(t) = 0 if not rejected, then we could define t conv such that In the definition above † † , there might exist t > t conv K such that h(t) = 1, which might put in question whether the convergence has been achieved.To avoid that, a statistically robust way to define t conv would be to take the distribution of h(t) in the time interval of consideration, repeat the experiment several times, and build the distribution of h(t) values for all those experiments, which can then be translated into a distribution of t values, with associated uncertainties.This resulting distribution should cluster around a value t that would be taken as t conv .
Both ways of estimating t conv are clearly dependent on the system of interest, as well as the initial condition and also the dynamical variable in question.Crucially, in practice, when dealing with computationally-generated distributions, such computation is also dependent on the size of the ensemble.There is not a unique way of doing it, and hence t conv is also dependent on the test used, as well as the significance level chosen.There are several ways to test this hypothesis [28].In this work, we use a two-sample Kolmogorov-Smirnov (KS) test [29].For two distributions P 1,n1 (x) and P 2,n2 (x) of sizes n 1 and n 2 respectively, the KS test is defined as For the KS test, null hypothesis H 0 should be rejected ‡ ‡ at significance at level α if D(P 1,n1 , P 2,n2 ) > C n1,n2,1−α , where C n1,n2,1−α can be found in [30].
We illustrate this approach by computing t conv for the spinup distribution shown in Fig- ures 1(a-c).To do so, we test H 0 with significance level α = 0.05, where the reference distribution is given by a 100,000 years single-trajectory solution starting from the same central IC (Supplementary Materials).This is presented in Figures 1(d-f), which shows that t conv is 90 years for salinity, 50 years for temperature, but only 19 years for atmosphere.The latter is substantially higher than what has been reported by Drotós et al. (2015) [17], which found a t conv of only 5 years for the L84 system, suggesting that the coupling with slow-mixing variables increases the relaxation period for the atmosphere variables in this context.
An alternative way to define a convergence time would be to assume that initially the statistic -in this case the KS statistic -D decays exponentially, such that D(t) ≈ D(t 0 ) exp −τ (t−t0) .In this case, such a convergence time could be taken as 1/τ .The characteristic decay exponent can be estimated by looking at the logarithm of D, which is presented in Figure S.1, and computing the angular coefficient of the straight line it approaches in the first few years of decay.This gives τ equals 0.0378, 0.264 and 0.1221 for T , S and X respectively.These correspond to estimated times of approximately 26 years, 38 years and 8 years respectively, which is roughly half the values † † Note that we opted to define tconv as normalised by K, so that the corresponding unit is year, instead of LTU.
‡ ‡ For convenience, in this work we use MATLAB's build-in function kstest2 instead.This function rejects the null hypothesis based on the p-value, and not by comparing the test statistic with a reference value.
of t conv estimated via Equation (3.2).Hence although quantitatively different, both approaches provide very similar information.

Caveats with the pullback attractor approach
The pullback approach has been proposed as an alternative way of defining climate: it gives a mathematically sound measure of the system's internal variability, and being time dependent, provides both a natural set of plausible states at each instant of time -the snapshot attractorand a natural probability distribution of events at each instant of time -the pullback invariant distribution.This has been discussed and illustrated by several authors [1,31], and has proven to be a more rigorous and useful definition of climate for long-term (e.g.IPCC-like) future scenarios.
This approach comes with some caveats though.By definition, the computation of such object requires an ensemble to be initialised in the infinite past, which is impractical from the computational point of view.In general, it is possible to approximately compute the attractor provided that the system is run for longer than t conv .But again, this is problematic, particularly in climate modelling: on one hand, some components of the Earth System evolve on long timescales of hundreds to thousands of years; on the other hand, anthropogenic, non-periodic forcing started only a couple of centuries ago.
Another caveat is that, while the pullback attractor represents all the internal variability of the mathematical model, it is known that only a few of these states can be representative of today's climate.Therefore, using the pullback attractor to measure "tomorrow's" climate might include a large number of unrealistic states -they are part of the internal dynamics of the model but not attainable within that time frame for a given initial condition.This will be discussed in the next section.

Micro initial condition ensembles and the evolution set
Although the pullback attractor provides a useful, mathematically sound definition for longterm climate (beyond the convergence time), it is less useful in quantifying the variability in the short-mid term (months to years, or even decades), when the intermittency of the dynamics is still dependent on the initial state of the system.This is because it overestimate the forecast uncertainty by allowing all possible states within the attractor, including those that do not reflect our knowledge of the present state of the system.
For example, considering the snapshot attractor for a given day (say "today"), it corresponds to a large range of possible values.But given sufficient information, it might be that only one of those states is possible (up to a certain level of residual uncertainty), so many of the states on the snapshot attractor are unrealistic given our knowledge of "today's" system.We also know that the climate today constrains the climate of tomorrow, in the range whose the system still carries the memory from the initial state -that therefore excludes a large portion of the pullback attractor.This means that any snapshots of the pullback attractor over-quantifies the variability and distorts the probability of events in the short and mid-term.This is illustrated in Figure 2, where we present the evolution of a micro ICE under climate change next to the evolution of the pullback ICE of Figure 1.This side-by-side comparison (see also Figure S.2 in the Supplementary Materials) shows that, in the first few decades, the pullback natural distribution, which is intrinsic to the mathematical system, over-represents climate uncertainty.Note that the evolution of the micro ICE is initially constrained to an smaller set, which is evolving over time, and seems to converge to the pullback attractor A(t) only after a few decades.For this reason, we name this as the evolution set E(t).For a given non-autonomous chaotic system, this set is solely dependent § § on the initial state X 0 , on the initial micro-uncertainty given by the variance σ X0 and on the initial time t 0 .Therefore, we shall denote the evolution set as E(t; X 0 , t 0 , σ X0 ).Some basic properties of this set are straightfowward.First, its existence is guaranteed by the existence and uniqueness for the IVP for (2.9).Second, by the definition, we have that E(t 0 ) = D X0 , the ICE set.Also by definition, we have that E(t; X 0 , t 0 , σ X0 ) −→ A(t) as t 0 −→ −∞.It also follows that, for an initial ensemble set within the pullback attractor, i.e.D X0 ⊆ A(t 0 ), we have that E(t) ⊆ A(t), for all t ≥ t 0 .In practice, when estimating both E(t) and A(t) numerically, these properties do not hold ipsis literis, and the design of the ensemble becomes quite important.We also note that, associated to E(t), this numerical example suggests the existence of a distribution µ E supported on this set, which we will assume to be true.Its relationship to the pullback invariant distribution µ A is not as clear though.
Climate modellers are familiar with the idea of exploring ICU using micro ICEs.Nevertheless, they are in general taken simply as an exploration of uncertainty, rather than the object we are trying to characterise.Here, we bring together the ideas of the pullback attractor with the methods applied in climate modelling and produce an attractor-like object which essentially represents future climate under climate change -which we called the evolution set.
The formalism above allow us to revisit the content of previous section, and reframe it in terms of "forward" convergence ¶ ¶ .There, the existence of a convergence time t conv might suggest that E(t) ≈ A(t) almost everywhere for t > t conv .Hence, the question is: does that really happen?Which conditions are necessary to prove that, for t >> t 0 : 1) E(t) and A(t) are sufficiently close * * * ; 2) µ A approximates µ E as t −→ t conv .If such statements are true, the pair (A, µ A ) would hold key mathematical information regarding future climate.
In the next section, we will explore some features of the evolution set by looking at its dependence on the initial conditions.

Macro initial condition uncertainty
Another issue related to the short-to-mid term climate prediction is the level of uncertainty of the actual state of the system in some variables.While small uncertainty can be covered by a micro IC ensemble, the uncertainty in the initial state of some variables might be of the same order of magnitude of the typical values for the variable itself, for instance if the initial state is based on a model spinup, or derived from the interpolation of sparse datasets, or even because of a lack of data.
From a climate prediction point of view, these are relevant, and macroscale variations in ocean quantities such as temperature and salinity, and atmospheric ICU have already been linked to decadal variations in regional climate in the Northern Hemisphere [7].The question is therefore how would such macro uncertainty impacts the evolution of the system, via its evolution E(t) set.§ § In practice, any numerical estimate of E will also depend on the size and shape of the initial ensemble.
¶ ¶ Here, we note again that, pullback attractors are not forward attractors in general.Although, under certain conditions, a pullback attractor could satisfy a weak form of forward convergence.The interested reader can find a detailed exposition of this in the section 9.5 of Kloeden and Yang (2020) [12].
* * * Note that two sets might be infinitely close but disjoint.For instance, the sets of rational and irrational points within the interval [0, 1] are disjoint but their closure (with respect to the standard topology) equals the full interval.

Macro ICU from a control simulation (single trajectory)
One of the sources of macro ICU is the potential to initiate climate ensembles from different states -including ocean states -from a long control run with an ESM.To illustrate this, we chose four different points in the attractor, all corresponding to a point in an existing trajectory after an initial 5,000 years spinup.For simplicity, we name these ICs by IC 1, IC 2, IC 3 and IC 4, with corresponding micro ICEs referred as ICE 1, ICE 2 and so on.Note that those ICs differ in all five dependent variables, and are illustrated in Figure 3 for the ocean variables.All ensemble distributions have the same variance, as noted in Section 2.2.
Figure 4 shows that, for the ocean variables, the dependence on the initial condition is significant.A first remark is that all four micro-initialised resulting distributions differ substantially from the pullback invariant distributions shown in Figures 2(a,b).Further to that, they are also different among themselves.For instance, in Figure 4(b), the micro ICE centred at IC 2 starting from a low temperature tends to decrease for a few years before increasing again, despite the monotonic increase in forcing.This is not followed by the micro ICE centred at the nearby IC 1, as shown in Figure 4(a) which spreads out very quickly after initialisation and is visually (increasing) monotonic from the beginning.In the case of IC 2, the decrease in temperature is accompanied by an initial increase in salinity as shown in Figure 4(f), which is then followed by a steady decrease.Nevertheless, in all four cases, the distributions seems to coincide after a few decades, becoming visually indistinguishable from each other and from the pullback invariant distribution (see also Figure S.3).
This macro ICU dependence has important consequences for climate prediction in seasonal to decadal time scales.A common practice in climate modelling is to start a simulation from initial conditions obtained from a spinup "control" run.This control run allows one to find the system's attractor, but does not resolve the uncertainty about where in the attractor one should start from.As we have seen, different micro ICEs could lead to different transient distributions, representing a different climate in the short-to-mid term -even if the initial condition is obtained from the same solution after spinup.

Macro ICU that reflect uncertainty in one variable
Another source of macro ICU is when initialising the model from observations, in which case the uncertainty in some variables could be orders of magnitude higher than others.As an example, if in-situ data is being used to initialise the model, it is possible that one might have measurements for one variable but not for others, for instance in case of defective equipment (e.g. via biofouling).In this case, the initial state of the variable is subjected to macro uncertainty.
This scenario is illustrated in Figure 5, where we highlighted four possible initial conditions, named by IC 5, IC 6, IC 7 and I 8, which are identical in the atmosphere variables X, Y, Z, but may differ in temperature and salinity (Supplementary Materials).For instance, IC 5 and IC 6 has identical temperature but differ in salinity; the converse is true for IC 6 and IC 7, and so on.
The sensitive to macro ICU with respect to a single variable is illustrated in Figure 6, which shows the results for micro ICEs starting from the ICs indicated in Figure 5.Note that macro uncertainty in salinity does not seen to alter the evolution set and its distribution, as indicated in Figure 6(a,b).On the other hand, macro uncertainty in temperature has a significant effect on salinity, as shown in Figures 6(e,g): the evolution set ans its distribution for salinity are significantly different, despite having ensembles around the same initial salinity state (see also Figure S.4).
This sensitivity of both E and µ E to macro ICU in a single, slow variable is remarkable, and suggests that a proper quantification of the uncertainty in future climate projections requires an assessment of macro ICU as well.

Convergence time and macro ICU
As macro ICU impact the evolution set and its distribution, one might asks whether the convergence time t conv is also affected by it.Here we revisit the concept of convergence time and show how it can vary with in a macro ICU scenario.We illustrate this by computing t conv , using Equation (3.2), for the eight micro ICEs shown in Figure 3 and Figure 5.The resulting t conv and corresponding evolution of the KS statistics are shown in Figure 7 for the ocean variables.
When starting from a control run trajectory (as per Figure 3), the resulting t conv can vary dramatically.This is shown in Figures 7(a,c).First, IC 1 provides a short t conv for both temperature and salinity, being of 14 years and 46 years, respectively.This t conv increases substantially from IC 1 to IC 2, being of 30 years for temperature and 70 for salinity.For IC 3 and IC 4, while the t conv for temperature remains of the same order (34 and 32 years, respectively), it still varies substantially for salinity, resulting in a t conv of 101 years for IC 3 and and 92 years for IC 4. We also note that the order of t conv is the same for both variables in this case: IC 1 has shortest t conv for both temperature and salinity, IC 2 is second, and so on.
When starting from chosen-values within the attractor (as per Figure 5), the results are rather different.This is shown in Figures 7(b,d).In particular, both variability and and order of t conv differs from those shown in Figures 7(a,c).For instance, the variability in t conv is 27 to 34 years for temperature (instead of 14 to 34 years) and 69 to 112 years (instead of 46 to 112) for salinity.Also, the shortest t conv for temperature (27 years) is given by IC 5, while the shortest for salinity is given by IC 6 (69 years).
The starkest contrast is observed when comparing IC 6 and IC 8.Note that, while both ICs have the same value of salinity, their respective micro ICEs have a t conv that differs by 43 years, highlighting the impact that macro ICU in a single variable (in this case ocean temperature) can have in other variables.
In Herein et al. ( 2016) [18], where the authors only looked at uncertainty using an intermediate-complexity model, they noted that t conv did not chance for micro ICEs starting at different instants of time.However, their micro ICE was generated perturbing only one variable (the surface pressure field), keeping the others equal for all ensemble members, while the results correspond to another variable, the annual mean surface temperature in a single grid point of the model (what they called a small scale) located within continental Europe.While treated there as a simple approximation, the results presented here for a much simpler model suggests that uncertainty in other variables can have a significant impact on the distribution and its convergence time, as the response from the initial uncertainty in slow variables could take a while to show.

How rate of change in forcing affects the uncertainty of climate predictions
In the context of climate change, the system is under an external forcing (e.g.change in temperature due to anthopogenic carbon dioxide emissions) that is both dynamic and uncertain.Those uncertainties are usually investigated via scenarios, which in the context of IPCC, have shown to dramatically affect the climatology predicted by CMIP models.In the context of this work, this external forcing uncertainty may also affect the evolution of an ICE as a distribution.We illustrate this by looking at the evolution of the micro ICE centred in IC 2 (shown in Figure 3) but under a slower rate of "climate" change regime.Here, we reduce the rate of change in forcing H by a quarter, from H = 0.01 to H = 0.0025, meaning that it now takes 400 years for the baseline forcing F m to increase by one unit.The resulting time series are shown in Figure 8 for ocean temperature and salinity, where we also included the H = 0.01 time series for reference.Changing, or in this case reducing, the speed of climate change has important effects on the resulting distributions.While the ICE distributions in Figure 8(c) shows a mildly monotonic decrease, Figure 8(a) shows that this behaviour, while kept, is much more pronounced under a weaker forcing.As this slowly changing distribution evolves, it again shows a distinct behaviour at around year 120: the distribution suddenly gets broader, with the temperature of several ensemble members decreasing sharply.About 40 to 50 years later, the ensemble narrows again and regain a shape akin to that of Figure 8(c).These behaviour are mirrored by the salinity distribution, as shown in Figures 8(b,d).
These curious behaviour, which is consequence solely of altering the rate of change in forcing, can be better seen when looking at the projection of the phase space onto the ocean variables subspace, as shown in Figure 9.At a faster climate change rate, shown in Figure 9

Conclusions
This article discussed several aspects related to the climate predictability in short-and mid-time scales, including annual to multi-decadal.To do so, we introduced the idea of an evolution set, where we combined the concepts of pullback attractor and micro ICU to produce an object lying within the system's pullback attractor whose shape is constrained by a more refined knowledge of the initial state of the system -via a micro ICE.While the evolution set is usually contained in the pullback attractor set, the latter is much larger, and their respective distributions, or climate projections, are different.
In addition to that, we attempted at defining a convergence time, as the time taken for an ICE distribution to become indistinguishable from the pullback invariant distribution.We also explored micro and macro ICU, revisited the concept of pullback attractor, and discussed the influence of those in the evolution set and the convergence time.We also discussed the effect of different rates of change in forcing in the evolution set.Given the significant differences produced, these results suggest that all these aspects should be considered when designing ensembles for chaotic, non-autonomous systems, in particular for ESMs in a climate-change scenario -i.e.under non-periodic external forcing.
Although the results obtained are dependent on the particular low-dimensional model used, the ideas are model-independent and should be applicable to any chaotic non-autonomous system.This includes the concepts of evolution set, micro and macro ICU and convergence time.
From a theoretical point of view, this work leaves many questions to be answered, which we believe to be of both mathematical and climate science relevance.The first set of questions relate to the evolution set E = E(t; t 0 , X 0 , σ X0 ): • Is it possible to prove rigorous results regarding the sensitivity and dependence of E to the central IC X 0 , initial time t 0 and variance σ X0 ?
• What is the relationship of E to the pullback attractor A? Is there any other relationship beyond E ⊆ A when D X0 ⊆ A(t 0 )?
• How many ensemble members are needed to characterise E for a given X 0 , t 0 and uncertainty as measured by σ X0 ?
• How dependent is E on the shape of the ICE?For instance, would a non-Gaussian distribution lead to a very different E?
• How does the distribution µ E relates to the pullback invariant distribution µ A of the pullback attractor?
Another important question is how does uncertainty in one variable propagates, or rather influence others?For example, we saw that macro ICU in temperature seems to greatly affect salinity, but the converse is not true.
A final but more ambitious question relates to the "size" of attractors, as illustrated in Figure 10.How large are attractors in ESMs?In other words: • Is it possible to estimate their shape without resourcing to brute force, given the computational limitations of running such models?
A final note on the evolution set E is that there can be many, depending on what observations ones uses to constrain the possible climate scenarios with.The same applies to the the evolution distribution µ E .In this practical sense, the central IC and variance used in the definition of E are just fudges to simulate the residual uncertainty after the information from the observations has been brought in.So the questions above, although generically formulated, might be asked in relation to an E constructed from assimilating some observation into a more realistic climate model for example.Nevertheless, answers to those questions would be a valuable resource in the design of relevant and influential climate model ensembles.

Experiments
Here we describe in detail the experimental design of each simulation.The following holds for all experiments, unless otherwise noted: • Length of simulation: 200 years, except experiments 1, 2, 3, 13 and 14.
The ICs vary across experiments, and are presented in detail below.

Figure 1 :
Figure 1: Left column: Pullback attractor and its natural distribution for L84-S61 computed from a 500 years micro ICE simulation, where green solid line shows the numerical solution starting from the central IC.Right column: Corresponding convergence time computed using Equation (3.2) and the KS statistic based on a 100,000 single trajectory simulation.(a,d) ocean temperature; (b,e) ocean salinity; (c,f) atmosphere variable X (intensity of westerly wind).

Figure 2 :
Figure 2: Comparing the pullback invariant distribution with the distribution generated by a micro ICE, with H = 0.01 in the first 100 years, and H = 0 in the remaining 100 years.Left column shows the evolution of an ensemble which initially covers the entire pullback attractor.Right column shows the evolution of a micro ICE.Panels (a-f) show: (a,d) ocean temperature; (b,e) ocean salinity; (c,f) atmosphere variable X (intensity of westerly wind).

Figure 3 :
Figure 3: Attractor for the system L84-S61 with H = 0 when F m = 7 (blue) and F m = 8 (red) projected on the ocean temperature-salinity (T, S) subspace.The black dots on the F m = 7 attractor indicated the location of ICs 1 to 4.

Figure 4 :
Figure 4: Macro ICU from a control run simulation: comparing the evolution set and distribution of the slow-mixing ocean variables for different micro ICEs in a macro ICU scenario, with H = 0.01 in the first 100 years, and H = 0 in the remaining 100 years.Left column shows the ocean temperature.Right column shows ocean salinity.Panels (a-f) show: (a,e) IC 1; (b,f) IC 2; (c,g) IC 3; and (d,h) IC 4.

Figure 5 :
Figure 5: Attractor for the system L84-S61 with H = 0 when F m = 7 (blue) and F m = 8 (red) projected on the ocean temperature-salinity (T, S) subspace.The black dots on the F m = 7 attractor indicated the location of ICs 5 to 8.

Figure 6 :
Figure 6: Macro ICU in a single variable: comparing the evolution set and distribution of the slow-mixing ocean variables for different micro ICEs in a macro ICU scenario, with H = 0.01 in the first 100 years, and H = 0 in the remaining 100 years.Left column shows the ocean temperature.Right column shows ocean salinity.Panels (a-f) show: (a,e) IC 5; (b,f) IC 6; (c,g) IC 7; and (d,h) IC 8.

Figure 7 :
Figure 7: Distance between the micro ICE distributions to the pullback invariant distribution, measured through the KS statistics (solid lines), and convergence time (dashed-dot lines) computed using Equation (3.2): (a,b) ocean temperature; (c,d) ocean salinity, for the micro ICEs centred at: ICs 1 to 4 (left column) as per Figure 4; ICs 5 to 8, as per Figure 6.

Figure 8 :
Figure 8: ICE distributions starting from IC 2 in Figure 3, for H = 0.01 (100 years of climate change, followed by 100 years of non-forced climate with F m = 8) shown in the upper panels, and H = 0.0025 (400 years of climate change, followed by 100 years of non-forced climate with F m = 8) in the bottom panels.Left column panels show temperature for (a) H = 0.01; (b) H = 0.0025.Right column panels show salinity for (c) H = 0.01; (d) H = 0.0025.
(b), the distribution seems to have less freedom to explore the phase space and has its way forced into the attractor F m = 8.At a slower climate change rate, presented in Figure 9(a), the ensemble members have now more freedom -and time -to explore the phase space and any intermediate attractors between those of F m = 7 and F m = 8.As suggested by Figure 10, one of those intermediate attractors is somehow broader (in the ocean variables) than the neighbour ones, and trajectories entering there might eventually reach (time allowing) lower values of temperature and higher values of salinity.

Figure 9 :
Figure 9: Projection of the phase space onto the (T, S) subspace, with a heatmap indicating the number of ensemble members that passes through each point at least once (no repetitions are counted): (a) H = 0.01; (b) H = 0.0025.These correspond to the joint distributions shown in Figure 8(a,c) and Figure 8(b,d), respectively.

Figure 10 :
Figure 10: Attractor for the non-forced L84-S61, projected over the ocean temperature-salinity (T, S) subspace, for several values of F m between 7 and 7.5.All attractors shown correspond to a single trajectory starting from the same IC (black dots).

Figure S. 2 :
Figure S.2: Comparing the pullback invariant distribution with the evolution distribution generated by a micro ICE, with H = 0.01 in the first 100 years, and H = 0 in the remaining 100 years.Solid line shows the ensemble mean, and the shade shows 1 standard deviation from the mean.In blue is shown an ensemble that initially covers the entire pullback attractor.In red is shown the evolution of a micro ICE.Panels correspond to individual variables: (a) atmosphere X; (b) atmosphere Y ; (c) atmosphere Z; (d) temperature; (e) salinity.

Figure S. 3 :
Figure S.3: Macro ICU from a control run simulation: comparing the evolution set and distribution of the slow-mixing ocean variables for different micro ICEs in a macro ICU scenario, with H = 0.01 in the first 100 years, and H = 0 in the remaining 100 years.Solid line shows the ensemble mean, and the shade shows 1 standard deviation from the mean.Results for IC 1, IC 2, IC 3 and IC 4 are presented in blue, red, yellow and magenta, respectively.Panels correspond to individual variables: (a) atmosphere X; (b) atmosphere Y ; (c) atmosphere Z; (d) temperature; (e) salinity.

Figure S. 4 :
Figure S.4: Macro ICU in a single variable: comparing the evolution set and distribution of the slow-mixing ocean variables for different micro ICEs in a macro ICU scenario, with H = 0.01 in the first 100 years, and H = 0 in the remaining 100 years.Solid line shows the ensemble mean, and the shade shows 1 standard deviation from the mean.Results for IC 5, IC 6, IC 7 and IC 8 are presented in blue, red, yellow and magenta, respectively.Panels correspond to individual variables: (a) atmosphere X; (b) atmosphere Y ; (c) atmosphere Z; (d) temperature; (e) salinity.

Table 1 :
Description of the parameters and their reference values in the forcing function F 0 (t), as per Daron andStainforth (2013) The table below contains all model parameters and their values.The values shown were the same in all simulations.Coefficient of internal diffusion in the ocean k a 1.8 • 10 −4 Coefficient of heat exchange between the ocean and atmosphere ω 1.3 • 10 −4 Coefficient derived from the linearised equation of state ϵ 1.1 • 10 −3 Coefficient derived from the linearised equation of state γ 0 7.8 • 10 −7 Coefficient for the atmospheric water transport γ 1 9.6 • 10 −8 Coupling parameter for the wind dependent atmospheric water transport Table S.1: Description of the parameters and their reference values used in the L84-S61 model, as per Daron and Stainforth (2013) [1] .