Mechanics and Thermodynamics of a New Minimal Model of the Atmosphere

The understanding of the fundamental properties of the climate system has long benefitted from the use of simple numerical models able to parsimoniously represent the essential ingredients of its processes. Here we introduce a new model for the atmosphere that is constructed by supplementing the now-classic Lorenz '96 one-dimensional lattice model with temperature-like variables. The model features an energy cycle that allows for conversion between the kinetic and potential forms and for introducing a notion of efficiency. The model's evolution is controlled by two contributions - a quasi-symplectic and a gradient one, which resemble (yet not conforming to) a metriplectic structure. After investigating the linear stability of the symmetric fixed point, we perform a systematic parametric investigation that allows us to define regions in the parameters space where at steady state stationary, quasi-periodic, and chaotic motions are realised, and study how the terms responsible for defining the energy budget of the system depend on the external forcing injecting energy in the kinetic and in the potential energy reservoirs. Finally, we find preliminary evidence that the model features extensive chaos. We also introduce a more complex version of the model that is able to accommodate for multiscale dynamics and that features an energy cycle that more closely mimics the one of the Earth's atmosphere.


Introduction
The climate is a non-equilibrium system whose dynamics is primarily driven by the uneven absorption of solar radiation, which is mainly absorbed near the surface and in the tropical latitudes, rather than aloft and in the mid-high latitudes, respectively. The system reacts to such an inhomogeneity in the local energy input through a complex set of instabilities and feedbacks affecting its dynamical processes and thermodynamic and radiative fluxes. Such processes lead to an overall reduction in the temperature gradients inside the the system and allow for the establishment of approximate steady state conditions [1,2].
An example of the re-equilibration mechanism can be described as follows. The large scale energy transport is mainly performed by atmospheric disturbances in the form of synoptic and (to a lesser extent) planetary eddies which are, in turn, fuelled by the baroclinic conversion of available potential energy into kinetic energy, which is then dissipated by friction. In turn, the presence of large low-high latitudes temperature gradients is responsible for the presence of a reservoir of available potential energy, which is continuously replenished thanks to the dishomogeneity of the radiative energy budget across the globe [1]. This is the core of the celebrated Lorenz energy cycle, which provides a powerful representation of the climate as an engine [3]. The thermodynamic viewpoint on the climate allows to define its efficiency and irreversibility [4,5,6,7,8,9].
The reconstruction, interpretation, and analysis of observative data (now facilitated by the recent advances in data science); analytical tools borrowed (mostly) from mathematics and physics; and numerical simulations all contribute to our understanding of the climate. This task is exceedingly demanding since the system features non-trivial variability on a vast range of temporal and spatial scales and, furthermore, our ability to observe it has changed enormously over time. Additionally, the presence of periodic as well as irregular fluctuations in the boundary conditions do not allow the climate to reach an exact steady state [10,11].
One of the features of the numerical investigation of the climate system is the reliance on hierarchies of models. In other terms, climate phenomena are investigated using a full range of models going from low dimensional ones to state-of-the-art Earth system models, able to represent with higher precision many aspects of the climate. See a discussion of the meaning and of the use of hierarchies of climate models in [12,13,14,11]. It is important to remark that, given the multiscale nature of the climate system, the heterogeneity of its subdomains, and the number of the active physical, chemical, and biological processes, the endeavour of construcing a model able to directly simulate all of them appears as a Sisyphean task, whilst, instead, the parametrization of the effect of the unresolved scales on those that are explicitly simulated is an essential component of any reasonable model of the Earth system [15,16,17].
Low-order models have played and still play today a very important role for improving our understanding of the geophysical flows. Apart from the landmark 3-dimensional model developed by Lorenz in 1963 [18] starting from the truncation of the equations describing the Rayleigh-Benard convection introduced by Saltzman [19], simple models have been key to scientific advances in oceanography [20,21,22], dynamical meteorology [23,24,25], climate dynamics [26,27,28,29], turbulence [30,31,32,33] and convection [34,35], among others. Additionally, low-order models supplemented by stochastic forcings have also provided the backbone of the stochastic theory of climate [36,37,38]. Aside from sheer mathematical-related curiosity, many of these models were created in order to shed some light on specific problems by using a physically-meaningful benchmark tool that is easier to analyse mathematically and faster to simulate numerically with respect to more complex models.

The Lorenz '96 model
Of special relevance for the present study is the now-celebrated Lorenz '96 model [25,39] (hereafter L96), whose structure is briefly recapitulated below. The model consists of a lattice of N gridpoints, whose state is described by a variable. The model has periodic boundary conditions, so that it can loosely be interpreted as describing the properties of the atmosphere along a latitudinal circle. The model features in an extremely simplified -almost metaphorical -way the main processes of the atmosphere: forcing, dissipation, and advection. Two versions of the model have been proposed: a one-level version, where the dynamics takes place on a single length scale, and a two-level version, where the lattice is augmented in order to describe dynamical processes on smaller spatial and temporal scales. The two-level version of the L96 model has the especially attractive property that the time-scale separation between the fast and the slow variables can be controlled by modulating one parameter.
The one-level L96 model can be written as: with boundary conditions: k = 1, . . . , N is the index of the gridpoints defining the lattice, the nonlinear term defines a nontrivial process of advection, F is an external forcing, and γ (usually taken with unitary value) modulates the dissipation. In the unforced and inviscid regime -i.e. setting F = γ = 0, the energy of the system, expressed as the sum of the squares of the variables, is conserved: If γ = 1 and N 1, the model's attractor is a fixed point for 0 ≤ F ≤ 8/9. As F is increased, the fixed point X 1 = X 2 = . . . = X K = F loses stability as the system undergoes bifurcations leading to a quasi-periodic behavior for moderate values of F and chaotic behaviour for F ≥ 4.5 [39]. The Lorenz '96 model is strongly chaotic when F ≥ 8 and its properties are extensive with the number of nodes N [59]. By and large, the mechanism of instability of the L96 model boils down to exchanges of energy between the symmetric state and the perturbations away from it, as particularly clear in the case of the linear instability analysis [39]. Nonetheless, the L96 model clearly features only one form of energy, which we may refer to as kinetic.

This Paper
In this paper we propose an extension of the L96 model whereby a second variable is attached to each gridpoint, representing, metaphorically, the local thermodynamical properties, which are advected by the dynamical variables of the L96 model, and undergo forcing and diffusion. The model proposed here features a meaningful definition of energy that includes the kinetic part already present in the L96 model plus a potential part associated to the fluctuations of the temperature in the domain. The fundamental advantage of the new model proposed here and analysed in detail in Sect. 2 is that it features an energy cycle that allows for conversion between the kinetic and potential forms of energy. Conceptually, this mirrors the change one has going from one a one-layer quasi-geostrophic model, where only barotropic processes are possible, to a two-layer quasi-geostrophic model, which, instead, features a coupling between dynamical and thermodynamic processes via baroclinic conversion [63]. Note that the Lorenz '63 model [18] and, more completely so, its extensions to higher order modes [64] do feature a nontrivial energetics where exchanges take place between potential and kinetic energy [65]. The energy of the system defines the symplectic component that contributes -together with the metric one -to defining the evolution of the model [45]. The rest of the paper is structured as follows. Section 2 provides a thorough introduction to the one-level version of the model presented here. The evolution equations are presented together with the rationale on which the model is based. Additionally, a detailed analysis of its mechanical and thermodynamic properties is carried out. Finally, the linear stability analysis of the symmetric fixed point is also described. In Sect. 3 we present the results of a large set of numerical integrations of the model, aimed at exploring its properties in a rather vast range of values of two main parameters, which control the input of energy in the kinetic and potential form. We discuss the thermodynamics of the model in terms of mean values and the fluctuations of the main terms describing the energetics of the model. Such a physical characterisation of the model is complemented by the analysis of how the first Lyapunov exponent depends on the two considered parameters, in order to be able to separate the regions where the asymptotic dynamics of the system takes place in a regular vs in a strange attractor [66,60] We will discover a non-trivial interplay between the two sources of forcings applied to the model. We then perform a preliminary analysis for assessing to what extent the system obeys extensive chaos. In Sect. 4 we summarise the main features of the model and the results obtained so far, and propose future lines of investigations. As in the case of original L96 model, the model introduced here can be formulated in a two-level fashion -see Appendix A -, with non trivial couplings among different levels and variables and with a fairly sophisticated energetics, which is conceptually rather similar to the one described by the Lorenz energy cycle in the atmosphere. The analysis of the properties of the two-level model will not be performed in this paper and will be the subject of future studies.

Model Formulation and Properties
We want to extend the standard L96 model presented in the introduction by adding a second set of variables for all the gridpoints k = 1, . . . , N. The goal is to construct a toy model able to describe in a very simple yet conceptually correct way the interaction between dynamical and thermodynamical processes of the atmosphere. The evolution equations of the model we propose in this contribution are the following: with k = 1, . . . , N. The variable θ k can be loosely interpreted as temperature at the grid-point k. The boundary conditions are defined as The variable θ k undergoes a constant forcing G, a linear dissipation term, and a nonlinear term representing, loosely speaking, the advection performed by the X variables. Additionally, θ k and X k are linearly coupled through a term proportional to α. The purpose of this coupling is to represent, in a very simplified way, the effect of correlated thermal and dynamical fluctuations on the dynamics, which allow for an exchange between kinetic and potential energy associated with thermal fluctuations, as discussed below. The introduction of a term proportional to α is the only -yet important -modification in the dynamics of the X variables for this model as compared to the classical L96 model, see Eq. 1. In what follows, we consider F, G, α, γ ≥ 0. The coupling between the X and the θ variables is constructed in such a way that in the unforced and inviscid limit (F = G = γ = 0) the total energy of the system given by the sum of its kinetic and potential components, is conserved: whereas, in general, K and P are not separately conserved. The dynamical role of the function E is explored in the next section.

Mechanics
By definition, the time derivative of any smooth observable Ψ (X 1 , . . . , X K , θ 1 , . . . , θ k ) is obtained by applying the generator of the Koopman operator to Ψ , as follows: We will now show that linear operator L can be written as the sum of a contribution coming from a symplectic (indeed, quasisimplectic, for the reasons detailed at the end of this section) term and a contribution coming from a gradient term. Indeed, we can write: where we use the Einstein convention for the indices. The function Γ defining the gradient contribution to the dynamics is: where Ξ = ∑ N k=1 X k and Θ = ∑ N k=1 θ k . it is clear that such a component describes the irreversible dynamics as it vanishes in the unforced, inviscid limit F = G = γ = 0.
We then discuss the symplectic term associated with the Poisson bracket. We have that where It is clear that the energy E is the generator of time translations according to the symplectic contribution and the antisymmetry of the Poisson brackets enforces the corresponding conservation law already discussed in Eq. 7. We remark that, in the unviscid and unforced limit the system is not Hamiltonian because the Poisson brackets do not fulfill the Jacobi identity {A, {B,C}} + {C, {A, B}} + {B, {C, A}} = 0. Because of this and of the fact that {Γ , E} = 0, the system given in Eqs. 4-5 is not metriplectic, i.e. the standard generalisation of Hamiltonian system to the dissipative case. [67]. Note that the dynamics of dissipative fluids is, instead, metriplectic [68], and so is the dynamics of the (extended) Lorenz '63 model, which is in fact derived from the Rayleigh-Bénard equations through systematic modal truncation [45]. The lack of an underlying Hamiltonian skeleton confirms the well-known fact that the L96 model cannot be easily related to any model of fluid flows.

Thermodynamics
Using Eqs. 8-9 we obtain the time evolution of the energy of the system: which implies that, at steady state 2γĒ = FΞ + GΘ , whereΨ is the long term average of the quantity Ψ . We next analyse the separate budget of the kinetic and potential energy. By inserting K and P in Eq. 8 we obtain: K P where I K (I P ) is the rate of input of kinetic (potential) energy, D K (D P ) is the dissipation rate of kinetic (potential) energy, and C P,K the conversion rate from potential to kinetic energy. Note that the components J and Y of the Poisson bracket given in Eqs. 10-11, which describe advection, do not give any net contribution. Equations 12-13 describe the energetics of the model presented in this work, which is represented by the diagram shown in Fig. 1. One can draw a parallel between the energetics of this model and the Lorenz energy cycle of the atmosphere, where, as well known, the input of energy comes almost entirely through the potential energy channel via baroclinic forcing associated with the differential heating of low versus high latitude regions [3,1]. The two-level version of the model introduced in this paper features an energetics that is conceptually closer to the one of the true atmosphere because it is able to describe energy cascades across scales on top of energy conversion processes, see Appendix A.
At steady state conditions one has 2γK = FΞ +C and 2γP = GΘ −C, which relate the size of the reservoirs of kinetic and potential energy to intensity of the acting forcings and energy exchange. We can also introduce a notion of efficiency of this model η =C/(FΞ + GΘ ) =C/(2γĒ), which relates the amount of energy exchanged between the two reservoirs of energy to the total energy input. Since X 2 k + θ 2 k ≥ 2|X k θ k | ∀k, we have that 2E ≥ |C|/α. Therefore, |η| ≤ α/γ, which provides a constraint on the efficiency of the system. Note that η is positive if, on the average, energy is converted from potential to kinetic, and negative otherwise.
Note that K ≥ 0 and P ≥ 0 by definition. As a result, if G, γ > 0 and F = 0, one hasC ≥ 0 (on the average we have an energy flux from potential to kinetic) 1 , whereas if F, γ > 0 and G = 0,C ≤ 0 (on the average we have an energy flux from kinetic to potential). If F = G = 0, instead, we have thatK =P =Ē = 0, where the origin is a stable fixed point for the system.

Linear Stability Analysis
We investigate the linear stability of the system analysed here around the fixed point corresponding to the symmetric solution X j = X k = const, 1, . . . , j, k, . . . N and θ j = θ k = const, 1, . . . , j, k, . . . N . By plugging this ansatz in Eqs. 4-5 one gets: Taking inspiration from [39], we then investigate the linear stability of this solution by substituting X k =X + A exp(σt) exp(iκk − iωt) and θ k =θ + B exp(σt) exp(iκk − iωt) in Eqs. 4-5, where A and B are complex constants, σ is a real number defining the growth rate (if positive) of the amplitude of the wave, whilst κ is the wavenumber and ω is the angular frequency of the wave. Neglecting terms that are quadratic in the wave amplitude, one obtains: We exclude the trivial solution A = B = 0 and, thanks to linearity, we set A = 1 (only the ration b = B/A is indeed relevant). We separate real and imaginary part in the previous equations and obtain: where Re{x} and Im{x} are the real and imaginary part of the complex number x, respectively. Solving the previous Eqs. 18-21 and finding the expression of σ and ω as a function of κ and of the parameters F, G, α, and γ gives the dispersion relation of the waves. Additionally, obtaining the real and imaginary part of b allows for understanding the relative amplitude of the waves in the X and θ variables. Note that If α = 0, γ = 1, and if we neglect the θ variables, we recover the result presented in [25]. Specifically, the condition σ = 0 defines the critical case of neutral wave and corresponds to the onset of a Hopf bifurcation (if, correspondingly, ω = 0).
In the classical case of the L96 model, it is possible to derive the minimal value of F such that the fixed point of the system loses stability and, correspondingly, obtain the wavelength and frequency of the emerging neutral wave. One finds that, taking a continuum approximation, F crit = 8/9, κ crit = arccos(1/4), and ω ≈= −1.291, As a result of the presence of the coupling between the X and θ variables, it is hard to find an explicit expression for the conditions supporting the presence of a neutral wave, also if one takes the special cases where one between F and G vanishes.

Results
Many are the possible scientific questions one can address regarding the model introduced above. Building on the large literature on the L96 model discussed in the introduction, and taking into account the extra features of the current model, we can mention the following lines of investigation: -Analysis of the bifurcations leading the system from fixed point to a periodic and quasi-periodic behaviour to a chaotic regime as the forcing is increased; -Systematic investigation of the predictability of the system -e.g. analysis of the finite-time and asymptotic Lyapunov exponents and the corresponding covariant Lyapunov vectors as a function of the two forcing parameters F and G; -Systematic investigation of the energetics of the system as a function of the two forcing parameters F and G; -Analysis of the signal propagation through waves in the quasi-periodic and weakly chaotic regime; -Definition of asymptotic scaling laws for the properties if the system for large values of F and G; -Detection and analysis of chaos extensivity as the number of gridpoints N → ∞; -Extension of the model to multiple scales and analysis of dynamics and of the energetics of scale-to-scale interaction.
Obviously, it is impossible to address with a high level of detail all these aspects in the present paper. Rather than focusing on one or few aspects among those above, since this is the first time this model is proposed to the scientific community, we will present some preliminary results that address partially each of the points above, in the hope of stimulating a reader into going in greater detail. Further studies by the authors that focus specifically in some of the aspects mentioned above will be reported elsewhere.
All the numerical integrations are performed using a Dormand-Prince method with adaptive time step and a spin time of 100 time units, with runs of 1000 time units. We make use of the Python module JiTCODE [69], an extension of SciPy's ODE that allows to numerically simulate ordinary differential equations, computing quantities of interest as Lyapunov exponents as well. All results have been double checked and confirmed using the MATLAB function ode45 where integrations are performed using the 4 th order Runge-Kutta integrator with adaptive time step [70].

Transition to Chaos and Predictability
A simple yet fundamentally correct way to characterise at qualitative level the dynamical properties of a system is to investigate to what extent its evolution is sensitive to its initial conditions. Roughly speaking, the first Lyapunov exponent of a system measures the asymptotic rate of growth or decay of the distance between two orbits which are initialised in the attractor of the system at infinitesimal distance from each other [71]. Similarly, one can define the sum of the first p Lyapunov exponents as defining the asymptotic average rate of growth or decay of the p-volume defined by p + 1 orbits that are initialised in the attractor of the system at infinitesimal distance from each other [66]. Indeed, for a Q dimensional continuous time dynamical system, it is possible to compute Q Lyapunov exponents λ 1 , . . . , λ Q , where the customary ordering is such that λ j ≤ λ k if k ≤ j of the [72]. If λ 1 < 0 the attractor is a fixed point, whilst if λ 1 = 0 the attractor is periodic or quasi-periodic. Finally, the presence of a positive first Lyapunov exponent is a significant evidence that the system is chaotic, and the value of such exponent determines quantitatively the rapidity with which two nearby trajectories diverge from each other. In this case, one has that there is at least one λ j = 0, j > 1, which corresponds to the direction of the flow [66,71]. Figure 2 shows the estimate of λ 1 for N = 36 as a function of F and G in the range 0 ≤ F, G ≤ 10. We remark that the Python JiTCODE module allows for the computation of the full spectrum of Lyapunov exponents using the algorithm proposed in [72]. The system has a negative λ 1 for small values of the forcings, as expected, see Fig. 2a). We remark that if G = 0 λ 1 ≤ 0 for F ≤ 1.4, whereas for the L96 model F crit = 8/9, indicating that presence of a mechanism of energy transfer between kinetic and potential energy and the presence of a new channel of dissipation (for potential energy) leads to higher stability for the system. We observe that λ 1 depends in a very nontrivial way on both F and G, as the system's behaviour depends delicately on how the energy is injected into it, because the dynamics of the X and θ variables is, in fact, quite distinct. It is extremely different to force the system through kinetic or the potential energy channel. We also observe that the theoretical prediction of F crit = √ 2 for F = G agrees with what shown in Fig. 2a).
Many other interesting features appear. Increasing F from zero to 3 while keeping G = 1.7, the asymptotic dynamics of the system changes first from quasi-periodic to a fixed point, then again to quasi-periodic, which then alternates with chaotic behaviour. Indeed, one can observe two complex tongue-like structures in Fig. 2a) for F ≥ 1, 7, G ≥ 0.5, which indicate the presence of a very nontrivial set of bifurcations for that regions of the parameters' space, defining the transition between the quasi-periodic behaviour -the green region -ad the chaotic regime -the yellowish region in Fig. 2a).
Zooming out towards a larger range of values for F and G the intuitive argument that increasing either F or G makes the system less predictable becomes more robust, even though there are regions where a destructive interference is clear (in terms of values of λ 1 between the two forms of forcing, compare the two throughs near the diagonal in Fig. 2b).

Energetics
It is useful to investigate the long-term average of the terms in Eqs.12-13 as a function of F and G, see Fig.3. The lack of equivalence between applying forcing to the X vs to the θ variables is extremely clear by looking at the conversion term (panel e).C is positive for large values of G and moderate values of F, and negative viceversa. The absolute value ofC increases with F (G) if G (F) is kept constant. The zero isoline strongly deviates from the diagonal and indicate that if F = G there is a net transfer of energy from kinetic to potential. The zero isoline ofC coincides with the ridge in the value of λ 1 shown in 2b), indicating that the condition of no net energy exchange between the two reservoirs of energy corresponds to a state where instabilities are rather strong. The zero-isoline of the efficiency η (see panel f), by definition, coincides with the one ofC. The absolute value of the efficiency grows with the asymmetry of the forcing, and peaks for moderate intensity of either F or G, suggesting -see Sect. 3.5 below -that the energy conversion becomes less efficient decreases when stronger forcings are considered.
The behaviour of the other thermodynamical quantities is somehow unsurprising, as we have that both input and dissipation of kinetic (potential) energy increase with F (G). We remark that, once again, the response of the system to the two individual forcings is quantitatively different. It should be noted that, when one considers F ≤ 2, for G ≥ 2 one has that the net input of kinetic energy is negative (with the dissipation of kinetic energy, being, by definition, positive). This indicates a very nontrivial impact of the thermodynamic variables on the dynamical ones, which are the only ones performing advection. As a result, there is an additional mechanism of energy loss for the system, whilst all the energy input takes place through the potential energy channel. Instead, when considering low values of G, the potential energy input is always positive -yet small.

Waves amidst Chaos
We highlight some qualitative features of the dynamics of the model that indicate the presence of wave-like structures amidst chaos in the regime of moderate forcing. Figure 4 shows some examples of evolution of the system of the system in the case of N = 36 sectors with F = 10, G = 0 (panel a); F = 10, G = 10 (panel b); and F = 0, G = 10 (panel c). If the forcing on the θ variables is switched off, the X variables behave similarly to the case of the L96 model, where, amidst chaos, the clear signature of a westward propagating phase velocity can be found. As can be guessed from the evolution equations, the θ variables feature weaker variability and similar pattern of the wave crests, as they are advected by the X variables and receive energy from them. The situations is qualitatively similar when both the X and θ variables are forced, but, quite naturally, the fluctuations of the θ variables are stronger than in the previous case. The situation is radically different if, instead, the forcing acts on the θ variables only: in this case, the wave crests move from the left to the right for the θ variables as well as for the X variables, which feature a much lower variability, as expected. We will further discuss in the following sections in more quantitative terms the differences emerging when forcing the X variables only, the θ variables only, or all variables.

Chaos Extensivity
Ruelle [73] proposed that systems with short-range interaction can feature extensive chaos, because large domains can be hierarchically partitioned into smaller, weakly interacting subdomains with similar properties. One way to test whether chaos extensivity is to analyse the finite-size scaling of the Lyapunov exponents. Specifically, one plots the obtained spectrum of Lyapunov exponents for different values of system size Q (in our case, Q = 2N) against the rescaled index x = ( j − 1/2)/Q and tests whether a universal curve is obtained in the limit of large values of Q [74,59]. We remark that chaos extensivity implies that the ratio between the Kaplan-Yorke dimension of the attractor [75], also referred to as Lyapunov dimension [60], and Q tends to a constant as Q → ∞.
In order to prove convincingly the extensive nature of chaos in the system analyzed here, one should test such property for all values of F and G. This is beyond the scope of this paper. Yet, preliminary results do confirm extensivity for the three reference cases F = 10, G = 0; F = G = 10; F = 0, G = 10 shown in Panels a)-c) of Fig. 4. Indeed, as shown in 4d), in these three cases, the Lyapunov exponents spectra seem to collapse to universal curves as N grows. This is especially encouraging in view of the clear evidence for chaos extensivity in the L96 model [55,59].

Scaling Laws for Strong Forcings
As thoroughly analyzed in [59], in the one-layer L96 model the average energy per unit site scales to a very high degree of approximation as F 1.33 for large values of F. The origin of such scaling law is still unknown. We report some preliminary results of scaling laws obtained for the current model in some special configurations of parameters. We have performed long integrations (1000 time units) at steady state for three set of experiments: 1. F = 2 j , j = 3, . . . , 14; G = 0 2. F = G = 2 j , j = 3, . . . , 14; 3. G = 2 j , j = 3, . . . , 14; F = 0 which correspond to applying a forcing of increasing strength on the X variables only, on both the X and the θ variables, or on the θ variables only, respectively. We obtain the following approximate asymptotic scaling laws, which are rather accurate when F and/or G are larger than 256: 1.¯,Ē ∝ F 1.33 ,Ξ ∝ F 0.33 ,Θ ∝ F −0.28 ,P = −C/2 ∝ F 0.71 , λ 1 ∝ F 0.66 ; 2.K,Ē,P ∝ F 1.33 ,P ≈ 0.7K,C < 0, |C| ∝ F 0.70 ,Ξ ,Θ ∝ F 0.33 ,Θ ≈ 0.7Ξ , f) η =C/(2γĒ) as a function of F and G.
3.P,Ē ∝ G 1.50 ,K =C/2 ∝ G 1.00 ,Θ ∝ G 0.50 ,Ξ ∝ G 0.00 , λ 1 ∝ G 0.50 ; where the uncertainty is 0.01 for all the numbers above. As clear from these scaling relations, and in agreement with what one could intuitively guess by looking at Fig. 4, it is rather different to force the system through the X or the θ variables, and the interplay between the two reservoirs of energy is far from trivial. If the forcing is applied to only one set of variables, the energy cycle is more enhanced, ceteris paribus, when the θ variables undergo the forcing. Indeed, the reservoir of total energy and the conversionC of energy between the two forms are larger than for corresponding case of forcing applied uniquely to the X variables. The behaviour of the quantitiesΞ andΘ is also extremely different in the two cases, implying a qualitatively different way the forcing impacts the spatially-coherent fluctuations of the variables. If F = G, the ratio of the average size of the two reservoirs of energy is a constant, withK being larger thatP (and, correspondingly,Ξ being larger thanΘ ), and the average flux of energy goes from kinetic to potential. In all cases, the amount of energy that is converted between the two forms becomes a negligible fraction of the total incoming energy in the limit of large forcing. In other terms, the efficiency of the model η =C/(2γĒ) tends to zero if either F or G tend to infinity, even ifC tends to infinity. It is then unsurprising that when considering the limit of large F, regardless of whether G is also increased, we obtain for the X variables results that are in agreement with what featured by the one-layer L96, compare with [59]. At this regard, a useful piece of information is obtained by looking at the properties of λ 1 in the large forcing limit. One obtains that in scenarios 1 and 2, λ 1 ∝ F 0.66 , which is again in excellent agreement with what obtained for the L96 model (including the pre-exponential factor). Scenarios 1) and 2) seem like featuring a rather similar dynamics, the main difference between the two being the strength of the fluctuations of the θ variables; compare with panels a) and b) of Fig. 4.
The growth of λ 1 with G is slower in the case F is set to zero. We can gain a qualitative understanding of the different impact on λ 1 of changes in the value of G vs F by comparing panels a), b) and c) of Fig. 4, which nonetheless describe weaker regimes of forcing (what follows stands also in the case of stronger applied forcing). We observe that the setup 3) where F vanishes is characterised by absolute instability, with little or no advection of anomalies, whereas setups 1) and 2) are characterised by convective instability, where anomalies are spatially advected [76].

Conclusions
Simple and conceptual models have proved extremely useful for better understanding the dynamics of climate as a whole as well as of its individual components. Indeed, their usefulness spans from being the testbeds for developing new methods in terms of data analysis, data assimulation, and model testing; to supporting the definition of new metrics for testing more complex models; to providing valuable insights in the basic active physical mechanisms and most prominent mathematical features.
The model presented in this paper goes in this direction and has been constructed in order to provide a new layer of physical complexity to the L96 model by adding a new variable to each gridpoint of the model. This variable can be loosely interpreted as a local temperature and allows for the establishment of a complex energetics for the system. We are also able to introduce a notion of efficiency, which is useful for studying the conversion of energy from one form to the other one. Extending previous analyses, we have provided a fairly complete analysis of the mechanics of the new model by separating a quasi-symplectic and a metric component to its dynamical structure.
We have then performed a preliminary analysis of some of the key aspects of the new model by investigating how its properties change as a function of the two parameters that control the input of kinetic and potential energy. The interplay between the two forcings is extremely non-trivial both in the weak forcing regime -where much needs to be explored regarding the transition from fixed point to quasi-periodic to chaotic asymptotic states -and in regimes associated with stronger forcing, where the system exhibits extensive chaos. We find clear evidence that forcing the system through the kinetic or potential energy channel has fundamental impacts on the nature of the flow, where absolute vs convective instability dominate if we forced the system through the potential energy vs kinetic energy channel, respectively. Similarly to the case of the L96 model, it is possible to obtain accurate power laws describing how some of the fundamental properties of the system scale with the forcing parameters in the limit of very strong forcings.
Finally, again along the lines of the L96 model, we have introduced a two-level version -see Appendix A -of the model, which allows for studying multi-scale dynamics and which features an energetics that resembles, conceptually, the one of the atmosphere, where the Lorenz energy cycle describes succinctly the input and output of energy in the kinetic and potential form as well as the conversion between the two forms and between energy compartments at small vs large scales.
The analysis presented here is only a first step in the direction of better understanding the properties of this model, which we believe has the potential of being of great interest for investigations in areas like statistical physics, nonlinear dynamics, data assimilation, mechanics, model reduction techniques, and extreme events.

A Two-level new model
Here we propose an extension of the model able to represent multiscale dynamics and energy exchanges across scales. A detailed investigation of the properties of this model will be discussed elsewhere. We present the evolution equation and discuss the energetics of the model. Mimicking the structure of the classical two-layer L96 model, we define the following evolution equations: dY j,k dt = cbY j+1,k (Y j−1,k −Y j+2,k ) − αcφ j,k − γcY j,k + f + hc b X k , dφ j,k dt = cb(Y j−1,k φ j−2,k −Y j+1,k φ j+2,k ) + αcY j,k − γcφ j,k + g + hc b θ k , with k = 1, ..., N; j = 1, ..., J, where Y j,k 's are j small-scale variables coupled with X k , similarly to the two-level L96 model, and φ j,k 's are j small-scale variables similarly coupled with θ k . Additionally, the variables Y j,k and φ j,k are also mutually coupled. The constant hc/b determines the strength of the coupling between variables at different scale, while c defines the time scale separation between the two levels and b controls the relative amplitude of the fluctuations between large and small scales. Finally, f (g) defines the forcing on the variables Y j,k (φ j,k ). The boundary conditions are the following: θ k−N = θ k+N = θ k , φ j,k−N = φ j,k+N = φ j,k , φ j−J,k = φ j,k−1 , φ j+J,k = φ j,k+1 .
Note that, choosing α = 0 and neglecting the θ and φ variables we obtain the classic two-layer L96 model. Below we show the equations for the fluxes of energy for X and θ alongside with those for Y and φ , as we have done in Section 2 for the one-level model.
The energy cycle of this model is depicted in Fig. 5 and is closely reminiscent of the Lorenz energy cycle of the atmosphere. Note that, if we add Eqs. 27 with 28, on the one hand, and Eqs. 29 with 30, on the other hand, we derive the energy budget for the total kinetic and total potential energy, respectively. Instead, if we add Eqs. 27 with 29, on the one hand, and Eqs. 28 with 30, on the other hand, we derive the budgets for the total energy at large and small scale, respectively.