The goal of this paper is to provide a complete picture of the long-term behavior of Lorenz’s maximum simplification equations along with the corresponding meteorological interpretation for all initial conditions and all values of the parameter.
In this paper our goal is to analyze the properties of equilibrium points and long-term behavior of the solution of the Lorenz’s (1960) maximum simplification equations (MSE). Lorenz (1960) introduced the maximum simplification equations and demonstrated that they could be used to understand certain atmospheric phenomena, such as the index cycle. The dynamical equations were simplified to a point where reasonable predictions of atmospheric conditions could not be expected, however, certain physical processes and features of interest were retained. This MSE consists of a set of three nonlinear and coupled system of three ordinary differential equations (ODE) that govern the evolution of the amplitudes of three modes arising from the double Fourier series expansion of the solution of the barotropic vorticity equation (BVE). The importance of MSE stems from the fact that while it is the simplest possible (minimal set) approximation to BVE, using it, Lorenz (1960) was able to capture many of the qualitative features of the flow including the exchange of kinetic energy between the zonal and meridional components. Motivated by this success, the basic spirit of this ground-breaking work of Lorenz (1960) namely approximating the solution of complex realistic problems using a low-order model (LOM) has become a standard arsenal in problem solving in geophysical domain. Refer to Saltzman (1962), Lorenz (1963), De Swart (1989), Wiin-Nielson (1992), Gluhovsky et al. (2002), and the references therein for details.
Thompson (1957) used a set of equations similar to MSE to explore the limits of deterministic predictability that stem from the uncertainty in the initial conditions. Epstein (1969) used MSE to study the stochastic dynamic prediction. Recently, Lakshmivarahan et al. (2003) used MSE to compare the performance of first-order and second-order methods for data assimilation. To our knowledge, long-term behavior of Lorenz’s MSE has not been explored. In this paper we examine the growth of small errors superimposed on equilibrium points of MSE leading to a set of interesting heteroclinic orbits. The necessary conditions for obtaining a solution that approaches an equilibrium point are established. In addition, we provide a complete catalog of the behavior of MSE starting from any arbitrary initial point. The impacts of changing various model parameters on the solutions with respect to the index cycle are also examined.
In section 2, MSE and their key properties are derived. The set of all equilibria for MSE are identified and a set of necessary conditions for the orbits to approach an equilibrium point are derived in section 3. A comprehensive analysis of stability of equilibria is covered in section 4. A complete picture of the evolution of small errors superimposed on equilibria is given in section 5. A complete classification of the behavior of the solution starting from an arbitrary point is covered in section 6. Meteorological implications of the above analysis are given in section 7. The concluding section 8 provides a summary.
2. The maximum simplification equations—The model
Our model equations are the Lorenz’s so-called maximum simplification equations given by
where A = A(t), F = F(t), and G = G(t) are the time varying amplitudes of the three spectral components obtained from the standard BVE and Ȧ, Ḟ, and Ġ refer to their time derivatives. See Lorenz (1960) for details of the derivation. The streamfunction ψ = ψ(x, y, t) and the vorticity ∇2ψ = ∇2ψ(x, y, t) corresponding to (2.1) are given by
where x and y are the usual space coordinates, ∇2 is the standard Laplacian operator, 2π/l = L is the distance between zonal wind maximum, and 2π/k = K is the wavelength of the disturbance. That is, K and L denote the spatial scales of the flow and the disturbance. As an example, L = 10 000 km and K = 5000 km, that is, l = 6.285 × 10−4 and k = 1.257 × 10−3. The expression for the easterly wind u = u(x, y, t) = −∂ψ/∂y and northerly wind υ = υ(x, y, t) = ∂ψ/∂x are given by
a. Dependence on the parameter r
It can be verified that the model (2.1) depends only on the ratio r = k/l = L/K and can be written as
A plot of the values of a = a(r), b = b(r), and c = c(r) as a function of r > 0 is given in Fig. 1. It can be verified that b = c when r = 41/3 = 0.7598 and a = c when r = 43 = 1.3161. Further, b → 0+ (from above) as r → 0 and a → 0− (from below) as r increases.
Let η = (A, F, G)T ∈ ℝ3 be a column (real) vector of three spectral components, where T denotes the transpose. That is, η1 = A, η2 = F, and η3 = G. Then (2.6) can be succinctly written as
where g: ℝ3 → ℝ3 is such that g(η) = [g1(η), g2(η), g3(η)]T with g1(η) = aFG, g2(η) = bAG, and g3(η) = cAF.
It should be interesting to note that the Eq. (2.6) is equivalent to the Euler’s equation of motion for a rigid body and is known as the Euler gyrostat (Obukhov 1973). Derivation and analysis of Euler gyrostat can be found in textbooks on classical mechanics (Landau and Lifshitz 1960; Wittenburg 1977; Arnol’d 1989). Recently Gluhovsky and his collaborators (Gluhovsky and Agee 1997; Gluhovsky and Tong 1999; Gluhovsky et al. 2002) have shown that a Volterra gyrostat (of which the Euler gyrostat is a special case) can be used as a building block for generating a wide variety of low-order models of interest in several application domains.
Implicit in the above description are the two sets of coordinate systems. First is the three-dimensional physical space consisting of the two space variables, x and y, and the time variable t, in which the physical quantities such as the streamfunction ψ, vorticity ∇2ψ, and the wind components u and υ are described. The second is the three-dimensional spectral space in which the spectral amplitudes A, F, and G evolve as functions of time according to (2.6) or (2.10). The goal is to characterize the long-term behavior of the physical quantities from those of the spectral components.
b. Special case of r = 1
It can be verified that a = −½, b = ½ and c = 0 when r = 1. In this case, Ġ = 0 and G(t) = G(0) for all t. Hence, if G(0) ≠ 0, (2.6) reduces to a linear harmonic oscillator
whose solution is periodic (Hirsch and Smale 1974). When G(0) = 0 and r = 1, then η(t) ≡ 0. In the following it is tacitly assumed that r ≠ 1.
c. Conservation of kinetic energy and vorticity
Let 𝗗 = diag(d1, d2, d3) be a diagonal matrix where di > 0, i = 1, 2, and 3. Consider the functional f: ℝ3 → ℝ defined by
It can be verified that the time derivative of this functional evaluated along the trajectory of (2.6) is given by
Thus, given any parameter vector p = (a, b, c)T of the model, for any d = (d1, d2, d3)T that is orthogonal to p, ḟ(η) ≡ 0 and for this choice of the vector d or matrix 𝗗, f (η) is conserved by the model. By direct computation it can be verified that the mean kinetic energy
and the mean square vorticity
d. Solution space and degree of freedom
Let η(0) = [A(0), F(0), G(0)]T be the given initial condition for (2.10). Then (2.15) implies that the solution η(t) of (2.10) is confined to the surface of the V ellipsoid, which is centered at the origin and is given by
Clearly, fV(η) depends only on η(0) and not on r. Similarly, from (2.14) it follows that the solution η(t) of (2.10) is confined to the surface of the E ellipsoid, which is also centered at the origin, and is given by
Combining these, it follows that the solution η(t) of (2.10) must lie in the set
where ∇fV(η) = (1/V(0))(A, F, 2G)T is the gradient of fV(η) and
where ∇fE(η) = [1/2E(0)][(A/l2), (F/k2), (2G/l2 + k2)]T is the gradient of fE(η). It can be verified that the gradient vectors ∇fV and ∇fE are linearly independent unless r = 1 and G = 0. Thus, the flow of the dynamics in (2.10) is always orthogonal to the plane defined by these two gradient vectors. Since ∇fV and ∇fE define the outward normal to the surface of the ellipsoids, the flow always lies along the tangent to the points in 𝒮 in (2.20).
Lorenz (1960) has shown that (2.6) can be solved analytically where the solution A(t), F(t), and G(t) are expressed as elliptic functions of time. In principle, we can use the tables of elliptic functions (Abramowitz and Stegun 1965) to characterize the long-term behavior of the solution of (2.6). However, following Lorenz (1960), we use the standard numerical methods to solve (2.6). While Lorenz (1960) used the central difference scheme that involves three time steps, we use the fourth-order adaptive Runge–Kutta (RK) method.
e. Boundedness of the trajectories
An immediate consequence of the property described in section 2b is that the trajectories η(t) of (2.10) remains bounded. Let η(0) = [A(0), F(0), G(0)]T be the initial state. Then from (2.14) and (2.15) we get
For example, when A(0) = 10−3 = G(0) and F(0) = 5, K = 5000 and L = 10 000, it can be verified that |A(t)| ≤ 2.5, |F(t)| ≤ 5.0, and |G(t)| ≤ 3.5355, for all t ≥ 0.
This property in turn implies that the model (2.1) is (externally) stable in the sense of Lyapunov (Hirsch and Smale 1974; Sparrow 1982). However, in view of the nonlinear interaction between modes, there is transfer of energy between modes leading to (internal) instability.
An example illustrating this transfer of energy between modes is given in Fig. 3 where the system (2.6) is solved using the fourth-order adaptive Runge–Kutta method. This figure provides a complete picture of the model trajectories, namely the variation of A(t), F(t), and G(t) with time; two dimensional plots of A(t) − F(t), A(t) − G(t) and F(t) − G(t), and the three-dimensional view of the solution.
3. Equilibrium points of the model
Recall that if = g(η) represents a dynamical system, then the set of points at which g(η) = 0 are called the equilibrium points of the system (Hirsch and Smale 1974; Devaney 1989; Sparrow 1982). Accordingly, we define three sets of equilibria for the system in (2.6). Let
It can be verified that ɛA, ɛF and ɛG constitute the set of all equilibrium points of (2.6). It can be seen that the origin is the only point common to all these three sets.
Given that the solution of the model (2.10) always lies in the set 𝒮 in (2.20), a question of interest at this point is: Which of the above equilibrium points lie arbitrarily close to the solution space 𝒮? Stated in other words, when does the model solution starting from an arbitrary initial point ever approach an equilibrium point?
In answering this question, first consider η* = (A*, 0, 0)T ∈ ɛA. It can be verified that η* ∈ 𝒮𝒜 when
that is, when
in which case A* = ±2V = ±2lE. Thus, (3.2) is a necessary condition for an orbit to include η* ∈ ɛA.
Considering η* = (0, F*, 0)T ∈ ɛF and repeating the above argument, it can be verified that
is a necessary condition for an orbit to include η* ∈ ɛF where F* = ±2V = ±2kE. Similarly, considering η* = (0, 0, G*)T ∈ ɛG, it readily follows that
is necessary for the orbit to include η* ∈ ɛG, where G* = ±2V = ±2(l2 + k2)E.
Stated in other words, the key to understanding the behavior of orbits of the model in (2.10) is the value of the ratio V/E. Lorenz (1960) recognized the importance of V/E (= 2k2 as the threshold for sub and supercritical disturbances, for example), but did not explain why it was significant. We see here that the necessary conditions for obtaining any solution including an equilibrium point as part of the orbit are set by specific values of the ratio V/E.
4. Stability of equilibria
Recall that the local stability properties of an equilibrium point, say η*, is governed by the distribution of the eigenvalues of the Jacobian, 𝗗g(η) of the mapping g(η) in (2.10) evaluated at η* (Hirsch and Smale 1974; Devaney 1989; Lewis et al. 2006, chapter 32). The Jacobian of g(η) is given by
a. Stability of equilibrium points in ɛA
Let η* ∈ ɛA. Then, the Jacobian (4.1) becomes
Let r < 1. Then, a < 0, b > 0, and c > 0. It can be verified that the eigenvalues of 𝗗g(η*) are given by
and the corresponding eigenvectors are given by
For latter reference, define the diagonal matrix Λ(1) of eigenvalues and 𝗩(1), the nonsingular matrix of the corresponding eigenvectors as
where (it follows from the definition)
Notice that while the eigenvalues λ(1)2 and λ(1)3 in (4.3) depend on A and the ratio r, the eigenvectors v(1)2 and v(1)3 are independent of A. That is, every point in ɛA shares the same set of eigenvectors that depends only on the value of the ratio r. For purposes of comparison, the eigenvalues and vectors of 𝗗g(η*) when η* belongs to each of the three equilibria are summarized in Table 1.
Let r > 1. Then a < 0, b > 0, and c < 0. It can be verified that the complex eigenvalues of 𝗗g(η*) are given by
and the corresponding complex eigenvectors are v(1)1, w(1)2 = u1 + iu2 and w(1)3 = u1 − iu2, respectively, where
We leave it to the reader to verify that this eigen decomposition in the complex domain can be equivalently expressed in real domain using the standard canonical form (Hirsch and Smale 1974) as
That is, the eigenvectors that are the columns of 𝗩(1) are given in (4.8). In this paper we use this latter form of eigen decomposition when the eigenvalues are complex. Refer to Table 1 for a summary.
Variation of bc and c/b as a function of r are given in Figs. 4a,b, respectively. It follows from (4.3)–(4.4) and (4.8)–(4.9) that the eigenvalues and eigenvectors depend on two factors: the value of the first coordinate A of the equilibrium point η* = (A, 0, 0)T and the value of the ratio r. For completeness, we consider several subcases.
1) Case 1: Let A > 0 and r < 1
In this case, it can be verified that λ(1)2 > 0 and λ(1)3 < 0, and the equilibrium point, η* = (A, 0, 0) is called a saddle point, with the eigenvector v(1)2 corresponding to the unstable (repelling) direction and the eigenvector v(1)3 to the stable (attracting) direction (Hirsch and Smale 1974). Refer to Table 2 for a summary. In verifying these claims, let η = η* + ξ where ξ is a small perturbation. Then, using the first-order Taylor series expansion, we obtain
as the local dynamics of evolution of the perturbation. Now change the variables from the standard (A, F, G) orthogonal coordinate system to oblique coordinate system defined by the eigenvectors v(1)1, v(1)2 and v(1)3 in (4.6). Define a new vector z = (z1, z2, z3)T where
and z2(t) = z2(0)eλ(1)t2 and z3(t) = z3(0)eλ(1)t3. Since λ(1)2 > 0, indeed z2 grows exponentially and λ(1)3 < 0 implies that z3 decreases exponentially as well. Refer to Fig. 5 for a geometric exposition of these ideas. Thus, any perturbation with a nonzero component along v(1)2 will grow. Hence, in this case η* = (A, 0, 0) is an unstable equilibrium point. Referring to Fig. 5, the angle θ that v(1)2 makes with the positive F axis (same as the angle v(1)3 makes with the negative F axis) is given by
Thus, as r increases from 0 to 1, θ decreases from 90° to 0°.
2) Case 2: Let A < 0 and r < 1
In this case since λ(1)2 < 0 and λ(1)3 > 0, the role of the eigenvalue–eigenvector pairs (λ(1)2, v(1)2) and (λ(1)3, v(1)3) are interchanged and the above analysis and conclusions hold. Thus, the equilibrium point η* = (A, 0, 0) is an unstable saddle point. Refer to Fig. 5.
(i) Effect of changing A
Thus, keeping r(<1) fixed, the force of attraction along v(1)3 and the force of repulsion along v(1)2 as measured by the magnitudes of the eigenvalues |λ(1)3| = Abc and |λ(1)2| = Abc, respectively, diminish to zero as A decreases from 1 to 0. As A increases in value in the negative direction, while the magnitudes of these forces increase, their role is reversed, with v(1)2 being the new direction of attraction and v(1)3 being the new direction of repulsion. This reversal of the role between v(1)2 and v(1)3 occurs when A = 0.
(ii) Effect of changing r
As r increases form 0+, the angle θ decreases from 90°. When r → 1−, (c/b) → 0+ and θ → 0+. Again we get a degenerate limit
When r = 1, c = 0 and the character of the equilibrium η* = (A, 0, 0)T changes from being a saddle point to a center. Refer to the next subsection for more details.
3) Case 3: Let A > 0 and r > 1 or A < 0 and r > 1
From (4.8) it follows that the nonzero eigenvalues are complex conjugate with zero real part. In this case the equilibrium point is called a center for the corresponding linearized problem (Hirsch and Smale 1974). A geometric view of the two new eigenvectors in (4.9) is shown in Fig. 6. It turns out that while the new eigenvector v(1)3 = (0, 1, 1)T is independent of r and makes a fixed 45° angle with the positive F axis, the new eigenvector v(1)2 = (0, −|(b/c)|, |(c/b)|)T makes an angle θ with the negative F axis where
As r increases from one to infinity, |c/b| increases form 0 to ½ and θ from 0 to 26.57°.
(i) Change of character of the equilibrium η* = (A, 0, 0)T
As r changes from 1− to 1+, it follows from (4.17) that the old eigenvectors v(1)2 and v(1)3 in (4.4) become degenerate and cease to exist. Simultaneously, a new set of eigenvectors v(1)2 and v(1)3 in (4.9) are created. Of these, v(1)3 = (0, 1, 1)T is fixed. When r → 1+, |c/b| → 0 and hence new v(1)2 → (0, −1, 0)T. Comparing this with (4.17), it looks as though the old v(1)3 is reincarnated as the new v(1)2 as r changes from 1− to 1+. Thus, as r changes from 1− to 1+, η* = (A, 0, 0)T changes its character from being a saddle to a center.
4) Case 4: Either A = 0 or r = 1
In this case, λ(1)2 = 0 = λ(1)3.
b. Stability of equilibrium points in ɛF
Let η* = (0, F, 0)T ∈ ɛF. Then, the Jacobian in (4.1) becomes
whose eigenvalues and eigenvectors are given by
with λ(2)1 = 0, λ(2)2 = Fac, λ(2)3 = −Fac and
The product ac and the ratio c/a are given by
Variation of ac and c/a as a function of r are given in Figs. 7a,b, respectively. From (4.20) and (4.22) it follows that the eigenvalues λ(2)2 and λ(2)3 depend on two factors: the value of the second coordinate F of the equilibrium point η* = (0, F, 0)T and the value of the ratio r. A summary of the properties of these eigenvalues for various combinations of values of F and r are given in Table 3.
Again it can be verified that η* = (0, F, 0)T is an unstable saddle point when r > 1 and F ≠ 0. However, when r < 1 and F ≠ 0, perturbations around this equilibrium will have a local periodic behavior. It is interesting to note that the behavior of the equilibria in ɛA for r < 1 are quite similar to those in ɛF for r > 1 and vice versa.
c. Stability of equilibrium points in ɛG
Let η* = (0, 0, G)T ∈ ɛG. Then the Jacobian in (4.1) becomes
whose complex eigenvalues are given by λ(3)1 = 0, λ(3)2 = iG|ab|, λ(3)3 = −iG|ab|. Using the standard canonical form, the eigen decomposition becomes
The product ab and ratio b/a are given by
From (4.25) it follows that the eigenvalues λ(3)2 and λ(3)3 are always complex conjugates (with zero real parts) for all values of r. Variation of ab and b/a as a function of r are given in Figs. 8a,b, respectively. Refer to Table 4 for a summary of the properties of these eigenvalues.
5. Small perturbation of equilibrium
In this section our goal is to quantify the growth of small errors superimposed on equilibrium points.
Let η* be an equilibrium point and ξ ∈ ℝ3 denote a (small) perturbation or error vector. Let η = η* + ξ be the perturbed initial condition for the model (2.6). Since * = g(η*) = 0 by definition, we get
as the exact dynamics of the perturbation. Following the developments in section 4, we consider three cases, one for each of the three sets of equilibria.
a. Small perturbation of equilibrium points in ɛA
Let η* = (A, 0, 0)T ∈ ɛA. Depending on the values of A and r, consistent with section 4, we again consider four subcases.
1) Case 1: Heteroclinic orbits
Let A > 0 and r < 1. In this case, η* = (A, 0, 0)T is a saddle point with v(1)2 = (0, 1, c/b)T denoting the local unstable direction and v(1)3 = (0, 1, −c/b)T denoting the local stable direction. Refer to Tables 1 and 2 for details.
These eigenvectors v(1)2 and v(1)3 partition the F–G plane (where A is constant) into four regions R1 to R4. Refer to Fig. 9 for details. It turns out that there is an intrinsic difference in the qualitative behavior of the trajectories starting from each of these four regions but close to the equilibrium point η*.
For concreteness, let A = 1.0, r = 0.5, and ξ(0) = (0.0, 0.001, 0.0)T ∈ ℛ1. Then, η* = (1, 0, 0)T and η = (1.0, 0.001, 0.0)T. A complete picture of the evolution of the model Eq. (5.1) starting from η is given in Fig. 10. A number of observations are in order.
(i) Heteroclinic orbits
The orbit or the trajectory starting from η ∈ ℛ1 that is close to the equilibrium point η* = (1, 0, 0)T moves away from it along the unstable direction v(1)2. In due course, this orbit bends with A decreasing but F remaining positive. It soon approaches the equilibrium point η** = (−1, 0, 0)T, which is the symmetric dual of η*, along v(1)2, which is the local direction of attraction for η**. As this trajectory gets very close to η**, it gets repelled along v(1)3, which is the local unstable direction for η**. Again it bends with A increasing but F still remaining positive. It now approaches η* along v(1)3, which is its direction of attraction. As it gets closer to η*, it again gets repelled along v(1)2 and this cycle repeats. Figure 10 shows the results of integration over a time period T = 200 units using the adaptive fourth-order Runge–Kutta scheme.
The class of orbits so generated is called a heteroclinic orbit (Sparrow 1982). A distinguishing feature of this orbit is that one part of it connecting η* and η** lies at the intersection of unstable manifold of η* and stable manifold of η**. The other part of the orbit, from η** to η*, lies at the intersection of the unstable manifold of η** and the stable manifold of η*.
(ii) Four distinct heteroclinic orbits
It turns out that each equilibrium point η* ∈ ɛA is endowed with four distinct classes heteroclinic orbits—one for each of the four regions ℛ1 to ℛ4. Figures 11 –13 contain illustrations of heteroclinic orbits starting from initial points in regions ℛ2 to ℛ4, respectively. The key distinguishing features of these four heteroclinic orbits are summarized in Table 5.
(iii) Effect of changing A
From Tables 1 and 2, it is clear that eigenvectors v(1)2 and v(1)3 of 𝗗g(η*) for η* ∈ ɛA are independent of A. But their associated eigenvalues λ(1)2 = Abc and λ(1)3 = −Abc, which are measures of the force of repulsion and attraction are directly proportional to A. Thus, as A increases in magnitude, the resulting increases in the force of repulsion and attraction, while do not alter the basic properties of the heteroclinic orbits, do affect the frequency of oscillations of the components of η. Examples of heteroclinic orbits for A = 5 and 10 with ξ = (0.0, 0.001, 0.0)T and r = 0.5 are given in Figs. 14 and 15.
(iv) Effect of changing r
From Fig. 4, it follows that c/b is a decreasing function of r with c/b being positive for 0 < r < 1 and negative for r > 1. We restrict this discussion to the case when 0 < r < 1 (r > 1 is considered in case 3 below). Referring to Fig. 9, it follows that (c/b) determines the slope of the eigenvectors v(1)2 and v(1)3 in the F–G plane. To understand the impact r, consider the example of the heteroclinic orbit in Fig. 10. In this example, since F remains positive, combining it Fig. 9, it can be verified that the heteroclinic orbit approaches η* = (1, 0, 0)T along v(1)3 (that separates region ℛ1 and ℛ4) and leaves of η* along v(1)2 (that separates region ℛ1 and ℛ2). That is, this orbit makes a turn (to the right) by an angle ϕ = 180 − 2θ where θ is given by (4.15). Clearly, this angle ϕ is small or large depending on if r is small or large.
(v) Effect of changing L and K
2) Case 2: A < 0 and r < 1
In this case, the equilibrium point η* = (−A, 0, 0)T ∈ ɛA is again a saddle point, with v(1)2 being the stable direction and v(1)3 the unstable direction. Hence all the comments and conclusions of case 1 carry over to this case.
3) Case 3: A > 0 and r > 1 or A < 0 and r > 1
In this case b > 0 but c < 0 and recall from Table 1 that the eigenvalues λ(1)2 = iA|bc| and λ(1)3 = −iA|bc|. This in turn implies that η* = (A, 0, 0)T is a center and the components F(t) and G(t) oscillate where the frequency of oscillation depends on both A and r. In verifying this, assume that A is nearly constant. In this case (2.6) becomes
Thus, T is a decreasing function of A and r.
The question now is: when will A be nearly a constant? Recall that when r is large, a < 0 but close to zero (a → 0−). Hence Ȧ is small in magnitude but its sign depends on those of F and G. Refer to Fig. 16. Thus for large r, A is very slowly varying and the orbits of (2.6) look quite similar to a merry-go-round where the cars go round and round in closed circuits but they reach a maximum when they cross the G axis and a minimum when they cross the F axis. Refer to Fig. 17 for an example of the trajectory in this case.
4) Case 4: A = 0 and r = 1
The special case of r = 1 is considered in detail in section 2. Consider the case when A = 0. In this case η* = (0, 0, 0)T and all the three eigenvalues of 𝗗g(η) are zero. Figure 18 provides examples of the evolution of perturbations for this case.
b. Small perturbation of equilibrium points in ɛF
Comparing entries in Table 1, it turns out that there is a one to one correspondence between the stability properties of equilibria in ɛA and ɛF. In particular, the properties of the eigenvalues and eigenvectors of the Jacobian of equilibria in ɛF for r > 1 (r < 1) are very similar to those of ɛA for r < 1 (r > 1). This in turn implies that the qualitative properties of the trajectories starting close to ɛF when r > 1 (r < 1) are quite similar to those starting close to ɛA when r < 1 (r > 1). In view of this similarity, to save space, we observe that all the remarks and conclusion of section 5a carry over to equilibria in ɛF.
c. Small perturbations of equilibrium points in ɛG
From (2.7) and (2.8) we have a < 0 and b > 0 for all r > 0. In view of this, from Table 1 it turns out that the eigenvalues of the Jacobian at the equilibrium point in ɛG are purely imaginary for all r > 0. This in turn implies that the qualitative properties of trajectories starting close to ɛG are quite similar to those starting close to ɛA when r > 1 or ɛF for r < 1. Hence all the relevant conclusions from section 5a carry over to this case as well.
6. Properties of orbits starting from an arbitrary point
a. Case 1
This ratio depends on the initial point η = (A, F, G)T and the value of r. Following Lorenz (1960) we define three types of initial conditions: R1(l) = 1 is called critical and R1(l) < (>)1 is called sub (sup) critical cases. From (6.1) it can be verified that
Since G2/F2 is nonnegative, we consider two subcases.
1) 0 < r < 1
In this case β is positive. Hence, 0 < R1(l) < 1 for all −β < δ < 0 and R1(l) > 1 when δ > 0 with R1(l) = 1 when δ = 0. As an example, if r = ½, then β = 15. Thus, R1(l) is subcritical for all −15 < δ < 0, and supercritical for all δ > 0.
The cross sections of the V and E ellipsoids by the condition (6.3) are given by
Combining these with (6.1) we get
which from (6.6) is consistent only when δ = 0. This argument leads to the following inescapable conclusion: orbits starting from initial points lying outside of the pair of straight lines given by
cannot approach points in ɛA.
A number of observations are in order.
- Critical points and eigenvectors of 𝗗g(η*), η* ∈ ɛA. Combining (4.25) with (6.8), the latter can be written as (c/b)F in the F–G plane (refer to Fig. 9) indeed corresponds to the eigenvectors v(1)2 = [0, 1, (c/b)]T and v(1)3 = [0, 1, −(c/b)]T. Thus, points along these two eigenvectors in the F–G plane in Fig. 9 correspond to the critical case. Similarly, the region ℛ1 and ℛ3 correspond to the subcritical and region ℛ2 and ℛ4 to the supercritical cases.
Critical points and heteroclinic orbits. Since the orbits starting from points along v(1)2 and v(1)3 approach points in ɛA, these orbits are necessarily heteroclinic orbits corresponding to the unstable case in section 5a. However, orbits starting from the sub and super critical regions can never approach an equilibrium point in ɛA.
2) r > 1
it follows from (6.1) that
Also, since β < 0, from (6.3) it follows that δ ≥ |β| > 0 for 2G2/F2 to be nonnegative. Hence this corresponds to the supercritical case and in this case orbits starting from any η = (A, F, G)T ∈ ℝ3 cannot reach an equilibrium point in ɛA. This is to be anticipated since when r > 1, any equilibrium point in ɛA is a center.
b. Case 2
Analogous to (6.2) we readily see that
Since the lhs of (6.13) is nonnegative, we consider two subcases.
1) r > 1
In this case, since α is positive, we get 0 < R2(k) < 1 (subcritical) when −α < δ < 0 and R2(k) > 1 (supercritical) when δ > 0. In the critical case (δ = 0), using (4.23), we obtain
Thus, the pair of straight lines corresponding to (6.14) are given by G = ±(c/a)A. These straight lines correspond to the eigenvectors v(2)2 = [1, 0, (c/a)]T and v(3)2 = [1, 0, −(c/a)]T. Thus, the orbits starting from points along v(2)2 and v(2)3 approach the equilibrium points η* = (0, F*, 0)T where F* = ±2υ = ±2kE. As in case 1, orbits starting from sub (super) critical regions cannot approach an equilibrium points in ɛF.
2) r < 1
it follows that R2(k) > 1 for all 0 < r < 1. Hence this corresponds to the supercritical case and orbits starting from any initial point η ∈ ℝ3 cannot approach an equilibrium point in ɛF. Again this is to be anticipated since η* ∈ ɛF is a center when r < 1.
c. Case 3
for all r, R3(l, k) < 1. Again it can be verified by following the above arguments that equilibrium η* ∈ ɛG can never be reached by orbits starting from any point η ∈ ℝ3.
7. Index cycle
Lorenz (1960) demonstrated that the set of three ODEs in (2.1) were capable of portraying the tendency of the average strength of the midlatitude westerly winds to undergo considerable variations in intensity, a phenomenon known as the index cycle (Namias 1950). The index cycle represents a cyclical transformation of the general circulation from a pattern of strong, zonal flow in the midlatitudes along with weak subtropical flow (high index) to strong westerly flow in the subtropical latitudes combined with weak midlatitude flow (low index) and back again, a variation lasting over an extended period of several weeks. The first term in the u component of the wind (2.4) is independent of x, representing the zonal flow. The latitudes of the zonal wind maxima are fixed (y = 3π/2l, 7π/2l, etc.) and the intensity (A/l) varies, therefore Lorenz (1960) noted that the variable A could represent the zonal index. Large positive values of A correspond to strong west-to-east flow at locations of u-wind maxima, while large negative values of A result in a shift in the location of the maximum zonal flow to the south by a distance π/l. An index cycle in this model will therefore be indicated by the value of A changing from Amax to −Amax and back to Amax. Physically plausible values of A are determined by assignment of an arbitrary standard unit of vorticity, for example, Lorenz (1960) set 1 unit of vorticity = 1/3 h−1. Given this, a value of A = 0.5 would correspond to a maximum zonal wind of approximately 75 m s−1.
a. Index cycles associated with heteroclinic orbits
In this work, we define the start of an index cycle to be the time when A first drops below 99% of Amax and the end of the cycle to be when A rises above 99% of Amax. The period of the model index cycle is defined as the length of time between the start and end of the cycle. Given this definition, we can examine the sensitivity of the index cycle period to changing initial conditions and values of r. Since we use the variable A to define the index cycle, we focus our analysis around the A equilibrium. The results of this sensitivity analysis are summarized in Tables 6 –8
An example of an index cycle portrayed by this model is provided in Fig. 10, given by a small F perturbation imposed upon ɛA. As noted in Fig. 9 this is a perturbation in the ℛ1 region of the F–G domain, and results in a heteroclinic orbit with properties summarized in Table 5. For several units of time, A remains near its initial value (A = +1) then at approximately t = 15, the F and G variables quickly increase while the value of A rapidly switches from +1 to −1. As A approaches −1, the values of F and G decrease to zero. The solution remains nearly the same for several units of time (around t = 50). The solution suddenly shifts again at approximately t = 85, when F increases and G decreases rapidly, and the value of A switches from −1 to +1. The period of the index cycle in this case is 69.8 units of time.
A similar example is provided in Fig. 11, given by a small G perturbation imposed upon ɛA. This is a perturbation in the ℛ2 region of the F–G domain, which results in a heteroclinic orbit with properties described in Table 5. As with the ℛ1 example, the evolution of A is consistent with the index cycle concept: A remains near its initial value (A = +1) for several units of time, switches quickly from +1 to −1, remains near the minimum value for some time, then quickly changes again from −1 to +1. However, the behavior of the F evolution in this example is similar in character to the G evolution in the previous (ℛ1) example, F changes sign as A switches from +1 to −1 and back again, while G remains positive throughout the evolution. In this case, the period of the index cycle is slightly longer (77.1 units of time) than in the previous example.
1) Effect of changing A
In the unstable situation (r < 1) near ɛA, the index cycle period decreases as the initial value of A increases (Table 6). It appears that relatively strong zonal flow (high initial A) will result in a shorter index cycle. This makes physical sense, in the cases with stronger initial zonal flow, more energy will be available for transfer to the perturbations, allowing for a faster transformation of the flow. This is also consistent with the stability analysis of ɛA described in section 4a(2). Larger values of A produce a larger force of attraction/repulsion about the saddle point, which is consistent with a short index cycle period.
2) Effect of changing r
The sensitivity of the index cycle period to variations in r is summarized in Table 7. Similarly to the previous example, for a given value of A, one expects the minimum index cycle period to occur when the magnitude of the eigenvalue is maximized. Since λ(1)2 = Abc, this would occur when the product bc is maximized. As shown in (4.10), the product bc is maximized for r = 0.6436. This is confirmed experimentally in Table 7 for values of r < 0.6 the index cycle period increases with decreasing r, and for values of r > 0.6 the index cycle increases with increasing r.
3) Effect of changing perturbation size
The size of the initial perturbation also affects the length of the index cycle (Table 8). Larger perturbations result in shorter index cycles. Since the evolution of F in the ℛ1 class of solutions is characterized as an oscillation between the small initial value and some much larger maximum, the initial perturbation fixes the minimum value of F. In order for A to remain near its maximum/minimum value (as in a long period index cycle), F and G must remain near zero for a relatively long time. This will not happen if the initial F perturbation is not near zero. The same is true for G perturbations (ℛ2 class of solutions). It is interesting to note that, in the unstable case, identical index periods are obtained for G perturbation =(c/b) × F perturbation (for r = 0.5, (c/b) = 2.7386).
b. Index cycles in the stable case
Lorenz (1960) stated that cycling of A was possible in the stable (r > 1) situation if the perturbations were large and V/E was supercritical. As shown in section 6, the critical point for r > 1 is determined by G2/A2 = (r4 − 1)/2. However, in both the subcritical and supercritical cases, the evolution of A is not consistent with the definition of an index cycle. An example of a subcritical case is shown in Fig. 3. While A varies considerably in this example, it remains positive and does not remain near its maximum or minimum value for any considerable period of time. On the other hand, the F evolution is consistent with the index cycle definition. In the supercritical situation, the behavior of A and G are switched, A varies from Amax to −Amax and G remains positive for initial G > 0 (negative for initial G < 0). While this is consistent with the index cycle definition, A remains at its maximum or minimum value for only an instant, therefore the A evolution is inconsistent with the index cycle definition in the supercritical case (for r > 1) as well.
In this paper we provide a complete portrayal of the long-term behavior of the solution of the Lorenz’s (1960) maximum simplification equations. While the solution is always periodic, it can be broadly classified into two groups—the unstable case leading to heteroclinic orbits and the stable case of periodic orbits. It turns out that this dichotomy critically depends on the value of the ratio V/E. The unstable case arise when V/E = 2l2 or 2k2 and the stable case for all other values of this ratio. We also examine the effect of various parameters on the index cycle for both the unstable and stable cases.
In the light of this analysis it should be interesting to refine the data assimilation experiment in Lakshmivarahan et al. (2003) and the properties of the moment dynamics in Epstein (1969) to include the unstable heteroclinic orbits and stable periodic orbits.
We wish to thank John Lewis for stimulating our interest in the lower-order models. When this was in progress, Dan Dobrovolschi of the National Meteorological Administration in Bucharest, Romania, suggested that interpretation of heteroclinic orbits is an interesting problem, which provided further impetus for this work. We wish to record our sincere thanks to Alex Gluhovsky, Yunheng Wang, and the two anonymous reviewers whose comments greatly improved the presentation.
Corresponding author address: S. Lakshmivarahan, School of Computer Science, University of Oklahoma, Norman, OK 73019. Email: E-mail:firstname.lastname@example.org