Nonlinear feedback in a six-dimensional Lorenz model : impact of an additional heating term

In this study, a six-dimensional Lorenz model (6DLM) is derived, based on a recent study using a fivedimensional (5-D) Lorenz model (LM), in order to examine the impact of an additional mode and its accompanying heating term on solution stability. The new mode added to improve the representation of the streamfunction is referred to as a secondary streamfunction mode, while the two additional modes, which appear in both the 6DLM and 5DLM but not in the original LM, are referred to as secondary temperature modes. Two energy conservation relationships of the 6DLM are first derived in the dissipationless limit. The impact of three additional modes on solution stability is examined by comparing numerical solutions and ensemble Lyapunov exponents of the 6DLM and 5DLM as well as the original LM. For the onset of chaos, the critical value of the normalized Rayleigh number (rc) is determined to be 41.1. The critical value is larger than that in the 3DLM (rc∼ 24.74), but slightly smaller than the one in the 5DLM (rc∼ 42.9). A stability analysis and numerical experiments obtained using generalized LMs, with or without simplifications, suggest the following: (1) negative nonlinear feedback in association with the secondary temperature modes, as first identified using the 5DLM, plays a dominant role in providing feedback for improving the solution’s stability of the 6DLM, (2) the additional heating term in association with the secondary streamfunction mode may destabilize the solution, and (3) overall feedback due to the secondary streamfunction mode is much smaller than the feedback due to the secondary temperature modes; therefore, the critical Rayleigh number of the 6DLM is comparable to that of the 5DLM. The 5DLM and 6DLM collectively suggest different roles for small-scale processes (i.e., stabilization vs. destabilization), consistent with the following statement by Lorenz (1972): “If the flap of a butterfly’s wings can be instrumental in generating a tornado, it can equally well be instrumental in preventing a tornado.” The implications of this and previous work, as well as future work, are also discussed.


Introduction
Fifty years have passed since Lorenz published his breakthrough modeling study (Lorenz, 1963) that changed our view regarding the predictability of weather and climate (e.g., IPCC, 2007;Pielke, 2008), laying the foundation for chaos theory (e.g., Gleick, 1987;Anthes, 2011).Since the degree of nonlinearity is finite in the original Lorenz model referred to as 3DLM, the impact of increased nonlinearity on systems' solutions and/or their stability has been studied using generalized Lorenz models (LMs) with additional Fourier modes (e.g., Curry, 1978;Curry et al., 1984;Franceschini and Tebaldi, 1985;Howard and Krishnamurti, 1986;Franceschini et al., 1988;Hermiz et al., 1995;Thiffeault and Horton, 1996;Musielak et al., 2005;Chen and Price, 2006;Roy and Musielak, 2007a, b, c;Lucarini and Fraedrich, 2009).However, such studies do not provide a definite answer regarding whether or not higher-order LMs lead to more stable solutions.
Lorenz demonstrated the association of the nonlinearity with the existence of non-trivial critical points and strange attractors in the 3DLM.Shen (2014a, denoted as Shen14) recently discussed the importance of nonlinearity in both producing new modes and enabling subsequent negative feedback to improve solution stability.The feedback loop of the Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
B.-W. Shen: A six-dimensional Lorenz model 3DLM was defined by Shen14 as a pair of downscale and upscale transfer processes associated with the Jacobian function (in Eq. 2).The feedback loop has been suggested to stabilize the solution for 1 < r < 24.74 within the 3DLM, as compared to the linearized 3DLM.Extending the nonlinear feedback loop in a five-dimensional LM (5DLM) can provide negative nonlinear feedback to produce non-trivial stable critical points when 1 < r < 42.9.The negative nonlinear feedback represents the collective impact of additional nonlinear terms and dissipative terms introduced by the two additional Fourier modes of the 5DLM.In this study (and in the previous study, Shen14), the two modes are added to improve the representation of the temperature perturbation, referred to here as secondary temperature modes.Improved stability with a higher critical Rayleigh parameter was verified by linearizing the 5DLM with respect to a non-trivial critical point and then performing a stability analysis over a wide range of values in parameters (σ , r).The outcome was possible due to the analytical solutions of the critical points in the 5DLM (e.g., Shen14).The role of the negative nonlinear feedback was further verified using the revised 3DLM that parameterizes the negative nonlinear feedback to suppress chaotic responses using a nonlinear eddy dissipation term.
In addition to the negative nonlinear feedback, Shen14 indicated that a conclusion derived from lower-dimensional LMs may not be applicable in all circumstances in a higherdimensional LM.For example, although the butterfly effect (of the first kind) with dependence of solutions on initial conditions appears in the 3DLM within the range between r = 25 and 40, it does not exist in the 5DLM.Therefore, to examine whether or not small perturbations can alter largescale structure (i.e., the butterfly effect of the second kind), a model containing proper representations of multiscale processes and their nonlinear interactions is required.As a result, it would require improving the degree of nonlinearity to address the question.
In a pioneering study using the generalized LM with a large number of Fourier modes, Curry et al. (1984) suggested that chaotic responses disappeared when sufficient modes were included.Shen14 hypothesized that the system's stability in the LMs, with a finite number of modes, can be improved with additional modes that provide negative nonlinear feedback associated with additional dissipative terms.However, since new modes can also introduce additional heating term(s), the competing role of the heating term(s) with nonlinear terms and/or with dissipative terms deserves to be examined so that the conditions under which solutions become more stable or chaotic can be better understood.Results obtained from work described here and the work of Shen14 are used to address the following question: for generalized LMs, under which conditions can the increased degree of nonlinearity improve solution stability?
To achieve the goal outlined above, the 3DLM to 5DLM was previously extended in Shen14 by including the two secondary temperature modes.In this study, the 5DLM is ex-tended to the 6DLM by adding an additional mode.The additional mode is included to improve the representation of the streamfunction (e.g., Eqs. 4 and 5), and is, therefore, referred to as the secondary streamfunction mode.While the secondary temperature modes of the 5DLM (as well as the 6DLM) introduce additional nonlinear terms and dissipative terms, which, in turn, provide negative nonlinear feedback, the secondary streamfunction mode of the 6DLM introduces additional nonlinear terms and adds a heating term.The approach, using incremental changes in the number of Fourier modes, can help trace their individual and/or collective impact on solution stability.For example, since the 6DLM also contains the negative nonlinear feedback in association with secondary temperature modes, it becomes feasible to examine the role of the additional heating term in the solution's stability and its competing impact with the negative nonlinear feedback.
The presented work is organized as follows.We describe the governing equations in Sect.2.1 and present the derivations of the 6DLM in Sect.2.2.We then discuss the energy conservation of the 6DLM in the dissipationless limit in Sect.2.3, and numerical approaches for integrations of the LMs and calculations of ensemble Lyapunov exponents in Sect.2.4.In Sect.3.1, we investigate the potential impact of the additional heating term on the solution's stability by performing stability analysis near the trivial critical point.We also illustrate how the feedback loop can be extended using the secondary streamfunction mode.In Sect.3.2, numerical results obtained from the 6DLM are provided and compared to results obtained from the 5DLM.To examine the role of the secondary streamfunction mode and to identify the major nonlinear feedback term, additional numerical experiments using the 6DLM and simplified 6DLMs are compared in Sect.3.3.Then, we discuss the dependence of the solution's stability on the Prandtl number (σ ) in Sect.3.4.Concluding remarks appear at the end.Mathematical derivations of the 5DLM and 6DLM are briefly summarized in the Supplement.

The governing equations
By assuming 2-D (x, z), incompressible and Boussinesq flow, the following equations were used by Saltzman (1962) and Lorenz (1963): Here ψ is the streamfunction that gives the u = −ψ z and w = ψ x , which, respectively, represent the horizontal and vertical velocities; θ is the temperature perturbation; and T represents the temperature difference at the bottom and top boundaries.The constants g, α, ν, and κ denote the acceleration of gravity, the coefficient of thermal expansion, the kinematic viscosity, and the thermal conductivity, respectively.The Jacobian of two arbitrary functions is defined as J (A, B) = (∂A/∂x) (∂B/∂z)−(∂A/∂z) (∂B/∂x).Additionally, Based on the above partial differential equations, Lorenz (1963) introduced a system of three ordinary differential equations to illustrate the characteristics of chaotic solutions.This system is a simplified version of the one derived by Saltzman (1962).For the reader's convenience, the same symbols as those in Saltzman (1962) and Lorenz (1963) are used here.

The 6-D Lorenz model (6DLM)
To generalize the original Lorenz model, we first use the following six Fourier modes (which are also listed in Table 1 of Shen14) to derive the 6DLM: Here l and m are defined as π a/H and π/H , representing the horizontal and vertical wavenumbers, respectively, and a is a ratio of the vertical scale of the convection cell to its horizontal scale, i.e., a = l/m.The term H is the domain height, and 2H /a represents the domain width.Using these modes, ψ and θ can be represented as follows: where C 1 and C 2 are constants, R a is the Rayleigh number and R c is its critical value for the free-slip Rayleigh-Benard problem.Using Eqs. ( 5) and ( 6), solutions within the 6DLM are represented by the six spatial modes M 1 to M 6 (Eqs.3, 4) and their corresponding time-varying amplitudes (X, Y , Z, X 1 , Y 1 , Z 1 ), respectively.By comparison, Eq. (3) was used to derived the 3DLM, and Eqs. ( 3) and (4) without M 4 were used to derive the 5DLM.While the 3DLM and 6DLM (5DLM) have one horizontal wavenumber, they contain two and four vertical wavenumbers, respectively.In the text below, to facilitate discussions, M 1 and M 4 are referred to as primary and secondary streamfunction modes, respectively, M 2 and M 3 are referred to as primary temperature modes, and M 5 and M 6 are referred to as secondary temperature modes.Here, the reader should note that an implicit limitation of this approach is that nonlinear interactions among the selected modes cannot generate (impact) any new (other) modes that are not pre-selected, suggesting limited (spatial) scale interactions.While the impact of the secondary temperature modes (i.e., Y 1 and Z 1 ) on the solution's stability was discussed by Shen14 with the 5DLM, the impact of the secondary streamfunction mode (i.e., X 1 ), which introduces a heating term (rX 1 ), is the focus of the 6DLM provided here.
To transform Eqs. ( 1) and (2) into the "phase" space, a major step is to calculate the nonlinear Jacobin functions.Calculations indicate that J (ψ, ∇ 2 ψ) in Eq. ( 1) does not lead to any explicit term in the final 6DLM, or the 3DLM or the 5DLM.Here, the Jacobian term of Eq. ( 2), which is written as follows, is discussed: Note that the 3DLM only contains the first two terms on the right-hand side of Eq. ( 7), namely XY J (M 1 , M 2 ) and −XZJ (M 1 , M 3 ), while the 5DLM includes the first four terms.
After derivations, we obtain the 6DLM with the following six equations: Here, τ = κ(1 + a 2 )(π/H ) 2 t (dimensionless time), σ = ν/κ (the Prandtl number), r = R a /R c (the normalized Rayleigh number, or the heating parameter), b = 4/(1 + a 2 ), and d o = (9 + a 2 )/(1 + a 2 ).After deriving the 6DLM in the fall of 2011, the 6DLM outlined here was compared with the work by Kennamer (1995); Musielak et al. (2005); Roy and Musielak (2007a)), who obtained the same 6DLM.A more detailed analysis regarding how the system conserves energy in the dissipationless limit, as well as a comparison with the 3DLM and 5DLM, is provided in the following discussion.
The 3DLM can be obtained from the 6DLM when terms that involve (X 1 , Y 1 , Z 1 ) are neglected.Alternatively, Eqs. ( 8)-( 10) can be viewed as a 3DLM with the feedback processes that result from the three additional modes.Therefore, the 6DLM can be viewed as a coupled system that consists of the 3DLM (Eqs.8-10) and a forced dissipative system with an additional heating term (e.g., Eqs.11-13).Here, and in Shen14, unless otherwise stated, the term "feedback" refers to the nonlinear process that involves the secondary modes, namely (X 1 , Y 1 , and/or Z 1 ).The 5DLM in Shen14 can be also obtained by ignoring the X 1 and dX 1 /dτ in the 6DLM.As a result, the 6DLM can be viewed as a coupled system that consists of the 5DLM and an additional equation (i.e., Eq. 11) that introduces nonlinear feedback associated with an additional heating term (i.e., Eq. 12).

Energy conservation in the 6-D non-dissipative LM
The domain-averaged kinetic energy (KE), available potential energy (APE), and potential energy (PE) are defined (e.g., Treve and Manley, 1982;Thiffeault and Horton, 1996;Blender and Lucarini, 2013;Shen, 2014a), as follows: Through straightforward derivations, we obtain the following equations: Here . KE p contains only a portion of the total KE of the 6DLM from the primary streamfunction mode X, but represents the total KE in the 5DLM and 3DLM.In a similar manner, as follows, Equations ( 17a) and ( 18) yield the following: while Eqs.( 17b) and ( 19) lead to the following: With Eqs. ( 8)-(13) in the dissipationless limit, the time derivatives of both Eqs. ( 20) and ( 21) are zero, so both C 3 and C 4 are constants.Therefore, Eqs. ( 20) and ( 21) indicate two energy conservation laws, including the conservation of the total KE and APE (i.e.,Eq. 20).However, it should be noted that, as follows, By comparison, the two energy conservation laws of the 5DLM are written as follows: It can been shown that both C 5 and C 6 are constants.Therefore, in the 5DLM, in addition to the conservation of the KE and APE, the KE and PE are also conserved.

Numerical approaches
Using the fourth-order Runge-Kutta scheme, the original and higher-order Lorenz models are integrated forward in time.
We vary the value of the heating parameter r but keep other parameters as constants, including σ = 10, a = 1/ √ 2, b = 8/3, d o = 19/3, and a minimum value for R c = 27π 4 /4.In Figs. 1, 2, 3 and 6, the initial conditions are given as follows: The dimensionless time interval ( τ ) is 0.0001.The total number of time steps (N) is 1 000 000 in Fig. 1 and 500 000 in Figs. 2, 3, and 6, yielding a total dimensionless time (τ ) of 100 and 50, respectively.In Figs. 2 and 6, the solutions of the 3DLM and 5DLM are rescaled by the analytical solutions of their critical points (i.e., Eqs.21 and 19 of Shen14).The solutions of the 6DLM are rescaled by the critical points of the 5DLM.In Sect.3.4, the dependence of solution stability on the Prandtl number (σ ) is discussed with selected values of σ .
To quantitatively evaluate whether or not the system is chaotic, we calculate the Lyapunov exponent (LE), a measure of the average separation speed of nearby trajectories on the critical point (e.g., Benettin et al., 1980;Froyland and Alfsen, 1984;Wolf et al., 1985;Nese, 1989;Zeng et al., 1991;Eckhardt and Yao, 1993;Christiansen and Rugh, 1997;Kazantsev, 1999;Sprott, 1997Sprott, , 2003;;Ding and Li, 2007; Li and Ding, 2011).In Shen14, the two methods implemented and tested are the trajectory separation (TS) method (e.g., Sprott, 1997Sprott, , 2003)); and the Gram-Schmidt reorthonormalization (GSR) procedure (e.g., Wolf et al., 1985;Christiansen and Rugh, 1997).Here, a brief summary of how LEs are calculated using the two methods is provided.Using given initial conditions (ICs) and a set of parameters in the LMs, the TS scheme calculates the largest LE, and the GSR scheme produces "n" LEs; here "n" is the dimension of the 5-D or 6-D LM.Calculations are conducted with τ = 0.0001 and N = 10 000 000, yielding τ = 1000.To minimize the dependence on the ICs, 10 000 ensemble (En = 10 000) runs with the same model configurations but different ICs are performed, and an ensemble averaged LE (eLE) is obtained from the average of the 10 000 LEs.A large N and En are used to understand the long-term average behavior of the solutions of the LMs and simplified LMs where some terms are ignored.While eLEs calculations using the above two methods were previously discussed and compared in Shen14, here, a calculation of the Kaplan-Yorke fractal dimension (Kaplan and Yorke, 1979) using the (three) leading eLEs from the GSR method is provided in Appendix A as an additional verification.Unless stated otherwise in the main text, the largest ensemble-averaged LE (eLE) for a given r is obtained from the TS method.
To examine the collective or individual impact of the nonlinear feedback terms and to identify the major feedback that can improve numerical predictability in the 5-D and 6-D LMs, we perform additional runs using the 6DLM with additional simplifications.The experiments, as listed in Table 1, include the following: (1) case 6DLMS1 where three nonlinear terms involving X 1 are neglected and only one feedback term (XY 1 ) is retained in Eqs. ( 9) and ( 10), (2) case 6DLMS2 where only XY 1 is ignored in Eq. (10), and (3) case 6DLMS3 where rX 1 is removed from Eq. ( 12).Results from these simplified 6DLMs are presented in Sect.3.3.

Discussion
In the following sections, we discuss the impact of additional modes on solution stability.In Sect.3.1, we illustrate the potential role of the M 4 mode by performing linear stability analysis at the trivial critical point.In Sects.3.2 and 3.3, we present and compare numerical results from the 6DLM with and without simplifications to identify the major feedback process.The dependence of solution stability on the Prandtl number (σ ) is discussed in Sect.3.4.

The impact of M 4 on linear stability
In this section, we first discuss the selection of M 4 and then its impact.As indicated in Shen14 and discussed in the Supplement, the inclusion of M 5 and M 6 modes is based on the analysis of the Jacobian term, J (ψ, θ ), and can improve the representations of the temperature perturbation and the nonlinear advection of temperature.The appearance of ∂M 5 /∂x associated with the linear term ∂θ/∂x of Eq. ( 1) requires the inclusion of an M 4 mode and the ∂M 4 /∂x associated with  8)-( 13) Ignoring terms that involve X 1 in Eqs. ( 9) and (10) 5, 6 42.3 N/A Same 6DLMS2 Eqs. ( 8)-( 13) Ignoring the term −XY 1 in Eq. ( 10) 5, 6 23.9 N/A Same 6DLMS3 Eqs. ( 8)-( 13) Ignoring the term rX 1 in Eq. ( 12) 5, 6 42.1 N/A Same 5D-NLM Eqs. ( 10)-( 14) of Shen14 Ignoring dissipative terms 1 N/A N/A N/A 6D-NLM Eqs. ( 8)-( 13) T ∂ψ/∂x of Eq. ( 2) provides feedback to the M 5 mode (in Table 1 of Shen14).The M 4 mode shares the same horizontal and vertical wavenumbers as the M 5 , but has a different phase (i.e., sin(lx) vs. cos(lx) in Eq. 4).Alternatively, via the ∂θ/∂x and T ∂ψ/∂x, the M 4 and M 5 modes are linked as follows: which can be derived by linearizing Eqs. ( 11) and ( 12) at the trivial critical point.The linearized equations are decoupled with the rest of the equations on the 6DLM, suggesting that the heating term (rX 1 ) can impact other modes as well as the stability of the nonlinear 6DLM via nonlinear feedback, as discussed below.The above equations are reduced to the following: By assuming the solution Y 1 ∝ exp(βτ ), we obtain the following two roots for β: Here, β + (β − ) represents the larger (smaller) root.An unstable normal mode with β + > 0 appears when r > d 3 o .When d o = 1, the result in Eq. ( 29) can be applied to the linearized 3DLM.As d o = 19/3 and r < d 3 o (∼ 254) in this study, both β + and β − are negative and ∂β/∂r is positive.The focus is β + because the corresponding mode dominates the solution as a result of a smaller decay rate as compared to β − .β + has a minimum (i.e., the largest decay rate) as r = 0, and increases as r increases (up to 254), leading to a decreasing decay rate.In the limit of r = 0 and σ ≥ 1, the minima of Eq. ( 29) can be written as follows: The β + = −d o provides the same decay rate as the one derived directly from Eq. ( 27) with r = 0 (i.e., the removal of rX 1 ).The simple analysis indicates that the inclusion of M 4 , as a result of β + < 0 and |β + (r = 0)| < |β + (r = 0)|, can lead to a solution component with a smaller decay rate.In other words, the inclusion of rX 1 effectively reduces the dissipative impact of −d o Y 1 in Eq. ( 27).Here, the reader should note that the relative impact of r with respect to σ can be estimated using the ratio between the first and second arguments of the radical in Eq. ( 29), written as 4σ (r − d 3 o )/(σ + 1) 2 /d 3 o .The result suggests that rX1 becomes less important when a larger σ is used.
The discussions provided above illustrate how the secondary streamfunction mode (M 4 ) may impact the growth rate of Y 1 via the linear heating term (rX 1 ).Additionally, M 4 can also provide its nonlinear feedback by extending the nonlinear feedback loop of the 5DLM (as well as the 3DLM), as follows (also see Table 2 of Shen 2014a): While Eqs. ( 31) and ( 32) form a feedback loop with M 2 → M 3 → M 2 , Eqs. ( 31) and ( 33) enable another feedback loop with M 2 → M 6 → M 2 .Equations ( 32) and (33) only contain the vertical advection of temperature due to ∂M 3 /∂x = ∂M 6 /∂x = 0.The two equations suggest that both M 3 and M 6 can provide upscaling feedback to M 2 through their interaction with M 4 , leading to two terms in Eq. ( 9), i.e., dY /dτ ∝ X 1 Z − 2X 1 Z 1 .When Z 1 is close to Z/2, their collective impact may become insignificant, X 1 (Z − 2Z 1 ) ∼ 0, as compared to the other terms in Eq. ( 9).Since the former criterion can be met near the stable critical points of the 5DLM (e.g., Eq. 20b of Shen14) and since the 6DLM shares some similarities with the 5DLM, X 1 Z and −2X 1 Z 1 are neglected in the 6DLMS1, whose results are discussed in Sect.3.3.In the next section, we first compare the numerical results of the 5DLM and 6DLM.

Numerical results of the 6DLM
In this section, we discuss the numerical results of the 6DLM beginning with energy conservation laws in the dissipationless limit.The non-dissipative version of the 6DLM (5DLM) is referred to as the 6D-NLM (5D-NLM).Figure 1 provides the time evolution of the total domain-averaged kinetic energy and available potential energy (KE + APE) for both the 6D-NLM (blue) and 5D-NLM (red).While the total domain-averaged kinetic energy and potential energy (KE + PE) is shown in pink for the 5D-NLM, the kinetic energy of the primary streamfunction mode and the potential energy (KE p + PE) is shown in green for the 6D-NLM.
Using the initial conditions in Eq. ( 25 Next, we compare the normalized solutions of (Y , Z) in the 3DLM, 5DLM, and 6DLM with two different values of r.Normalization scales are defined by the critical points listed in Table 1. Figure 2a and b display the solutions from the 3DLM and 6DLM with r = 35.Although the critical value (r c ) for the onset of chaos is r c = 24.74 in the 3DLM (Lorenz, 1963), a larger value is chosen for comparison with the 6DLM.The solution of the 3DLM never reaches a steady state but oscillates irregularly with time surrounding the nontrivial critical points.In contrast, as indicated by the converged trajectory that approaches a critical point that is close to (Y /Y c , Z/Z c ) = (−1, 1), the 6DLM yields a steady-state solution.Note that the normalization scales, Y c and Z c , are the critical points of the 5DLM, because it is difficult to obtain the analytical solution of the critical points in the 6DLM and the former and latter share similarities as discussed later.The 6DLM continues to generate steady-state solutions until r is beyond 41.1 (as discussed in Fig. 4).With an r value of 42.0, the 6DLM leads to a chaotic solution with a "butterfly" pattern in Y -Z space (Fig. 2d), while the 5DLM still produces a stable solution (Fig. 2c).
In the following, we discuss the time evolution of the solutions for the 5DLM and 6DLM to examine the impact of the secondary modes on solution's stability and to identify the major feedback associated with these modes.First, we analyze the dZ/dτ (e.g., Eq. ( 10) for the 6DLM and Eq. ( 12) of Shen14 for the 5DLM) for the cases using r = 35 that have steady-state solutions.Figure 3 indicates that all of the terms with the exception of X 1 Y , in the dZ/dτ of the 6DLM, yield comparable results to their counterparts in the 5DLM, indicating that XY 1 also plays an important role in stabilizing the solution of the 6DLM as compared to the 5DLM.While the negative feedback by XY 1 was verified by parameterizing its impact as a nonlinear eddy dissipation term into the 3DLM in Shen14, further verification using the 6DLM is provided in the following section.Due to a small value of X 1 , the X 1 Y is small as compared to other terms.A small value of X 1 could also be inferred from the steady-state solution to Eq. ( 11), giving X 1 = Y 1 /d 2 Y 1 as d o = 19/3.Additionally, the time evolution of the XY suggests that a steady state in the 5DLM is reached earlier than it is in the 6DLM, consistent with the decay rate analysis in Sect.3.1.
Figure 4 provides the analysis, used to determine the critical value of r for the onset of chaos for both the 5DLM and 6DLM, of the eLEs as a function of the normalized Rayleigh parameter r.Both models produce similar distributions of the eLEs for 35 ≤ r ≤ 50, with the following features: (1) within the stable region (as eLEs < 0), the magnitude of the eLEs is relatively smaller in the 6DLM; (2) the 6DLM requires a slightly smaller r (r c ∼ 41.1) for the onset of chaos than the 5DLM (r c ∼ 42.9); and (3) in fully chaotic regions (e.g., r > 44), the eLEs of the 5DLM and 6DLM are in good agreement, with very small differences.The first two results are consistent with the stability analysis provided in Sect.3.1, suggesting that inclusion of the M 4 mode in the 6DLM may reduce the dissipative impact associated with the M 5 mode.

Numerical results of the simplified 6DLMs
In this section, we analyze the eLEs of the 6DLM with or without additional approximations to identify the major feedback term and the impact of M 4 in the 6DLM.While the 6DLM has four nonlinear feedback terms (X 1 Z and −2X 1 Z 1 in Eq. 9, and −XY 1 and −X 1 Y in Eq. 10), the 5DLM only has one term, −XY 1 .Nonlinear feedback terms are defined as the nonlinear terms involving the secondary modes (X 1 , Y 1 , and Z 1 ).Therefore, comparable eLEs between these two LMs suggest that −XY 1 may play the most significant role in providing feedback for stabilizing solutions in the 6DLM.To verify this hypothesis, additional experiments are performed with the following simplified 6DLMs: 6DLMS1, 6DLMS2 and 6DLMS3, as introduced in Sect.2.4 and listed in Table 1.While the 6DLMS1 case retains only one nonlinear feedback term, XY 1 , the 6DLMS2 case only neglects this term.By comparison, the 6DLMS3 case is designed to examine the role of the linear heating term (rX 1 ) in Eq. ( 12).The corresponding eLEs are shown in Fig. 5.The eLEs of the 6DLMS2 resemble those of the 3DLM (Fig. 5a) with the exception of  the window regions, indirectly indicating the importance of XY 1 in stabilizing the solutions in the 6DLM.With the exception of the transition regions from eLEs < 0 to eLEs > 0 over a small range of r (i.e., r ∼ 41-43), the eLEs of the 6DLMS1 and 6DLMS3 are close to those in the 6DLM and 5DLM.The r c of these two cases are determined to be 42.3 and 42.1, respectively, which are slightly larger (smaller) than r c = 41.1 (r c = 42.9) for the 6DLM (5DLM), as shown in Fig. 5b.In addition, the magnitudes of the LEs in the stable regions are determined to be relatively larger (smaller) than those in the 6DLM (5DLM).Since the 6DLMS1 ignores the nonlinear feedback terms associated with the X 1 and since the 6DLMS3 neglects the rX 1 term, the features of the 6DLMS1 and 6DLMS3, as compared to the 6DLM, also indicate that the impact of the M 4 may slightly destabilize solutions.
The eLEs represent the averaged behavior of the model's solutions over a very large timescale, so N = 10 000 000 and T = N t = 1000 (e.g., the T in Eq. ( 23) of Shen14 should approach infinity) are used.Since some of the terms in the simplified LMs (e.g., 6DLMS1-3) are ignored, it is important to check the time evolution of the solutions on a finite timescale in order to understand whether and how the solutions approach a stable critical point, or oscillate rapidly between (unstable) non-trivial critical points.To this end, we examine the r time diagram of the normalized solutions in Fig. 6, which displays the primary mode, −Y /Y c , and secondary mode, −Y 1 /Y 1c , from the 6DLM, 6DLMS1, and 6DLMS3.Here, Y c and Y 1c are the analytical solutions of the critical points from the 5DLM.Using this approach, the deviation of the normalized solutions from 1 (i.e., −Y /Y c − 1) indicates the impact of the M 4 mode that is missing in the 5DLM.In Fig. 6, the sharp gradient of the solutions with dense contour lines near the constant value of r = 43 (in black) roughly indicates the critical value of r for the onset of chaos, consistent with the analysis of the eLEs in Fig. 5 (see Table 1).In stable regions, the primary mode, −Y /Y c , evolves with time and comes within 1 ± 0.01 in each of the three cases (Fig. 6a, c, e).For the 6DLMS1 that only includes one nonlinear feedback term (XY 1 ), the values of the secondary mode, −Y 1 /Y 1c , in stable regions are also within 1 ± 0.01 (Fig. 6d).By comparison, the normalized solutions (−Y 1 /Y 1c ) for the 6DLM and 6DLMS3 are within 1 and 0.9 in the steady state, suggesting a deviation within 10 % from the corresponding critical point of the 5DLM.If we view the stable solutions of the 5DLM as the results of the control run, the 6DLM provides approximate steady-state solutions that have derivations of only around 1 % in Y and approximately 10 % in Y 1 .The above results indicate that the nonlinear terms associated with the X 1 (i.e., M 4 mode) may produce larger relative deviations in the secondary mode Y 1 (a highwavenumber mode) than in the primary mode Y (a low wavenumber mode).By comparing the 3DLM and 5DLM, Shen14 suggested that the stability of solutions in the 3DLM can be im- proved by the negative nonlinear feedback through the term (−XY 1 ), enabled by the secondary temperature modes (Y 1 and Z 1 ) in the 5DLM.The result motivated an examination of whether or not a higher-dimensional model is more stable or less chaotic (i.e., a larger critical value of r) than a lower-dimensional model.In this study, the comparison of the 5DLM and 6DLM indicates that the additional mode (M 4 ) in the 6DLM does not help increase but slightly decreases the critical value of r for the onset of chaos.In other words, the inclusion of M 4 provides positive feedback that destabilizes the solutions through the heating term (e.g., rX 1 in Eq. 12) and/or through its nonlinear interaction with other modes.Based on the results obtained from the 5DLM and 6DLM, we have demonstrated the roles of secondary modes (i.e., small-scale processes) in stabilizing and destabilizing system solutions.In addition, the collective impact of these secondary modes on the improvement of solution stability has been examined.Since the aforementioned results are obtained from the LMs with a fixed value of σ = 10, the dependence of the stability in the 6DLM on various values of σ is discussed in the next section.

Dependence of stability on σ
Previous sections discussed the stability problem only by varying the heating parameter, r.Here, we examine the dependence of solution stability on the parameter σ , and address the question of whether or not the 6DLM still requires a smaller (larger) r c for the onset of chaos than the 5DLM (3DLM) when different values of σ are used.To efficiently achieve the goal, we conduct the eLE analysis for the 6DLM using selected values of σ , and compare it with that from the 5DLM.The dependence of the 5DLM's stability on σ was previously examined by Shen14 by performing both linear stability and eLE analyses.
For comparisons, the results obtained from the stability analysis of the 5DLM and 3DLM in Shen14 are briefly summarized as follows: in Fig. 7, pink and black lines indicate the contour lines of the Re(λ) = 0 in the (σ , r) space for the linearized 3DLM and 5DLM, respectively.Since λ is the largest eigenvalue, each line describes the critical value r l c as a function of σ , where the superscript "l" of r l c indicates the local (or linear) analysis.Following each of the contour lines in the direction of increasing σ , its right (or left) hand side contains areas with negative (or positive) values of Re(λ), suggesting stable (or unstable) solutions.Therefore, unstable solutions (Re(λ) > 0) appear as r l c < r.Solid circles with the same color scheme indicate the r c determined using the eLE analysis with selected values of σ : σ = 10, 13, 16, 19, 22 and 25.Given a σ , r c is, in general, smaller than r l c in both the 3DLM and 5DLM, as previously documented (see Shen14 for additional details).
The r c of the 6DLM, with the eLE analysis, is shown in Fig. 7 with blue multiplication signs.For all of the selected cases, the critical value r c in the 6DLM is larger than that in the 3DLM, suggesting that over the range between σ = 10 ∼ 25, the 6DLM requires a larger r for the onset of chaos than the 3DLM.By comparison, in each of the selected cases with σ = 10, 13, 16, and 19, the critical value (r c ) in the 6DLM is (slightly) smaller than the one in the 5DLM.As a result, the 6DLM is less stable than the 5DLM as 10 ≤ r < 22.However, for the case with σ = 22 (or σ = 25), the r c of the 6DLM is comparable (or slightly larger), as compared to that of the 5DLM.The results may indicate a different role for the M 4 mode between σ < 22 and σ > 22, or suggest the importance of increasing the ensemble members and/or increasing the coverage of the initial conditions for the calculations of the eLEs, all of which are subject to future study.

Concluding remarks
Five-and six-dimensional Lorenz models (5DLM and 6DLM) were derived here and in Shen14 to examine the impact of additional modes on solution's stability.The 5DLM includes two new Fourier modes (i.e., the secondary temperature modes M 5 and M 6 ) that introduce the additional nonlinear and dissipative terms.The 6DLM is a super set of the 5DLM, and contains one more Fourier mode (i.e., the secondary streamfunction mode M 4 ) that introduces additional nonlinear terms and adds a heating term.The individual and collective impacts of these terms on solution stability were investigated.The 5DLM and 6DLM have comparable critical Rayleigh parameters for the onset of the chaos, and the parameters are larger than that of the 3DLM.Based on the calculations of the ensemble-averaged Lyapunov exponents (eLEs), the critical value r c for the 6DLM (5DLM) with σ = 10 is approximately 41.1 (42.9).Therefore, while the solution of the 3DLM becomes chaotic when r ranges from 25 to 40, the 6DLM (5DLM) still produces stable steady-state solutions, suggesting that predictability can be improved by the increased degree of nonlinearity.
A quantitative comparison of the eLEs from the generalized LMs with or without additional simplifications suggests the following: (1) the negative nonlinear feedback, first identified in the 5DLM and represented by XY 1 in both the 5DLM and 6DLM, plays a dominant role in providing feedback for stabilizing the solution in the 6DLM, and (2) the additional heating term (rX 1 ) associated with the M 4 mode may destabilize the solution in the 6DLM, which has a smaller r c as compared to the 5DLM.The stability analysis provided in Sect.3.1 indicates that the heating term rX 1 may effectively reduce the dissipative effect associated with the M 5 mode, and, in turn, provides effective "positive" feedback through the nonlinear feedback loop; (3) as a result of much smaller values in the X 1 , the induced destabilization (by the additional heating term) is much smaller than the induced stabilization (by the negative nonlinear feedback term).Additionally, two nonlinear feedback terms associated with M 4 nearly cancel one another (e.g., Eqs. 32 and 33).Therefore, the r c of the 6DLM is only slightly smaller than that of the 5DLM.The 5DLM and 6DLM collectively illustrate the different roles of various high-wavenumber modes in stablizing or destabilizing a system's solutions.Additional analyses of mathematical derivations and numerical results are summarized below.
As compared to the 3-D and 5-D LMs in the dissipationless limit, the 6-D non-dissipative LM also poses two energy conservation relations.One states the conservation of the total domain-averaged kinetic energy (KE) and available potential energy (APE), enabling the transfer between KE and APE.The result is consistent with the result in the 3-D and 5-D non-dissipative LMs.In contrast, the additional conservation law only provides the conservation of the domainaveraged kinetic energy associated with the primary streamfunction mode (KE p ) and the total domain-averaged potential energy (PE), instead of the total KE and PE, as compared to the 3DLM and 5DLM.The two conservations do pose constraints on all six modes of the 6DLM.However, the potential issues (e.g., whether inconsistent forcing may exist) are beyond the scope of the present study.
The competing impact of the nonlinearities and the dissipation and heating terms can be illustrated using Eq. ( 10) of the 6DLM, as follows: The first nonlinear term (XY ) and the linear term (bZ) can act as a forcing and dissipative term, respectively, in the 3-D, 5-D, and 6-D LMs.The second and third nonlinear terms (XY 1 and X 1 Y ) are introduced as additional dissipative terms by the new modes.X 1 Y is much smaller than the other terms, and XY 1 can help reach a balance with XY and bZ to stabilize the solution.The negative nonlinear feedback by XY 1 was first illustrated by Shen14 for the 5DLM.However, the feedback by XY 1 in the 6DLM may be (slightly) different from that in the 5DLM.Specifically, while XY 1 of the 5DLM includes the feedback associated with additional nonlinear and dissipative terms, XY 1 of the 6DLM includes the feedback from the additional nonlinear and heating terms such as rX 1 .
The above results provide different impacts associated with various secondary modes, consistent with Lorenz's statement in 1972, as follows: "If the flap of a butterfly's wings can be instrumental in generating a tornado, it can equally well be instrumental in preventing a tornado."The quote suggests the appearance of both positive and negative feedbacks (i.e., stabilization and destabilization) in as-sociation with various "small-scale" processes.Since mode truncation is unavoidable in finite-resolution models, the answer to the question of whether or not the feedback by new modes is positive or negative should be made in the proper context.The approach outlined here may help us understand why some generalized LMs have a larger r c while others have a smaller r c as compared to the 3DLM.For example, among the five different generalized LMs in Tables 1  and 2 of Roy and Musielak (2007c), the two LMs that include M 5 and M 6 have a r c of ∼ 40-42, comparable to the r c in the 5DLM (6DLM) outlined here.The 2 (1, 3) and 2 (0, 4) modes in Roy and Musielak (2007c) are the same as the M 5 and M 6 modes in this study.In addition, the 14-D LM, with a comparable r c (r c ∼ 43.48) described by Curry (1978), also includes these two modes 2 (1, 3) and 2 (0, 4), and does not have a vertical wavenumber higher than that of 2 (0, 4).In contrast, the 5-D LM of Roy and Musielak (2007b), which has a smaller r c (r c ∼ 22.5), does include an additional heating term, although the two additional modes are different from the secondary modes of the 5DLM and 6DLM in this study.Although preliminary analyses seem encouraging, however, detailed comparisons with other generalized LMs (e.g.,Howard and Krishnamurti, 1986;Hermiz et al., 1995;Thiffeault and Horton, 1996) are still required.In addition, the further extension of the nonlinear feedback loop is being studied with M 7 − M 9 modes, here M 7 = √ 2 sin(lx) sin(5mz), M 8 = √ 2 cos(lx) sin(5mz), and M 9 = sin(6mz).Preliminary results indicates that a larger r c is required for the onset of chaos (e.g., r c =116.9 for the 7DLM with M 8 and M 9 modes).Using a 3-D non-dissipative Lorenz model, which is shown to be a conservative system, we discussed the collective and competing impact of the nonlinear feedback loop and heating term on the energy cycle with four different regimes (e.g., Shen, 2014b).We will further analyze the energy cycle in the higher-order dissipative or non-dissipative Lorenz models using the same approach and compare the results with those using a different approach (e.g., Pelino et al., 2014).
The 5DLM and 6DLM share some similarities regarding the system's stability, but the 6DLM has one additional model.To further our understanding of the dynamics of chaos, it is required to address whether and where additional critical points may appear and impact solution's stability in the 6DLM.Due to increasing difficulties in obtaining the analytical solutions of the critical points for the 6DLM, it becomes more challenging to perform an analysis near the critical points.In addition to the analysis for examining the competing impact between the additional dissipative and heating terms, the dependence of solution's stability on the timescale (i.e., duration) of the "forcing" terms deserves additional attention.Results obtained in this study indicate eLE dependence on the number of modes (i.e., different resolutions) and resolved processes (i.e., dissipative terms or heating term).To improve our confidence in the model's long-term climate projections using high-resolution global weather or climate models (Shen et al., 2006(Shen et al., , 2012(Shen et al., , 2013)), it is important to understand whether and how the long-term stability (eLE) in the global models may be influenced by the change in a model's grid spacing as well as the resolved "forcing" associated with different physics parameterizations.Achieving this goal requires the extension or revision of the TS method for eLE calculations in the global models.Among a variety of numerical methods that are for the calculations of LEs, the TS method does not require the variational equation, which is often difficult to obtain as a result of complicated nonlinearity in physics parameterizations and/or other model components.Therefore, the TS method, which has been tested with revised 3DLMs that contain complicated nonlinear terms to parameterize the impact of negative nonlinear feedback (e.g., Shen, 2015), will be implemented in our global model to examine the impact of the model's changes (e.g., grid spacing) on the solution's stability.

Appendix A: Fractal dimension of the 6DLM
Various methods are available for calculating fractal dimensions.There are several mathematical definitions of different types of fractal dimension (Grassberger and Procaccia, 1983;Nese et al., 1987;Ruelle, 1989;Zeng et al., 1992).In this study, we only discuss the method for calculating the so-called Kaplan-Yorke dimension (D ky ), which requires the calculation of Lyapunov exponents (LEs) and thus can be used for the verification of LE calculation.The Kaplan-Yorke dimension is defined as follows (Kaplan and Yorke, 1979;Nese et al., 1987): where LE i is the ith Lyapunov exponent, and K(≤ n) is the largest integer for which K i=1 LE i ≥ 0. D ky = 0 as LE 1 < 0 and D ky = n as n i=1 LE i > 0. In this study, "n" ensemble-averaged Lyapunov exponents (eLEs), which are produced using the GSR method (e.g., Shen14), are used to estimate the corresponding D ky .The summation of all eLEs is provided in Fig. A1a, where −13.667, −30.667, and −94 are the values for the 3DLM, 5DLM and 6DLM, respectively, and are consistent with the stability analysis.For example, in the 6DLM, the summation of all eLEs should be equal to −(σ + 1 + b + d o σ + d o + 4b).The three leading eLEs for the 3DLM, 5DLM and 6DLM are provided in Fig. A1b.The corresponding fractal dimension obtained using the eLEs is provided in Fig. A2.For r = 28, the leading eLEs of the 3DLM are (0.892743 × 10 +0 , −0.701148 × 10 −3 , −0.145587 × 10 +2 ), which results in D ky = 2.06127208.The value is very close to the value of 2.063 documented in Nese et al. (1987Nese et al. ( , p. 1957) and the value of 2.062 reported by Sprott (http://sprott.physics.wisc.edu/chaos/lorenzle.htm).Here, the reader should note that the second eLE is very small but not exactly equal to zero, indicating the impact of the 10 000 different initial conditions and/or the "finite" integration time (T = 1000) in this study.The Supplement related to this article is available online at doi:10.5194/npg-22-749-2015-supplement.

Figure 1 .
Figure 1.Time evolution of energy conservation laws from the 5D-NLM and 6D-NLM.(KE + PE) and (KE + APE) are displayed for the 5D-NLM, while (KE p + PE) and (KE + APE) are shown for the 6D-NLM.(a) and (b) are for r = 25 and r = 45, respectively.All fields are normalized using the constant C o = π 2 κ 2 1+a 2 a ), the initial values of the normalized KE + APE for the 6D-NLM (Eq.20) and the 5D-NLM (Eq.23) are given as C 3 /C o and C 5 /C o , respectively, and equal to −σ/2r.C 3 /C o (or C 5 /C o ) is −0.2 for r = 25 and −0.11 for r = 45.The initial values of the normalized KE p + PE for the 6D-NLM (Eq.21) and the KE + PE for the 5D-NLM (Eq.24) are given as C 4 /C o and C 6 /C o , respectively, and both are zero.To effectively illustrate the conservation properties of the four quantities above, the time evolution of their deviations from the corresponding initial values produce four lines when plotted.Each of the lines may be shifted by a constant.For example, while the red line in Fig. 1 represents the time evolution of the deviation for KE + APE in the 5D-NLM, (i.e., KE 5-D (τ ) + APE 5-D (τ ) − KE 5-D (0) − APE 5-D (0)), the blue line with a constant shift of 0.02 represents the time evolution of the deviation for KE + PE in the 6D-NLM (i.e., KE(τ ) + APE(τ ) − KE(0) − APE(0) + 0.02).As indicated in Fig. 1, each of the four quantities is conservative.

Figure 2 .
Figure 2. (Y , Z) plots in the 3DLM (a) and 6DLM (b) with r = 35, and 5DLM (c) and 6DLM (d) with r = 42.Lorenz strange attractors appear in (a) and (d).All of the solutions are normalized by the corresponding critical points, namely, Eq. (21) of Shen14 for the 3DLM and Eq.(19) of Shen14 for the 5DLM and 6DLM.

Figure 3 .Figure 4 .
Figure 3. Forcing terms of dZ/dτ with r = 35, which are from Eq. (12) for the 5DLM of Shen (2014a) (a) and Eq.(10) for the 6DLM (b), respectively.The black and orange lines represent XY and bZ, respectively, while the blue and red lines represent XY 1 and 5X 1 Y , respectively.

Figure 5 .
Figure 5. Same as Fig. 4 except for (a) the 3DLM (in pink) and the 6DLMS2 (in orange), and (b) the 6DLMS1 (in red) and 6DLMS3 (in green).

Figure 6 .
Figure 6.The r-time diagram of numerical solutions from the 6DLM (a, b), 6DLMS1 (c, d), and 6DLMS3 (e, f).r ranges from 25 to 50 with r = 0.5.(a, c, e) show −Y/Y c , and (b, d, f) show −Y 1 /Y 1c .Y c and Y 1c are the critical points of the 5DLM as defined in Eq. (19) of Shen14.The black line indicates a constant value of r = 43.

Figure 7 .
Figure 7.The r c of the 6DLM as a function of σ .The r c , shown by blue multiplication signs (X), is determined by the eLEs of the nonlinear 6DLM.The pink and black lines indicate a constant contour of Re(λ) = 0 for the linear 3DLM and 5DLM, respectively, indicating the corresponding r c , based on a linear stability analysis.Solid circles with the same color scheme indicate a r c , determined by the eLE analysis with r = 0.1 in the corresponding nonlinear LM.

Figure A1 .Figure A2 .
Figure A1.The summation of all ensemble-averaged Lyapunov exponents (eLEs) in the LMs (a) and three leading ensembleaveraged Lyapunov exponents (eLEs) as a function of the normalized Rayleigh number (r) (b).The pink, black, and blue lines indicate the eLEs for the 3-D, 5-D and 6-D LMs, respectively.The solid, dotted, and dashed lines display the first, second and third eLEs, respectively.In (a), the pink, black, and blue lines are shifted with a constant value of 13.667, 30.667 + 0.02 and 94.0 + 0.04, respectively.To save computational resources, the eLEs of the 5-D and 6-D LMs are calculated over a shorter range of values for r (i.e., 35 < r < 50).

Table 1 .
A list of numerical experiments for different Lorenz models.The column "Modifications" indicates additional changes in the "Equations".The r c and r l c are determined based on the eLE analyses and the linear stability analysis, respectively.Solutions in "Figures" are rescaled using the factors listed in the "Scaling factors".* For the 3DLM, the ensemble averaged LE is 1.2 × 10 −2 at r = 23.7, and becomes 0.26 at r = 24.The 5-D and 6-D non-dissipative Lorenz models (5D-NLM and 6D-NLM) are used to examine the energy conservation properties.