Parameterization of Negative Nonlinear Feedback using a Five-dimensional Lorenz Model

Over 50 years ago, Professor Edwards Lorenz of MIT illustrated the sensitive dependence of numerical solutions on initial conditions (ICs), known as chaos or a butterfly effect of the first kind. The equations he used have been referred to as the three-dimensional Lorenz model (3DLM). While the 3DLM has been extensively examined to illustrate the role of nonlinearity in producing chaos, higher-dimensional Lorenz models have been studied in order to understand the impact of increased degrees of nonlinearity on system stability. Compared to the 3DLM, the 5D or 6D Lorenz model (5DLM or 6DLM) requires a larger normalized Rayleigh parameter (r) for the onset of chaos (e.g., Shen, 2014; 2015). Mathematical analysis suggests that negative nonlinear feedback associated with additional small-scale modes in the 5D or 6DLM can effectively stabilize solutions. Here, to aid understanding regarding whether nonlinear feedback can be incorporated with one or two additional (parameterized) terms, a revised 3DLM is proposed for emulating the negative nonlinear feedback resolved by the 5DLM. For chaotic solutions, the critical value of the normalized Rayleigh parameter (rc) in the revised 3DLM (~52.1) is larger than that of the 3DLM (rc~24.74); and is comparable to but larger than the rc of the corresponding 5DLM (~42.9). The result suggests that the stability of the revised 3DLM is overestimated between 43<r<52 as compared to the 5DLM. For stable solutions, the revised 3DLM produces the same steady-state (equilibrium-state) solutions as those in the 5DLM, when parameters are kept the same in both models (e.g., r=35). The revised 3DLM is unable to accurately depict the evolution of transient solutions in both amplitudes and phases, as compared to the 5DLM, but the trend mode of its transient solution is comparable to that of the 5DLM. The trend mode is determined as the non-oscillatory component of the solution using the parallel version of the ensemble Empirical Model Decomposition (EMD) algorithm. In addition, with the revised 3DLM, the dependence of the equilibrium-state solution on the implementation of the parameterized term suggests the importance of emulating the nonlinear (horizontal and vertical) advection of temperature into the parameterized term. Characteristics of the solution (such as overestimated stability, lack of detailed evolution, and a comparable trend mode) also appear in a different revised 3DLM with parameterized terms using the 6DLM. Correspondence to: Bo-Wen Shen, Associate Professor, Department of Mathematics and Statistics, Center for Climate and Sustainability Studies, Computational Science Research Center, San Diego State University, 5500 Campanile Drive San Diego, CA 92182-7720, USA, E-mail: bshen@mail.sdsu.edu Received: July 05, 2015; Accepted: August 15, 2015; Published: August 19, 2015


Introduction
In his pioneering modeling paper published in 1963, Professor Edwards Lorenz of MIT illustrated the sensitive dependence of numerical solutions on initial conditions (ICs) using a very elegant set of nonlinear ordinary differential equations [1]. The numerical phenomenon of a solution's dependance on ICs, appearing as the normalized Rayleigh parameter (r) exceeds a critical value of r c , is now known as chaos or a butterfly effect of the first kind [2]. The equations are referred to as the (three-dimensional) Lorenz model (3DLM). Since the theory's introduction by Lorenz, researchers have suggested that the source of chaos in the 3DLM is nonlinearity [3][4][5][6]. Various researchers have inferred that (1) small-scale processes can have a huge impact on large-scale processes (say, leading to the generation of a tornado) and (2) increased degrees of nonlinearity may further increase the degree of chaotic behavior. In regards to the generation of large-scale systems (e.g., a tornado), the process is referred to as a butterfly effect of the second kind [2], while the sensitive dependence of numerical solutions on ICs is called a butterfly effect of the first kind. However, Pielke [4] indicated that the inferred statements may not be accurate, and research shows that nonlinear feedbacks associated with newly added modes can be positive or negative [2,7]. In other words, nonlinear feedbacks associated with small-scale modes can either stabilize or destabilize solutions, consistent with view of Lorenz [8] on the role of small-scale processes: 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 3DLM was derived from governing equations for 2D Rayleigh-Benard convection, which have two nonlinear advection terms, 2 ( , ) J ψ ψ ∇ and ( , ) J ψ θ , here, ψ is the streamfunction, and θ is the temperature perturbation. However, only the impact of the latter is included in the 3DLM. In other words, nonlinearity in the 3DLM involves only nonlinear advection of temperature. The nonlinear interaction of two wave modes, which is expressed by the Jacobian term, can generate or impact a third wave mode through a downscale (or upscale) transfer process. The third wave mode can provide feedback to the incipient wave mode(s) via its subsequent upscale (or downscale) transfer process. Downscale and upscale transfer processes associated with the third mode form a nonlinear feedback loop that can be continuously extended when new modes are continuously generated. Appendix B of Shen [2] summarizes the successive downscale and upscale transfer processes via the Jacobian term, ( , ) J ψ θ . Related discussions suggest that an ideal numerical model should include an infinite number of Fourier modes. However, practically, all available numerical models have a finite number of modes. Thus the extension of their nonlinear feedback loop is finite (and incomplete). As a result, mode truncation may impact the degree of nonlinearity and related feedback. Since chaos does appear to be the result of the inclusion of nonlinearity with a limited number of Fourier modes (i.e., a limited degree of nonlinearity), it is important to understand the impact of the increased degrees of nonlinearity in thermodynamic processes (i.e., ( , )) J ψ θ on a solution's stability. To illustrate this process, the approach, as presented here, was to derive and examine the five-and six-dimensional Lorenz models (5DLM and 6DLM). Newly introduced "nonlinear thermodynamic processes" involve the nonlinear interaction of existing modes (in the 3DLM) and new modes that are associated with either dissipation terms or a heating term. While high-dimensional Lorenz models have been derived and examined in order to understand the impact of mode truncations on a solution's stability over a period of several years [9][10][11][12], efforts have been made to discuss (i) the extension of the nonlinear feedback loop in the 5D or 6DLM; (ii) the role of the nonlinear feedback in stabilizing or destabilizing solutions; (iii) clarifications of the first and second kind of butterfly effect; and (iv) the feasibility of parameterizing negative (or positive) feedbacks into a revised 3DLM with the aim of improving system stability without the need to use higher dimensional LMs, and while understanding how parameterization may change time-varying solutions and the equilibrium state solution.
With the original 3DLM and 5DLM, the impact of improved nonlinear advection (associated with additional high-wavenumber modes) on solution stability was previously discussed [2]. Then, the 5DLM is viewed as the 3DLM with the nonlinear feedback explicitly resolved by two additional high wavenumber modes. The feedback includes nonlinear advection as well as dissipative terms associated with the new modes. In this work, an outline of how the process can be emulated in the 3D Lorenz model is provided. In Sections 2.1 and 2.2, discussion regarding the procedures required for parameterization of the nonlinear feedback using the 5DLM or 6DLM is provided. Mathematical analysis for the dissipative terms is provided in Section 2.3 and numerical methods for calculations of the solutions and ensemble Lyapunov exponents are discussed in Section 2.4. Results are presented in Section 3, followed by concluding remarks.

Revised Lorenz Models and Numerical Approaches
The revised three-dimensional Lorenz Model (3DLM) includes two general parameterization terms (e.g., Q 1 and Q 2 ), as follows: here, σ and r are the Prandtl number the normalized Rayleigh number, respectively [1,2]. Note that r is also known as the heating parameter.
( , , ) X Y Z represent the amplitudes of the solutions resolved by the Fourier modes of Lorenz [1]. In this study, the three modes in the 3DLM are referred to as primary modes. Q 1 and Q 2 represent the nonlinear feedback from small-scale processes, derived analytically from a higher-dimensional Lorenz model, e.g., 5DLM or 6DLM. As illustrated in Shen [2], nonlinear terms of the 3DLM, as well as 5DLM and 6DLM, come solely from the advection of temperature perturbations. Thus, only including additional nonlinear terms in Eqs.

3DLM with parameterization using the 5DLM (3DLM-P5d)
In the following, determination of Q 1 and Q 2 using the 5DLM, which consists of the following equations [2], is discussed: here, both b and d 0 are constants (b=8/3 and d o =19/3) and are defined in Shen (2014) and Shen (2015) [2,7]. As compared to Y and Z, 1 Y and 1 Z represent the amplitudes of temperature perturbation modes that are resolved using two additional Fourier modes with high wave-numbers [2,7]. The two modes are referred to as secondary temperature modes. As defined in the previous studies [2,7], the term "feedback" refers to the nonlinear processes that involve the secondary modes (i.e., Y 1 and/ or Z 1 ). Therefore, Eqs. (4)(5)(6) in the 5DLM can be viewed as a 3DLM with the feedback processes (i.e., XY 1 ) that result from the secondary modes. Physically, the term XY 1 represents the nonlinear advection of secondary temperature mode that provides feedback to the primary temperature mode.
By comparing the revised 3DLM (i.e., Eqs. 1-3) with Eqs. (4-6) of 5DLM, Q 1 =0 and Q 2 =XY 1 are first obtained. In Appendix A of Shen [2], XY 1 was shown to be responsible for conversion of domain averaged potential energy between Z and Z 1 mode. Thus the negative nonlinear feedback of XY 1 is also associated with the Z 1 . Next, it is beneficial to express the secondary mode (i.e., Y 1 ) in terms of the primary modes (i.e., X,Y,Z). To achieve this, Y 1 is solved by assuming a steady-state solution for secondary modes in Eqs. (7)(8), i.e., dY 1 /dτ=0 and 1 / 0 dZ dτ = , as follows: In Eq. (10), the nonlinear advection of secondary temperature mode is represented by the primary modes (X,Y,Z). With Q 1 =0 and Q 2 in Eq. 10, the revised 3DLM (Eqs. 1-3) is referred to as the 3DLM-P5d or 3DLMP5d. Note that Eq. (10) can be simplified, as follows: (1) When X is small and Z is positive, Q 2 =-qX 2 , where q is a tunable parameter. This type of parameterization with ~0.17 0.19 q − was first included in a revised 3DLM by Shen [2], who indicated that a larger r c is required for the onset of chaos. (2) When X is large (i.e., X 2 >>bd o ), Q 2 =-bZ and is the same as the dissipative term (-bZ) in Eq. (6). Although it is challenging to linearize the revised 3DLM that includes a complicated nonlinear Q 2 term (i.e. Eq. 10), the two simplified cases make it feasible to linearize Eq. (3) and to, thus, calculate Lyapunov vectors using the Gram-Schmidt reorthonormalization procedure discussed in Section 2.4. Additionally, these two simplified versions can be used to illustrate how changes in parameterizations may lead to different equilibriumstate solutions, as discussed below.
In Appendix A of Shen [2], domain-averaged kinetic energy ( ) are written, as follows: By multiplying Eqs. 1 and 2 by X and σ − , respectively, one obtains: Adding the two questions above yields: In Eq. (15), the nonlinear term (XY) associated with the primary modes is implicit and is internally responsible for energy conversion, while the nonlinear feedback term involving the secondary modes is represented by plays a role similar to bZ, which adds negative potential energy into the system (Eq. 12) and provides negative thermodynamic feedback. A steady-state solution is reached when the right hand side of Eq. (15) equals zero. The steadystate solution with 2 0 Q = is the critical point solution of the 3DLM (i.e., c c X bZ = ± ). With the general form of 2 Q in Eq. 10, a steady-state solution appears when the primary modes ( , , X Y Z ) are replaced by the analytical solutions of the critical points in the 5DLM (Eqs. 19a-c of Shen, 2014), written, as follows: Eq. (16) suggests that steady-state solutions in the 3DLMP5d are the same as those of the 5DLM. Since parameterization of nonlinear feedback is based on the assumption of steady-state solutions for the secondary modes in the 5DLM, the result is anticipated. As discussed earlier, the parameterized term ( 2 Q ) is dominated by when X is relatively small and is dominated by bZ − when X is relatively large. The two simplified cases lead to the following different steady-state solutions: 1) , respectively. Note that Z c has the same form as Eq. (16a) for the 3DLM, the revised 3DLMs, and the 5DLM. Therefore, by choosing a different Q 2 , we can change the equilibrium-state solutions for the revised 3DLM (Eqs. 1-3) (e.g., from the fixed-point solutions of the 3DLM, when Q 2 =0, to fixedpoint solutions of the 5DLM, when Q 2 is described by Eq. 10). Although proper selection of q (=1/2) in the the first simplified case may lead to the same equilibrium-state solutions as those of the second simplified case, their time-varying solutions are still different. In fact, when X is small as compared to bd o in Eq. (10), this condition poses an upper bound on r (or q) as a result of The parameterized term (Q 2 ) emulates the nonlinear feedback associated with the nonlinear advection of temperature and dissipative processes due to the secondary modes. Thus, any simplifications in Q 2 indicate changes in the nonlinear advection and/or dissipation. The above discussions suggest the dependence for both transient and steadstate solutions on the detailed implementation of a parameterization, including nonlinear advection and/or dissipative processes. Further illustration with numerical results are provided in Figure 4 and discussed in Section 3.

3DLM with parameterization using the 6DLM (3DLM-P6d)
To emulate the nonlinear feedback associated with the secondary modes, the above approach assumes a steady state solution for the secondary modes and, thus, represents the secondary modes in terms of primary modes. To further assure the validity of the approach, a test is performed with the 6DLM. Incremental changes in the number of Fourier modes from the 3DLM, to the 5DLM and 6DLM can help trace the individual and/or collective impact on solution stability. While the secondary temperature modes (Y 1 and Z 1 ) of the 5DLM, as well as the 6DLM, introduce additional nonlinear and dissipative terms, which, in turn, provide negative nonlinear feedback, the secondary streamfunction mode (X 1 ) of the 6DLM (Eqs. 4-5 of Shen (2015) [7]) introduces additional nonlinear terms and adds a heating term. Therefore, these LMs were used to examine the competing impact between an additional heating term and the negative nonlinear feedback [7]. Research has indicated the following: 1) negative nonlinear feedback in association with the secondary temperature modes plays a dominant role in providing feedback for improving the solution stability of the 6DLM; 2) the additional heating term in association with the secondary streamfunction mode may destabilize the solution, but its impact is much smaller than that of the negative nonlinear feedback. Thus, the 6DLM has a comparable but slightly smaller c r (~41.1) for the onset of chaos, as compared to the 5DLM with 42.9 c r = ( Table 1). The 5DLM and 6DLM collectively suggest different roles for smallscale processes (i.e., stabilization vs. destabilization).

Dissipation vs. Numerical Smoothing
The above discussions emphasize the role of negative nonlinear thermodynamic feedback in improving solution stability. The relationship of the nonlinear feedback to the additional nonlinear terms (i.e., J(ψ,θ)) and the dissipative term (i.e., 2 κ θ ∇ ) was discussed in Shen (2014) [2]. An important question to be answered is how these processes may be identified in real world models such as climate models and/or weather models [13][14][15]. While it takes a significant effort to thoroughly address this question, an initial attempt is made in the following discussion. The nonlinear term, J(ψ,θ), representing the horizontal and vertical advection of temperature, is commonly included in numerical models (explicitly or implicitly). While the dissipative term may be physically interpreted as an eddy dissipative process, it can be numerically linked to an averaging operator, as discussed below. For example, the dissipative term 2 θ ∇ can be simplified, as follows: and the rate of change in θ can be approximated by Here, x ∆ and ∆t represent a grid spacing and time interval, The above, which is approximated from the linear thermodynamic equation with the dissipative term (i.e., t xx θ κθ ≈ ), qualitatively suggests that the averaging operator plays a role numerically similar to the Laplacian operator ( 2 ∇ ) in contributing to dissipative processes.

Numerical approaches
Using the 4th order Runge-Kutta scheme, the original and revised Lorenz models are integrated forward in time. The value of the heating parameter, r, is varied, while other parameters, including 10 and q=0.17 [2,7]. Below, the ensemble Lyapunov exponents are first discussed then followed by time-dependent solutions obtained from the 5DLM, 3DLMP5d, and 3DLMP6d. Numerical approaches, as discussed by Shen [2], are also briefly summarized.
For the time-dependent solutions, the initial conditions are, as follows: The dimensionless time interval ( τ ∆ ) is 0.0001. The total number of time steps (N) is 500,000. 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 [16][17][18][19][20][21][22][23][24][25]. In Shen [2], the two methods implemented and tested are the trajectory separation (TS) method [23] and the Gram-Schmidt reorthonormalization (GSR) procedure [17,21]. Using given initial conditions (ICs) and a set of parameters in the LMs, the TS scheme calculates the largest LE while the GSR scheme produces "n" LEs. Here, "n" is the dimension of the 5D or 6DLM. Obtaining the n LEs using the GSR scheme requires the derivation of the so-called variational equation. Calculations are conducted with 0.0001 τ ∆ = and 10, 000, 000 N = , yielding 1, 000 To minimize dependence on the ICs, 10,000 ensemble ( 10, 000 En = ) 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 E n are used to understand the long-term average behavior of the solutions in the LMs. The eLEs calculations of the 3DLM, 5DLM and 6DLM, obtained using the two methods above, were previously discussed and compared in Shen [2,7]. A calculation of the Kaplana and Yorke fractal dimension [26] using the (three) leading eLEs from the GSR method was provided in Appendix A of Shen [7], as an additional verification. The results of the 3DLM, obtained using both methods, are presented for a comparison. Only the TS method is applied in the 3DLMP5d and 3DLMP6d because both models introduce complicated nonlinear terms, making it difficult to derive the variational equation for the GSR scheme.
For stable solutions, as shown by numerical solutions to be discussed in Section 3, the 3DLMP5d and 3DLMP6d produce the same steady-state solutions as those in the 5DLM and 6DLM, respectively. The result is expected due to the assumption of the parameterization, as discussed in Sections 2.1 and 2.2. As nonlinear feedback processes are emulated with diagnostic equations in the revised 3DLMs and as nonlinear feedback processes are resolved by prognostic equations in the 5DLM and 6DLM, the revised 3DLMs are unable to accurately depict the evolution of transient solutions in both phases and amplitudes prior to the steady state, as compared to solutions in the 5DLM or 6DLM. Therefore, using a different methodology for comparing the time-varying solutions of the revised 3DLMs and the corresponding high-dimensional LMs is important.
A common practice is to compare the time averaged quantities of two fields in place of a direct point-to-point comparison. However, it is not trivial to determine a time scale for averaging, in particular when a solution evolves at multiple scales. Therefore, to achieve the goal of comparing solutions from different LMs, the trend mode, defined as the non-oscillatory component of the solution, is calculated and compared. The core technology for this calculation is the empirical mode decomposition (EMD) [27,28]. The EMD, as first developed by Huang et al. [27], extracts a finite number (N) of the oscillatory intrinsic mode functions (IMFs) from the local extrema of the data. The residual, which represents the difference between the raw data and a sum of N IMFs, is called the trend mode, which is non-oscillatory. The ensemble EMD (EEMD) method was proposed to effectively resolve the issue of mode mixing in the EMD [28]. Although the method is promising, it requires tremendous computing resources, as computational resources are linearly proportional to the number of ensemble trials. To reduce the time to obtain solutions, we have developed a parallel version of EEMD (PEEMD) with a three-level parallelism using message-passing interface (MPI) tasks for the first two levels and with OpenMP threads in the third level [29][30][31]. The PEEMD leads to a parallel speedup of 720 using 200 eight-core processors, and is used to calculate the trend modes in Figure 4c.  function of the normalized Rayleigh parameter (r) for the 3DLM, 5DLM, and revised 3DLMs (i.e., 3DLMP5d and 3DLMP6d). Table 1 summaries the critical value (r c ) for the onset of chaos in these models. As first indicated in Shen [2], the 5DLM requires a larger r for the onset of chaos than the 3DLM. The critical value (r c ) of the 5DLM is 42.9 while the r c of the 3DLM is 24.74. It was shown that the negative nonlinear feedback via the 1 XY term can help suppress chaotic responses. 1 XY is present in the 5DLM but not in the 3DLM. Since the 3DLMP5d includes the impact of 1 XY by assuming a steady-state solution for 1 Y and 1 Z in Eqs. (7)(8), the critical value (r c ) is expected to be comparable to that of the 5DLM and larger than that of the 3DLM, as shown in Figure 1 and Table 1. Interestingly, the r c of the 3DLMP5d is 52.1, which is larger than that of the 5DLM. The result suggests that the solution of 3DLMP5d is less chaotic than that of the 5DLM, and that the transient evolution of small-scale processes by the secondary modes, which appear associated with nonzero 1 / dY dτ and 1 / dZ dτ in the 5DLM, may destabilize the solution of the 5DLM, as further discussed using Figure 4.

Numerical Results
The 6DLM was derived to illustrate the competing impact of an additional heating term and negative nonlinear feedback [7]. It was shown that the negative nonlinear feedback plays a dominant role in providing feedback, and that the 6DLM has a comparable but slightly smaller r c (41.1) as compared to the 5DLM. The negative nonlinear feedback and the impact of an additional heating term in the 6DLM are emulated in the 3DLMP6d using steady-state solutions for the secondary modes. Therefore, the 3DLMP6d with two parameterized terms is used to ensure the validity of the nonlinear thermodynamic parameterization. In Figure 1 and Table 1, the r c of the 3DLMP6d is 51.4, which is comparable, but slightly smaller, than that of the 3DLMP5d (r c =52.1). Similarly, the 6DLM has a smaller r c than that of the 5DLM. Figure 2 displays Y-Z plots with r=35 from the 3DLM, 5DLM, 3DLMP5d, and 3DLMP6d, indicating a chaotic solution for the 3DLM but stable solutions for the rest of the LMs.
To illustrate how the revised 3DLMs, as well as the 5DLM, reach a steady-state solution with r=35, we analyze the terms on the right hand side of Eq. (6), including XY, bZ, and XY 1 . Figure 3 indicates that the first two terms are balanced with the XY 1 in order to have a steadystate solution. When the feedback of XY 1 is added into the 3DLMP5d, solutions for the XY and bZ (Figure 3b) are comparable to those of the 5DLM, indicating that XY and bZ are balanced with Q 2 in the 3DLMP5d. The corresponding solutions of the 3DLMP6 (Figure 3c) are also comparable to those of the 5DLM and 3DLMP5, indicating the validity of the approach in determining and emulating the nonlinear feedback using the 6DLM.
A more detailed comparison is provided in Figure 4 where the transient solutions for r=35 and r=40 from the 5DLM, 3DLMP5 and 3DLMP6 are examined. As indicated in panels (a)-(b), all of the solutions eventually become steady, while the solution of the 5DLM displays a smaller decay rate. In the zoomed-in view ( Figure  4c), the oscillatory solution of 5DLM displays larger amplitudes as compared to those of the 3DLMP5 and 3DLMP6. A comparison of the solutions obtained from the 3DLMP5d and 3DLMP6d indicates that the latter, with an additional term (i.e., X 1 ), produces a solution with a comparable amplitude but a different phase, as compared to the former solution. The result is consistent with that of Shen [7] in that the positive feedback associated with the secondary streamfunction mode is small as compared to the negative feedback largely associated with the secondary temperature modes. To compare the evolution of different solutions, Figure 4c displays the trend mode of the solutions at the early stage (i.e., 0 5 τ ≤ ≤ ). While transient solutions from the revised 3DLMs with parameterizations (3DLMP5 and 3DLMP6) produce large differences in both phases and amplitudes, their trend modes are very close to that of the 5DLM. The result implies that the revised 3DLMs with parameterizations, which are emulated from the feedback of secondary modes, can produce comparable time averaged solutions as compared to the 5DLM or 6DLM, and lead to the same steady-state solutions when r is smaller than the r c of the 5DLM (or 6DLM).
The above discussions demonstrate that nonlinear feedback (i.e., XY 1 ) in the 5DLM can be analytically expressed in terms of original Fourier modes within the 3DLM, and that feedbacks can be emulated by introducing an additional term (i.e., 2 2 ( , , ) Q Q X Y Z = ) into the 3DLM. The term, which emulates the collective impact of the horizontal advection of temperature perturbation and dissipations associated with the high-wavenumber modes, can be analytically determined by solving for a steady-state solution for the high-wavenumber modes within the high-dimensional LM (e.g., 5DLM), as discussed in section 2.1. Using this approach, the revised 3DLM is viewed as the 3DLM with nonlinear feedback parameterized by the additional term. This approach is referred to as the parameterization of nonlinear thermodynamic feedback.

Concluding Remarks
The sensitive dependence of numerical solutions on initial conditions (ICs), appearing in the 3DLM when the normalized Rayleigh parameter (r) exceeds a critical value (~24.74 c r ), is known as chaos or butterfly effect of the first kind. Although chaotic solutions are associated with the inclusion of nonlinearity (as indicated in a comparison between the nonlinear 3DLM and linearized 3DLM), increased degrees of nonlinearity in higher-dimensional Lorenz    52.1 , which is larger than that of the 3DLM (~24.74 c r ). The value is comparable to but larger than the r c of the corresponding 5DLM (~42.9 ). Similar characteristics are observed when comparing the 3DLMP6d and 6DLM, which yield r c values of 51.4 and 41.1, respectively. The above features suggest that the stability of the 3DMLP5d (3DLMP6d) is overestimated between 43 52 r ≤ ≤ ( 41.1 51 r < < ), as compared to the 5DLM (or 6DLM). For stable solutions, the revised 3DLM produces the same steady-state (equilibrium-state) solutions as those in the 5DLM (or 6DLM) when parameters are set to be the same in the models (e.g., 35 r = ). Although the revised 3DLM (3DLMP5d or 3DLMP6d) is unable to accurately represent the evolution of transient solutions in both amplitudes and phases, as compared to the 5DLM, the trend mode of its transient solution is comparable to that of the 5DLM. The trend mode is determined using the newly developed parallel ensemble empirical mode decomposition (PEEMD) method. The characteristics of solutions such as overestimated stability, lack of detailed evolution, and a comparable trend mode are often found in the comparisons of coarse-resolution climate model and fine-resolution climate (or weather) model simulations.
This study and previous studies [2,7] have discussed the role of negative nonlinear thermodynamic feedback in improving solution stability and showed its relationship to the additional nonlinear terms (i.e., J(ψ,θ)) and the dissipative term (i.e., 2 κ θ ∇ ). A nonlinear feedback loop consists of a pair of downscale and upscale transfer processes via the Jacobian term [2], and can be continuously extended as long as new modes can be continuously generated. The extended nonlinear feedback loop may further stabilize or destabilize system solutions. The mathematical analysis provided in section 2.3 suggests that an averaging operator may contribute to dissipative processes, which play a role numerically similar to the Laplacian operator ( 2 ∇ ). As nonlinear advection terms and averaging operators are often included in climate models, we hypothesize that the nonlinear feedback of numerical averaging/smoothing may suppress chaotic responses, leading to more stable solutions. This hypothesis will be examined using higherdimensional Lorenz models and climate/weather models with the TS method in a future study. As numerical averaging may impact the phase and amplitude of transient solutions, posing a challenge for verifying time-varying numerical solutions against observations (or, comparing the transient solution obtained from a low-and high-dimensional models), an alternative method such as the PEEMD for comparing their trend modes is required. Detailed analyses and comparisons of the different Lorenz models and their application to climate model simulations are currently being planned.