## Abstract

The theory of chaotic dynamical systems gives many tools that can be used in climate studies. The widely used ones are the Lyapunov exponents, the Kolmogorov entropy and the attractor dimension characterizing global quantities of a system. Another potentially useful tool from dynamical system theory arises from the fact that the local analysis of a system probability distribution function (PDF) can be accomplished by using a procedure that involves an expansion in terms of unstable periodic orbits (UPOs). The system measure (or its statistical characteristics) is approximated as a weighted sum over the orbits. The weights are inversely proportional to the orbit instability characteristics so that the least unstable orbits make larger contributions to the PDF. Consequently, one can expect some relationship between the least unstable orbits and the local maxima of the system PDF. As a result, the most probable system trajectories (or ‘circulation regimes’ in some sense) may be explained in terms of orbits. For the special classes of chaotic dynamical systems, there is a strict theory guaranteeing the accuracy of this approach. However, a typical atmospheric model may not qualify for these theorems. In our study, we will try to apply the idea of UPO expansion to the simple atmospheric system based on the barotropic vorticity equation of the sphere. We will check how well orbits approximate the system attractor, its statistical characteristics and PDF. The connection of the most probable states of the system with the least unstable periodic orbits will also be analysed.

## 1. Introduction

After Lorenz’s discovery of deterministic chaos [1], dynamical system theory has played an increasing role in the theory of climate. The widely accepted fact now is that a typical atmospheric system, as the famous Lorenz-63 model, is dissipative (i.e. it is contracting in phase space) and chaotic (system trajectories are sensitive with respect to initial conditions). The system energy is bounded and it evolves in a bounded set. These two facts guarantee (for finite-dimensional systems) the existence of a non-trivial invariant attracting a fractal set inside the system phase space. This set is called ‘the attractor’ and it was introduced to dynamical system theory by Ladyzhenskaya [2]. The attractor is one of the most important properties of the system since it attracts all system trajectories and the system motion occurs inside the attractor neighbourhood. This type of behaviour is well known for the Lorenz-63 system. Many atmospheric systems, including barotropic [3], two-layer quasi-geostrophic [4] and primitive atmospheric systems [5,6] also have this attractor existence property.

An important property of an attractor is that a stationary probability distribution function (PDF) could be defined on it. This PDF allows calculations of averages along system trajectories as well as statistical characteristics of the system, giving a strong theoretical basis for climate modelling. Dynamical system theory provides additional important information about such systems. Lyapunov exponents (local and global) give the measure of instability of system trajectories; attractor dimension is a general measure for the complexity of system dynamics; information entropy is a useful tool for analysing PDFs on the system attractors, etc.

The next logical step in understanding and quantifying attractors of atmospheric models would be construction of the corresponding PDF, giving all information about system statistical characteristics. Unfortunately, these attractors are fractals with large dimensionality, so it is not possible to perform this sort of calculation. The only possible way to describe such a complex object should be some kind of an approximation method. Some regions on the attractor could be visited by the system more frequently than others, and the PDF will have a local maximum in these areas. It is natural to try to approximate the PDF using these features as building blocks. Ideally, this approximation procedure should be based on the analysis of the system equations rather than on *a posteriori* analysis of the (complex) system PDF in multi-dimensional phase space. The first natural question here is what objects are good candidates for the local PDF maxima?

In 1979, Charney & De Vore [7] proposed the idea that these building blocks could be unstable stationary points of the system. Indeed, in the vicinity of the stationary point, the right-hand side of the system is close to zero, the system evolves slowly, so that the PDF of the system could have a local maximum. This approach was widely accepted in 1980–2000 [8–11]. These studies showed that inside the attractors of some simplified atmospheric models, there exists a small number of slightly unstable stationary points. The system trajectory spends a lot of time in the vicinity of these points. As a result, the so-called quasi-stationary ‘weather regimes’ could be defined. The average lifetime in the neighbourhood of the stationary point is determined by the stability characteristics of this point. This fact allows one to classify the states on the attractor of the system and to connect the predictability of the system to the stability of the obtained stationary solutions [12,13]. All of the above results were obtained for weakly chaotic atmospheric models with attractor dimensions not exceeding values of 3–6.

Unfortunately, extension of this theory to more complex systems with multi-dimensional attractors was not very successful [14]. The problem is that in the case of a strongly chaotic system, its PDF becomes close to a normal PDF and the stationary points become highly unstable and generally do not define any regimes. As an example, we can refer to results of papers [15,16] with the analysis of the regime behaviour for a three-level quasi geostrophic (QG) system. In the completely baroclinic case, the authors were able to calculate only one stationary solution that does not determine any regime. A bimodal PDF with a large number of stationary solutions existed only for unrealistically high dissipation values.

There have been many published studies devoted to the determination of weather regimes for general circulation models (GCMs) and the real climate system. The situation here is of course not as clear as for simple atmospheric models. From our point of view, it is very probable that such regimes exist (see [17–19] for discussion), but it looks like there are no sharp maxima in the PDF of GCMs. Moreover, though Berner & Branstator [18] have related features in GCM PDFs to features in the GCM’s average trajectories, we do not know of any studies that have related PDF maxima to the stationary solutions of a GCM. Finally, investigators have identified only a few regimes, so that it is impossible to approximate the attractors of modern GCMs by steady regimes with much accuracy. The possible explanation for this problem could be as follows. The multi-dimensional highly chaotic atmospheric systems may have just a few (or zero) slightly unstable stationary points. Consequently, it is impossible to approximate the system attractor using them.

There exists another method of approximating chaotic attractors based on the results of dynamical system theory, namely the possibility to reconstruct a system PDF using unstable periodic orbits (UPOs) of the system [20–22]. The strong chaoticity (i.e. hyperbolicity or Axiom A property) of the system is usually required for this method to be theoretically justified (see [22,23] for further details). For such systems, the set of periodic orbits is dense on the attractor of the system and the explicit expression for the system PDF via UPO expansion could be obtained. As a result, any arbitrary trajectory of the system can be approximated with any prescribed accuracy by these periodic orbits, and all system statistical characteristics (including the probability to find the system trajectory near a given state) could be calculated using UPOs. First application of this UPO expansion strategy was performed in [24] for a very simple dynamical system.

For systems of atmospheric dynamics, it is not possible to verify conditions guaranteeing correctness of the UPO expansion procedure. Atmospheric models are likely to have non-zero Lyapunov exponents, and the generalization of the UPO expansion theory for systems with non-zero Lyapunov exponents has not yet been obtained. Following the ideas of Gallavotti (‘chaotic hypothesis’ [25]), we will suppose that for typical multi-dimensional chaotic atmospheric systems, UPO expansion may still work, at least if we are going to calculate system macroscopic (average) characteristics, and apply the UPO expansion numerically. This approach is supported by the experimental results for the Lorenz system [26], as well as for some more sophisticated hydrodynamical systems [27–29].

In some sense, it is much more advantageous to use UPOs instead of stationary points for approximation purposes. First of all, the number of periodic solutions is much larger than the number of stationary points (a stationary point is just a special case of a periodic solution). With UPOs, one can also approximate time-dependent circulation regimes and oscillations as well as transitions between quasi-stationary states (if any) of the system. For example, in the study of Selten & Branstator [17], the authors raise the idea that the transitions between the regimes in a quasi-geostrophic atmospheric model occur along the unstable periodic solutions of the model. Gritsun [30] has shown that the leading complex empirical orthogonal functions (EOFs) in the T12 barotropic system are defined by several of the least UPOs of the system.

An important consequence of the UPO expansion formula is that the least unstable orbits make larger contributions to the PDF. As a result, the probability to see the system near a given state of the phase space has increased values in the vicinity of these orbits. Consequently, the most probable states (MPSs) of the system attractor could be located near the least unstable UPOs (provided, of course, the UPO expansion formula is valid for the system). If the system slows down near one of the MPSs, we may expect to obtain a corresponding ‘circulation regime’ in the traditional sense.

Taking all these facts into account, one can conclude that it is reasonable to expect that UPOs could play an important role in quantifying and simplifying the complex dynamics of highly chaotic and highly nonlinear multi-dimensional atmospheric systems. In this paper, we will try to develop this approach for examples of two atmospheric models based on the barotropic vorticity equation of the sphere. The paper is structured as follows. In §2, we formulate our models and present their basic statistical characteristics. Section 3 is devoted to the details of numerical methods for calculating UPOs, and some basic information about the UPOs of barotropic models are provided. The UPO-related results from dynamical system theory are discussed in §4. We will also present results on how the UPOs could be used to approximate model statistical characteristics and model circulation regimes and MPSs here. Our results are summarized in §5 of the paper.

## 2. Barotropic models of the atmosphere

In this section, we will describe two atmospheric models and present their basic statistical and instability characteristics. The problem of existence of the model circulation regimes will also be studied.

Both models are based on the barotropic vorticity equation of the rotating sphere,
2.1
Here, *ψ*=*ψ*(*Θ*,*φ*) is the dimensionless streamfunction (the length scale is the Earth’s radius and the time scale is the inverse Earth’s angular velocity, *Θ* is latitude and *φ* is longitude), the Coriolis parameter is *l*=*l*(*Θ*), the orography is *H*, the orography normalization is *k* (the orography scale is 10 000 m), the external time-independent forcing is *f*_{ext} and the turbulent viscosity and boundary layer friction coefficients are *μ* and *α*, respectively. The Laplace and Jacobi operators on the sphere are *Δ* and *J*.

We use the Galerkin method for spatial approximation of (2.1) with the set of basic functions being equatorially antisymmetric spherical harmonics. For the first model, we use T12 spatial truncation (i.e. the azimuth indices of spherical harmonics do not exceed 12). The system dimension equals 78. The second model has T21 truncation, its phase space dimension is 231. The Matsuno scheme was used for time approximation in both cases.

Model parameters were set in the following way. We used the real orography of the Earth’s Northern Hemisphere as *H* truncated to T12 or T21, respectively (*k*=0.14 in both cases). The methodology for constructing the external forcing for the system *f*_{ext} duplicated the iterative procedure from [31] and included the substitution of 300 hPa streamfunction January daily data from the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis archive into the right-hand side of (2.1) and calculation of the residual term (see [31] for details). As a result, the models should be able to reproduce the overall average state and variability, but cannot be expected to match reality in any more detail. Boundary-layer damping *α* was set to 1/25 days in dimension units, whereas *μ*=7×10^{−5}. (For the shortest modes in the T21 system, this corresponds to a characteristic damping scale of 5 days.)

After spinup of 10^{4} days from zero initial conditions, the models were integrated for 10^{6} days in order to obtain their statistical characteristics. In figure 1, we show the average states and standard deviations for T12 (figure 1*c*,*f*) and T21 (figure 1*b*,*e*) models together with the average daily January–February 300 hPa streamfunction from the NCEP/NCAR reanalysis data (years 1948–2010 were used in this calculation). It can be seen that both models reproduce the average state of the climate system with good accuracy. The standard deviation of the T21 system is also close to that of the data (with 30% weaker variability over the oceans), while the T12 variability does not have the correct structure and is about twice as weak as it should be.

Two leading modes of models/data variability patterns (leading EOFs) are shown in figure 2. In the case of the T21 system (figure 2*b*,*e*), we have a good agreement in the amount of variability explained by two leading EOFs (21/6% for the data and 19/10% for the T21 model), though the structures of the EOFs do not agree particularly well. The T12 model (figure 2*c*,*f*) is not very successful in reproducing leading EOFs. However, both its first and second EOFs are physically reasonable, representing the Arctic oscillation and the Pacific/North America pattern that appear as leading EOFs for 500 hPa NCEP/NCAR streamfunction fields [32]. Note also that the structure of variability of our T12 system is in fact very close to that of the T21 model from [31] and the quasi-geostrophic model T21QGL3 of [17].

It is worthwhile mentioning that both models are chaotic. To demonstrate this fact, we calculated the Lyapunov spectrum for both models with a technique based on the Oseledets theorem (described in [33]). The T21 system is obviously much more chaotic. The largest Lyapunov exponent decorrelation time (Lyapunov time) is 6 days, the model has 27 positive Lyapunov exponents *λk* and its Kaplan–Yorke attractor dimension equals 65. The smaller system has just six positive Lyapunov exponents (the Lyapunov time is 16.5 days), with the value of attractor dimension 13.5.

Two-dimensional PDF projections on leading EOF planes calculated with the box-counting algorithm are shown in figures 3 and 4 (T12 model, figure 3*a*–*d*; T21 model, figure 4*a*–*d*). None of these projections has any evidence of the ‘multiple regimes’ behaviour. The two-dimensional PDF projections for the T12 system look less Gaussian like, as could be expected for less chaotic systems. Note that EOF1–EOF2 PDF projection of the T12 model is close to that of the Selten & Branstator [17] model (note that our EOF1 and EOF2 are opposite to Selten & Branstator [17]).

As mentioned in §1, some places of the system attractor may be more visited by the system trajectory than others (i.e. the probability of the system to stay near these places is high). Hereafter, we will call these points MPSs of the system. In the same manner, one can define the most probable trajectory (MPT) segment as the part of the system trajectory with its *ε*-neighbourhood containing most of the system states (we will use *L*_{2} streamfunction metrics everywhere in the paper to calculate distances in the phase space). The leading MPT plays the role of the preferred path of the system evolution.

A possible way to find MPSs is to directly calculate the probability that the system trajectory visits test regions of the system phase space (ideally, the collection of the test regions should cover all the system attractor). The regions with the highest probability will be the MPSs. The natural way to choose a set of test regions is to calculate a long enough system trajectory (‘test trajectory’) and consider as many of the system states on it (say every 24 h) as potential MPSs. More precisely, one considers an *ε*-neighbourhood centred on each of these states and calculates the number of moments when the system trajectory is inside it, deriving a corresponding probability. The test state with the highest probability is MPS1.

Following this approach, we calculated the 3×10^{5} day trajectory of the T12 model (after 1×10^{4} day spinup, the system is supposed to be on the attractor) and take 3×10^{5} points on it serving us as test regions. Next, we calculated another 1×10^{6} day trajectory (‘dataset1’). For each of 3×10^{5} points, we calculated the number of states from ‘dataset1’ being inside its *ε*-neighbourhood. The point with the highest probability is MPS1. The *K*th MPS is the point with the highest probability separated from all previously chosen *K*−1 MPSs by at least *ε* (since we would like to have MPSs separated from each other). The stability of the calculated MPSs was verified by recalculation of the above probabilities using another independent 1×10^{6} day trajectory showing that at least the five leading MPSs change only slightly. Results of the calculations are stable in a wide range of *ε* as well.

The MPSs obtained by this procedure are not the same as traditional ‘circulation regimes’ since the system may not slow down in the vicinity of MPSs, and MPSs may not be a strict maximum of the PDF (except the leading MPS). On the other hand, calculation of the MPSs could be accomplished in the full-dimensional phase space while ‘circulation regimes’ are usually calculated in its low-dimensional subspace.

To calculate the leading 30 day MPT, we need to split the ‘test trajectory’ into 30 day pieces and use them (their *ε*-neighbourhoods) as the new test regions. Repeating the above calculations with these new 10 000 test regions (in our analysis, we actually used 30 000 from the extended model run) and ‘dataset1’, one can obtain MPT1.

Results of our calculations show that MPSs and MPTs exist for both models in spite of the fact that PDF projections (i.e. figures 3*a*–*d* and 4*a*–*d*) suggest no ‘regime-like’ behaviour. With *ε*^{2}=2×10^{−6} (i.e. about 15% of the attractor diameter), the leading five MPSs of the T12 system case are responsible for 1.8, 1.4, 1.3, 1.2 and 0.95% of all system states (for leading MPTs, corresponding values are 7, 4.5, 3.7, 3.2 and 2.8%). Close values were obtained for the T21 model. As a result, the density of system states near both MPT1 and MPS1 is 20 times larger than the average. It is also interesting to see that this density is very sensitive with respect to the position in the phase space (i.e. MPS1 captures twice as many states as MPS5).

The evolution of MPT1 and the position of 12 leading MPSs of the T12 model projected onto the EOF1–EOF2 plane are shown in figure 5*a*. Leading MPSs are concentrated in an area of PDF maximum (figure 3*a*). The MPT1 looks like a system motion through the area populated by these leading MPSs connecting some of them. Five leading MPSs in the physical space are shown in figure 5*b*–*f*. One can see that leading MPSs have many common features with circulation regimes obtained in [17] and [34] (figure 5*b*–*f* should be rotated 90^{°} clockwise and zoomed to be compared with fig. 9*b* of [17]). MPS1 has negative values in the Arctic with a minimum over Alaska and positive anomalies over the US, Pacific and Europe, and has close structure with regime R of [17]. The structure of MPS3 reflects that of regime R of [34]. The MPS5 (minimum over Greenland and dipole structure in the North Pacific) could be identified as regime G from [17] and [34]. MPS4 has several common features with regime A from [17]. It should be pointed out however that the model studied here and the definition of MPSs are different from those of [17] and [34], and it is not reasonable to expect a close match between MPSs and regimes.

The dynamics of MPT1 in the physical space is shown in figure 6*a*–*d*,*i*,*j*. System states at day 4 (figure 6*b*), day 10 (figure 6*c*) and day 20 (figure 6*i*) are close to MPS4, MPS1 and MPS5, respectively, suggesting that MPT1 connects MPS4 with MPS1 and MPS5. The transition time from MPS4 (regime A) to MPS1 (regime R) is about 6 days, while it takes 10 days for the system to get to MPS5 (regime G) from MPS1. This is close to the results of [17], who estimated most probable A–R and R–G transition times to be 6 days, and suggested that the transition from A to G should pass through the R regime.

In §5, we will try to relate this behaviour to the UPOs and stationary points of the models.

## 3. Unstable periodic orbits of barotropic atmospheric models

In this section, we will briefly describe numerical methods used for calculating UPOs of T12 and T21 systems (most of the technical details can be found in appendix A). Some general results on UPOs in the T12 and T21 systems will be presented.

By definition, a periodic orbit with period *T* is a system trajectory returning into the same position in the phase space after time *T*. Write the system equation in the resolved form *ψ*(*t*)=*S*(*t*,*ψ*_{0}). As a result, we get the following equation (the system is supposed to be autonomous) for UPO:
3.1
The numerical procedure for finding UPOs is based on equation (3.1), which is a system of nonlinear equations with respect to the initial position of the orbit and its period (for autonomous systems such as (2.1), an initial position of the system and an integration time uniquely define the system trajectory).

There are, of course, many different numerical methods for solving a system of nonlinear equations. Based on our early experience of calculating periodic orbits for the T12 system [35], we used the damped Newton method with two-dimensional tensor correction terms most of the time (details are described in appendix A), however, other techniques (i.e. minimization procedures based on limited memory Broyden–Fletcher–Goldfarb–Shanno method or generalized minimal residual method schemes [36–39]) could also be very effective [28,27].

A 600 Intel Xeon E5450 processor of the cluster computer MVS-100K of the Joint Supercomputer Center, Russian Academy of Sciences, was used in our calculations. As a result of huge numerical work, we found 2322 UPOs and 50 stationary points for the T12 system and 1554 UPOs and 320 stationary points for the T21 system. While the calculation of orbits for the T12 model is relatively simple, with the T21 system, the situation is much more problematic since the T21 system is much more unstable and has three times larger dimensionality. These difficulties were also mentioned in [40], namely that the system-enhanced instability and large dimensionality make it harder to find good initial guesses from the system trajectory. Another important reason contributing to the overall performance of the numerical algorithm is the fact that the T21 model has much more near-to-zero Lyupunov exponents (slowing down convergence of the Newton methods even more). As a result of all these numerical problems, we failed to calculate enough UPOs for the T21 system, especially for those having large periods.

Let us first look at the T12 orbits. Distribution of the system trajectory (light grey colour), the 20 least unstable UPOs (black colour) and the five least unstable stationary points projected onto the EOF1 and EOF2 planes are shown in figure 7*b*. Several least unstable orbits are situated in the area rarely visited by the system trajectory (i.e. these orbits are invisible from the observer standpoint). But most of the least unstable orbits and one stationary point are situated in the region near the PDF maximum (cf. figure 3*a*). The 150 shortest orbits and the 110 most unstable orbits of the T12 system are shown in figure 7*a*,*c*, respectively (numbers of orbits were chosen to match the total length of the least unstable UPOs in figure 7*b*). In both cases, the region of PDF maximum is rarely visited by these UPOs. It is clear from figure 7 that different parts of the attractor could have very different instability characteristics. To obtain decent PDF approximations, one needs to use some long orbits in addition to the shortest ones; choosing a few least unstable orbits may not approximate the attractor so well.

The set of 500 least unstable orbits ‘visually’ covers all system trajectories giving the hope that UPOs can, indeed, approximate the system attractor. Also, as suggested by figure 7*b*, some of the least unstable orbits are outside the trajectory attractor and should not be used in the attractor approximation procedure.

Distributions for the number of positive Lyapunov exponents *λ*^{+}, orbit Lyapunov (Kaplan–Yorke) dimensions and the growth rate for the number of periodic points with given periods are shown in figure 8. The T12 orbits have from 2 to 18 positive exponents (all orbits and steady states are unstable). The distribution peaks at 7, which is close to the number of *λ*^{+} of the T12 system (equals to six). A similar correspondence could be seen for the orbit Lyapunov dimension (the system attractor dimension is 13.5 and distribution maximum is 14).

The growth rate for the number of periodic points with given period *Γ*(*T*+*Δ*)/*Γ*(*T*) calculated for the T12 model using the set of 2322 UPOs is shown in figure 8*c*. Here, *Γ*(*T*+*Δ*) is the number of periodic points with periods inside the interval (*T*,*T*+*Δ*), *Γ*(*T*) the same quantity but for (*T*−*Δ*,*T*). The value of *Δ* was taken here as half of the inverse Kolmogorov entropy of the system (i.e. is approximately equal to 4.5 days.

In an important case of Axiom A systems (where one can guarantee the correctness of the UPO expansion procedure), the number of periodic points with period less than or equal to *T*(*Γ*(*T*)) is exponential with respect to the period and can be estimated as
where *H*_{top} is the topological entropy of the system [23]. As a result (i.e. if this theory is valid for our model) for the growth rate, we may expect (assuming that ). The line is shown as the upper horizontal line in figure 8*c* (the lower horizontal line is for *e*^{1/6} (i.e. the theoretical growth rate relaxed by a factor of 3). Of course, it is difficult to estimate true growth rate of periodic points with given period from figure 8*c* (the graph is very noisy), but we think it is two to three times lower than the theoretical estimate (with *H*_{top} estimated as ). As a result, we may expect some difficulties in using the UPO expansion theory in the case of our models.

The T21 model orbits are pictured in figure 9 ((*a*) all orbits, (*b*) 100 shortest orbits, (*c*) 20 least unstable orbits). Now, even using all the 1554 orbits found, we do not exactly match the attractor. Some of the least unstable orbits fall outside the trajectory attractor of the system again. Orbits are much more unstable compared with the T12 model case (average number of positive Lyapunov exponents is 27).

## 4. Unstable periodic orbits and model circulation regimes and statistics

In this section, we will describe the UPO expansion procedure. We will also present results on how, with the set of UPOs, one can reconstruct the system statistical characteristics and system circulation regimes for the T21 and T12 barotropic models.

When one wants to calculate the average value of some state-dependent quantity of interest, the standard strategy is to calculate a long enough system trajectory and evaluate the quantity at every trajectory point. Next, one has to calculate the average value of the obtained data. As a result, we arrive at the traditional formula for the average along the trajectory,
4.1
Here, *F* is a quantity of interest, 〈*F*〉 its average value, *S*(*n*,*ψ*_{0}) a system state at time *n*, *N* the number of measurements assumed to be large enough, *v*_{n} the weight of the measurement at time *n* (in most typical cases, all weights are equal to one) and *V* the total weight of the measurements (*V* =*N* in the previous case).

The UPOs of the system are embedded into the system attractor, and the system trajectory evolves in the phase space passing from one orbit to another. As a result, UPOs may approximate the system trajectory to some extent and could be used to approximate (4.1) as well (provided of course that many UPOs of the system are known). The idea of the approximation is to replace the averaging along the trajectory by the averaging along the orbits. From the physical point of view, it is clear that the orbits visited by a system trajectory more frequently must have larger weights in this approximation. Consequently, we obtain the approximation formula (4.2) that looks almost the same as (4.1),
4.2
Here, *ψ*^{i} is a periodic point of the system with period *T*′ less than or equal to some value *T* (i.e. *ψ*^{i}=*S*(*T*′,*ψ*^{i}); the UPO with period *T*′ has *T*′ periodic points; a stationary point of a system is a UPO with period 1), *I*(*T*) is the total number of periodic points of the system with periods less than or equal to *T*,*w*_{i} is the weight of a periodic point depending on orbit instability characteristics (all periodic points of the same UPO have the same weights). Note that only the orbits from the ‘trajectory attractor’ should be used in (4.2), and *I*(*T*) should not be small to obtain a decent approximation.

For the special class of the so-called Axiom A systems, this procedure could be theoretically justified [22,23]. For these systems, the set of periodic orbits is indeed dense on the attractor of the system, any arbitrary trajectory of the system can be approximated with any prescribed accuracy by the UPOs, and the equivalence of the trajectory averaging (4.1) and the UPO expansion (4.2) could be obtained [20,21,22,23]. As a result, all system statistical characteristics (including the probability to find the system trajectory near a given state) could be calculated using UPOs. According to this theory [22], weights *w*_{i} must be calculated as
4.3
where are positive Lyapunov exponents of the orbit containing the *i*th periodic point. The meaning of (4.3) consists of the fact that the orbit contribution into the average (4.2) should be inversely proportional to the volume expansion rate along the orbit unstable directions. As a result, the least unstable orbits make more contributions to the average.

System (2.1), studied here is not likely to be an Axiom A system. Following the ideas of Gallavotty [25], we will suppose that for typical multi-dimensional chaotic atmospheric systems, the UPO expansion may still work and apply the UPO expansion formulae (4.2) and (4.3) to the T12 and T21 models numerically.

First, we are going to check the accuracy of (4.2) and (4.3). As clear from figure 7, the number of UPOs does not grow with period as it should according to the theory, but it grows much slower. In the same time, it is clear intuitively that the weight formula should reflect somehow the UPO growth rate since the limit (4.2) must be non-zero. Consequently, one may expect that some modification of (4.3) may work better in practice.

First, we may relax (4.3) by some factor *k*_{r} accounting for the slower dependence of *Γ*(*T*) from period *T*,
4.4
Analysis of figure 8 shows that *k*_{r} could be somewhere from 2 to 4. Zoldi & Greenside [27] and Kazantsev [28] used another weight formula, namely,
4.5
In (4.5), the weight of a periodic point is proportional to the ‘exit time’ from its vicinity. Kazantsev [28,29] applied (4.2) and (4.5) to approximate statistical characteristics of the barotropic ocean model, and obtained very promising result. Another possibility to relax (4.2) was suggested in [41] as
4.6
where a value smoothed over some neighbourhood of the orbit is used instead of the exact value.

First, the numerical test of the UPO expansion method will be devoted to the reconstruction of the system PDF in the vicinity of the calculated UPOs. Let us calculate the probability for the system trajectory to visit the *ε*-neighbourhood of the *k*th orbit of the system. The orbit can be viewed as some set *Mk* in the system phase space. Define *Θ*_{A}(*ψ*) to be a characteristic function of the set *A*, being zero for *ψ* outside the set and one otherwise. Denote orbit *ε*-neighbourhood as *O*_{e}(*M*_{k}) and the probability for the system trajectory to be inside *O*_{e}(*M*_{k}) as *Ξ*_{Mk,ε}. *Ξ*_{Mk,ε} could be calculated directly (as we did in §2 of the paper) using the 1×10^{6} day trajectory. We just need to calculate the number of states from the dataset falling inside *O*_{e}(*M*_{k}). This could be formally rewritten using (4.1) with *F*=*F*(*ψ*)≡*Θ*_{A}(*ψ*) and *A*≡*O*_{δ}(*M*_{k}) as
4.7
Equation (4.7) simply means that in order to estimate probability *Ξ*_{Mk,ε}, we should find all points from ‘dataset1’ inside *O*_{e}(*M*_{k}) and every such point gives a contribution of *ν*_{n}/*V* ≡1/*N* into the *Ξ*_{Mk,ε} (*N*=1×10^{6} in our case).

On the other hand, according to the UPO expansion formalism, we can also estimate *Ξ*_{Mk,ε} using (4.2), i.e. as
4.8
In this case, in order to estimate probability *Ξ*_{Mk,ε} we should find all periodic points inside *O*_{e}(*M*_{k}) and each of these points gives a contribution of *w*_{i}/*W* into the *Ξ*_{Mk,ε}. For the T12 model, this calculation was done using all available periodic points with periods less than *T*=150 days (i.e. 2250 orbits were used in (4.8)). The value of *ε* was taken as 10^{−3} (the results are stable if *ε* is increased or decreased by factors up to 5). For the T21 model, all 1554 orbits were used in the calculations with *ε*^{2}=4×10^{−6}.

Using both the direct method (i.e. trajectory averaging) and the UPO expansion formula, we estimated probabiliies *Ξ*_{Mk,ε} and corresponding normalized probabilities for every UPO found for the system.

In figure 10*a*, the values of *η*_{k} calculated directly (*y*-axis, log10 scale) are shown against the ones obtained via the UPO expansion with theoretical weights (4.3) (*x*-axis, log10 scale) for each of 2322 orbits (i.e. *k*=1…2322). Ideally, we should expect *x*=*y* dependence. However, this is not the case, and the UPO expansion formula with weights (4.3) does not always correctly estimate *η*_{k}. In particular, this happens for the six slightly unstable orbits of the system (they are marked out in figure 10*a*). The true value of *η*_{k} for these orbits is several orders lower than its estimate with (4.3) and (4.8). Further analysis shows that these orbits are the ones outside the trajectory attractor (figure 7*b*). Consequently, these orbits should be removed form (4.8) to obtain correct representation of the system statistics. This helps a little, but the theory still underestimates *η*_{k} for very unstable orbits. It could be either due to the fact that the theoretic weight formula is incorrect for our system or simply because we missed some important orbits of the system. In both cases, the use of the approximation formulae (4.4)–(4.6) could help. Indeed, they give much better estimates for the distribution of *η*_{k} (figure 10*b* corresponds to the use of (4.6), figure 10*c* and *d* corespond to (4.5) and (4.4), respectively, and in figure 10*d*, we use *k*_{r}=4). The best results were obtained with (4.4) and (4.6).

The results of the same calculation but for the T21 model are presented in figure 10*e*–*h*. Now the best performance was obtained with (4.5) (figure 10*e*), probably because for T21, we have even less orbits so that (4.3) must be relaxed even more. Next in quality is (4.4) with large relaxation *k*_{r}=16 (figure 10*h*). Worst (but still qualitatively correct) are (4.4) with *k*_{r}=4 (figure 10*g*) and (4.6) (figure 10*f*), respectively.

The major conclusion that could be made from figure 10 is that the UPO expansion procedure is working, but for the best results, the theoretic weight formula (4.3) must be relaxed considerably.

The second test of the UPO expansion procedure is to reproduce the PDF of the system projected onto its leading EOFs (figure 3*a*–*d*). The numerical method is analogue to the previous case and is based on the box-counting algorithm. We need to approximate the probability of the system to be inside two-dimensional boxes of the selected EOF plane using the UPO expansion. This can be done with the help of (4.8), where *Θ*_{A}(*ψ*_{0}) is now a characteristic function of the selected two-dimensional box (and in (4.8), we are counting all periodic points with the appropriate two-dimensional coordinates). The calculation with the best weight function from the previous experiment (i.e. (4.6)) was performed, and the results were shown in figure 3*e*–*h*. Comparing corresponding pictures, one can see that the UPO-based approximations of the PDF projections have reasonable quality. Not only the shape of the PDF is reproduced, but also its orientation, ridges and local features. The difference between corresponding distributions is never more than 10 per cent in the L2 norm.

In addition, we repeat the same calculation but for the vorticity PDF projected onto the planes of its leading (vorticity) EOFs. PDF reconstructions (figure 11*e*–*h*) reflect most of the features of the true PDF projections (figure 11*a*–*d*).

Other basic statistical characteristics of the T12 model, including its average state, standard deviation, leading EOFs and lagged covariances were also reconstructed with the UPO expansion procedure (not shown). They are almost the same as the ones calculated using the model output, which it is expected taking into account the good quality of the reconstructed PDF.

In the T21 model case, we cannot expect the same results from (4.8) because of its much worse performance in the previous experiment. The reconstructions of corresponding PDF projections with (4.5) and (4.8) (that gave the best results in the first experiment) were presented in figure 4*e*–*h*. Though the positions of the PDF maximum are always correct and PDF orientation is almost correct, differences with figure 4*a*–*d* are apparent. The reconstructed PDFs have much sharper maxima with almost zero tails. The most evident reason is that we missed a huge number of orbits on the sides of the attractor.

For the T21 model, we also calculated its mean state (figure 12*a*), standard deviation (figure 12*b*) and leading EOFs. The mean state is reconstructed almost exactly (cf. figures 1*b* and 12*a*); the standard deviation has correct structure and somewhat low amplitude (figures 1*f* and 12*b*). The variability patterns are different, mainly because we were unable to calculate enough orbits with periods in a 14–30 day range (only 300 long orbits were found). As a result, short orbits dominate variability, giving incorrect small-scale high-frequency patterns as leading EOFs of the system. Leading EOFs of the low-passed orbits are much closer to reality.

Now let us discuss the relationship between MPSs and MPT of the system and UPOs. The situation could be best viewed in the example of the T12 system where we know enough long-period orbits. In figure 13*a*–*c*, the dynamics of MPT1 is shown in the EOF1–EOF2 plane (figure 13*a*), the EOF1–EOF3 plane (figure 13*b*) and the EOF2–EOF3 plane (figure 13*c*) as black curves. The black circle marks placed over the trajectory indicates system position at every 4 days of its evolution. Projections of 15 leading MPSs are also shown in figure 13*a*–*c* as black dots. MPS1 is marked additionally by a circle; MPS2, square; MPS3, diamond; MPS4, triangle; MPS5, circle with bar. Suggested MPT1 evolution from MPS4 (at about day 4) to MPS1 (at day 10) and finally to MPS5 (at day 20) is clearly seen here.

Next, we calculated the distances between MPT1 and all UPOs. It was found that MPT1 is situated in the region filled by several least unstable orbits (orbits NN 5, 9, 18, 21 are the closest ones, all having MPT1 in their *ε*-neighbourhoods with *ε*^{2}=1×10^{−6}). Part of the UPO5 reflecting the motion of MPT1 is presented in figure 13*a*–*c* as a grey curve with open circles, whereas UPO21 is shown as a grey curve with closed circles. The second least unstable stationary point is also shown here as a cross. Obviously, the UPO5 (and UPO21) follow closely MPT1 and all evolve around a stationary point on the EOF1–EOF2 and EOF2–EOF3 planes. The nearest MPS to the stationary point is MPS2. Evolutions of MPT1 and UPO21 in physical space are shown in figure 6 (MPT1 states at days 0, 4, 10, 15, 20, 25 are shown in figure 6*a*,*b*,*c*,*d*,*i*,*j*, respectively). As expected, UPO21 (figure 6*e*–*h*,*j*,*k*) closely follows MPT1. As a result, we can conclude that least unstable UPOs approximate most probable system trajectories. MPSs in turns are situated in the vicinity of these orbits.

For the T21 model, we can also find some connection between the least unstable orbits and model regimes. As an example in figure 14, we present the evolution of MPT1 of the T21 system (at days 0, 4, 8 and 12) and the closest orbit (19th least unstable orbit of the model to be specific). One can also see the orbit tracing MPT1, but it is not as close as for T12. One should point out however that the attractor dimension of the T21 system is five times larger than that of T12, and many more orbits are required for comparable accuracy.

The final question studied in this paper is the reproduction of the probability for the system to visit the vicinity of MPTs. The only difference from the first test (where we calculated the probability for the system to visit the UPO neighbourhood) is that we should calculate probabilities for any of the 10 000 MPTs found in §2 instead of the set of UPOs as in the first test. The result of this experiment is shown in figure 15*a* (T12 model) and figure 15*b* (T21 model). Here, for the T12 and T21 models, we used weight function formulae (4.6) and (4.5). In the T12 case, UPOs give reasonably close estimates for the regime probabilities (note that true regimes are of course the ones with largest probabilities) with the error being less than 50 per cent in most cases. The somewhat wrong slope of the probability ratios could be due to the fact that we still do not know some orbits in unpopulated areas of the T12 attractor. For the T21 model, UPOs give systematically higher values for the regime probabilities (by a factor of 6 on average). As mentioned earlier, this is mostly because of the fact that many UPOs are missed by numerical methods. Relative probabilities are more or less correct though.

## 5. Conclusion

The theory of chaotic dynamical systems gives many tools that can be used in climate studies. The widely used ones are the Lyapunov exponents, the Kolmogorov and information entropies, the attractor dimension, etc. Most of them are global quantities describing the system as a whole. But the local properties of the system are at least as important. Some regions of the phase space could be more populated than the others, resulting in PDF local maxima. Complexity of the system behaviour and its predictability could change from one region of the phase space to another, as well as its sensitivity to the external forcings. Fortunately, the theory gives a useful tool to study the local properties of the system attractor. It is based on the system trajectory approximation with the help of its UPOs. The system PDF and its statistical characteristics could be obtained throughout the UPO expansion as a weighted sum over the orbits. Orbit weight is inversely proportional to the orbit instability characteristics, so less unstable orbits make larger contributions to the results. For the special classes of chaotic dynamical systems, there is a strict theory guaranteeing the accuracy of this approach. The major disadvantage of the proposed technique is its numerical complexity (the search of UPOs is very expensive numerically, and for large systems, one has to calculate many orbits). In addition, typical atmospheric systems do not have the properties required by the theorems beneath the UPO expansion method so numerical verification of the idea is required.

In our paper, we tried to use UPOs for approximating statistical characteristics of the relatively simple barotropic model of the atmosphere. For the T12 resolution system, we were able to calculate more than 2000 UPOs, which was enough to reproduce most system characteristics. The UPO-based approximations of the system PDF have a reasonable quality. Not only the shape of the PDF is reproduced, but also its orientation, ridges and local features. The difference between corresponding distributions is never more than 10 per cent in the L2 norm. The same results were obtained for the vorticity PDF of the model. As a consequence, basic statistical characteristics of the model, including its average state, standard deviation, leading EOFs and lagged covariances were automatically approximated by the UPO expansion procedure with excellent accuracy. It should be noted however that the theoretic weight formula must be relaxed considerably for better performance of the UPO expansion procedure. It was also found that some of the least unstable orbits of the system make much smaller contributions to the system measures than they were supposed to (taken into account the large weights), and do not define any notable MPS (i.e. they are invisible for the observer). One may speculate that these orbits belong to the so-called invisible part of the attractor [42].

Another useful consequence of the UPO expansion is the fact that the probability to see the system near a given state of the phase space has increased values in the vicinity of the least unstable orbits. This is why one can expect to find system MPSs near the least unstable UPOs. The leading MPS calculated for the T12 model is close to the regime R obtained in [17]. MPS4 and MPS5 could be also indentified as regimes A and G, respectively, from [17]. The most probable system trajectory connects MPS4 with MPS1 and MPS5 in a way that it takes about 6 days for the system to reach MPS1 (regime R) from MPS4 (regime A) and another 10 days to get to MPS5 (regime G). This is close to the results of [17], who estimated most probable A–R and R–G transition times to be 6 days and suggested that the transition from A to G should pass through the R regime. Further analysis showed that several least unstable UPOs do in fact closely trace this MPT and connecting regimes A–G and R, as was suggested by [17]. This connection is provided by 20 day pieces of 100 day orbits rather than a 20 day orbit, as was suggested by [17]. It should be pointed out, however, that correspondence of the results with the [17] model may be somewhat fortuitous (the models are different and both are simple enough).

In the T21 case, 1554 orbits that were found are clearly not enough to obtain all the information about the system. Nevertheless, the trajectory prefers to stay near the least unstable orbits of the model (however, invisible orbits do exist). Again, as in the T12 model case, there is a connection between the least unstable orbits and model MPT segments. MPT1 of the T21 system passes through the region densely populated by the least unstable UPOs and is traced by one of them. Basic statistical characteristics of the system could also be recovered via the UPO expansion procedure. As for the fine properties of the system PDF and regime behaviour, they could not be reconstructed correctly due to the lack of UPOs found for the system.

Taking into account all our results, we think that the UPO expansion strategy could be a very useful (though numerically expensive) tool for a broad range of problems.

## Acknowledgements

This research was supported by RFBR project 11-05-01119 and the Ministry of Education and Science of Russian Federation projects 8326 and 8671. Two anonymous reviewers provided helpful suggestions and comments to improve the paper significantly. The author is grateful to G. Branstator for very useful discussions.

## Appendix A

Let us briefly describe the numerical method used to calculate UPOs in our work. Rewrite (3.1) as follows:
This is the system of *N* nonlinear equations (*N* is the dimension of the model phase space) with respect to *N*+1 unknowns *ψ*_{0} (the *N*-dimensional vector of the orbit’s initial position) and *T* (the orbit period). The ambiguity of the system is caused by the fact that the same orbit can be determined by different initial conditions obtained from each other by a shift along a periodic trajectory.

Let *x*^{i} and *T*^{i} be the *i*th approximations for *ψ*_{0} and *T*. Denote *x*^{i+1}=*x*^{i}+*h*^{i}, *T*^{i+1}=*T*^{i}+*τ*^{i} and *x*^{i}(*T*^{i})=*S*(*T*^{i},*x*^{i}). The standard iteration of the Newton method is based on the expansion of the equation
into a Taylor series with respect to *h*^{i} and *τ*^{i} up to the first-order terms
A 1
With ∂*F*(*x*^{i},*T*^{i})/∂*x*^{i}≡∂(*S*(*T*^{i},*x*^{i})−*x*^{i})/∂*x*^{i}=*L*(*T*^{i},*x*^{i})−*E* and ∂*F*(*x*^{i},*T*^{i})/∂*T*^{i}≡∂(*S*(*T*^{i},*x*^{i})−*x*^{i})/∂*T*^{i}=*f*(*S*(*T*^{i},*x*^{i})), we obtain a linear system of *N* equations for *N*+1 unknowns *h*^{i} and *τ*^{i},
A 2
To remove the uncertainty, one has to complement the system by the so-called phase condition. We used the simplest one by requiring the orthogonality between the search direction and system trajectory direction (i.e. system (2.1) right-hand side) at current state as
A 3
From (A2) and (A3), one can determine the search direction (*h*^{i},*τ*^{i}) and calculate the next iteration as *x*^{i+1}=*x*^{i}+*ηh*^{i}, *T*^{i+1}=*T*^{i}+*ητ*^{i}. The step damping factor *η* equals the one for the standard Newton method, but for the ill-conditioned linear problem, the step could become large, making the first-order Taylor expansion incorrect. In our work, we used the ‘line search procedure’ [43] for *η* all the time to improve the convergence of the algorithm. For numerical implementation of the Newton method, we have to calculate the operator of the tangent linear system ∂(*S*(*T*^{i},*x*^{i}))/∂*x*^{i}=*L*(*T*^{i},*x*^{i}). As a result, one iteration of the Newton method costs approximately the same as a run of the full tangent linear system for the current value of *T*^{i}. When systems (A2) and (A3) become very ill-conditioned, we use a second-order tensor correction to the Taylor expansion (A1). Resulting changes in the numerical algorithm were implemented according to [44,45], and are described in [35].

The convergence of the Newton algorithm crucially depends on the quality of the initial conditions. In our work, we used two strategies. First, we scanned a long system trajectory to find a pair of close points and used the first point from the pair as an initial guess for *ψ*_{0} and the time in between as *T*_{0}. Second, we randomly took the initial position at the model attractor and made initial guesses for the period of the orbit from the uniform distribution *k*=1…*K*_{0} for each of those points. With the first approach, the method tends to find orbits in most visited parts of the attractor. With the second one, we were able to obtain the orbits far from the attractor core.

## Footnotes

One contribution of 13 to a Theme Issue ‘Mathematics applied to the climate system’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.