Eddy saturation in a reduced two-level model of the atmosphere

Eddy saturation describes the nonlinear mechanism in geophysical flows whereby, when average conditions are considered, direct forcing of the zonal flow increases the eddy kinetic energy, while the energy associated with the zonal flow does not increase. Here we present a minimal baroclinic model that exhibits complete eddy saturation. Starting from Phillips' classical quasi-geostrophic two-level model on the beta channel of the mid-latitudes, we derive a reduced order model comprising of six ordinary differential equations including parameterised eddies. This model features two physically realisable steady state solutions, one a purely zonal flow and one where, additionally, finite eddy motions are present. As the baroclinic forcing in the form of diabatic heating is increased, the zonal solution loses stability and the eddy solution becomes attracting. After this bifurcation, the zonal components of the solution are independent of the baroclinic forcing, and the excess of heat in the low latitudes is efficiently transported northwards by finite eddies, in the spirit of baroclinic adjustment.


Introduction
The equilibrium volume transport of the Antarctic Circumpolar Current (ACC) is known to be balanced by three contributions, namely the input of momentum at the ocean surface by the wind, the downward transport of this momentum by eddies and the bottom drag as the opposite force to the wind (Munk and Palmén 1951;Nadeau and Ferrari 2015). Straub (1993) was first to suggest that the ACC volume transport is independent of the wind stress at the ocean surface although the wind was assumed to be sufficiently strong. Since then more studies, using resolved turbulent ocean eddies, have confirmed this finding (e.g. Munday, Johnson, and Marshall (2013); Nadeau and Ferrari (2015)). According to those studies, an increase in wind stress yields an increase of the vertical shear in the channel and the flow becomes baroclinically unstable, thus generating stronger eddies instead of a higher volume transport, a process called eddy saturation. Marshall et al. (2017), inspired by a previous model of variability in atmospheric storm tracks (Ambaum and Novak 2014), explained the physical principles of eddy saturation using a simple model with just three ingredients, including a zonal mo-mentum budget, a closure relation between the eddy form stress and eddy energy, and an eddy energy budget. In their model, both the vertical shear and the volume transport of the ACC are predicted to be independent of the forcing by the surface wind and instead controlled by the requirement of a sufficiently unstable vertical shear to overcome the stabilizing role of the eddy energy dissipation. Moreover, the model explains the increase of eddy energy with wind stress and they conclude an analogy to the interaction between wave activity and baroclinicity in the original model of atmospheric storm tracks (Ambaum and Novak 2014).
However, the process driving the variability of storm tracks in the mid-latitude atmosphere differs from the transport in the ocean in several ways. While mechanical stress dominates the oceanic processes described above, the baroclinic instability in the atmosphere is a result of the horizontal temperature gradient between the equator and poles, which is primarily due to the presence of a stronger radiative forcing at low rather then high latitudes. At all levels of the atmosphere, the meridional temperature gradient is proportional to the vertical shear of the mean flow by thermal wind balance and is the source of available potential energy for the eddies which are associated with the storm tracks. The corresponding poleward eddy heat fluxes in the storm tracks act to weaken the baroclinicity and therefore the mean flow (e.g. Pedlosky (1979); Holton and Hakim (2013)). The latter process is often referred to as baroclinic adjustment (Stone 1978). The Lorenz energy cycle provides a comprehensive view of energetics of the climate system that takes into account forcing, dissipation, and exchange of energy between available potential and kinetic form, and between energy pertaining to the mean flow and energy pertaining to eddy motions (Lorenz 1967;Peixoto and Oort 1992;Lucarini et al. 2014).
The heuristic model proposed by Ambaum and Novak (2014) describes this baroclinic interaction between mean flow and eddy activity. In its steady state, the model predicts a two-way equilibration of storm tracks to extratropical thermal forcing and eddy friction: baroclinicity is independent of the thermal forcing but proportional to the eddy dissipation, whereas storm track activity is independent of eddy dissipation but proportional to thermal forcing of large-scale baroclinicity and is therefore reminiscent of the eddy saturation phenomenon in the ACC (Novak, Ambaum, and Harvey 2018).
Within the Lorenz energy cycle framework, one can see the eddy saturation mechanism as that a forcing acting on the zonal fields results, on average, in an increase in the reservoir of eddy energy (in both kinetic and available potential form), whereas the reservoir of zonal energy stays largely unaffected. Additionally, one can view the eddy saturation mechanism as an extreme form of the baroclinic adjustment process; see also discussion in Lucarini, Speranza, and Vitolo (2007).
The simplest model that can incorporate baroclinic processes together with diabatic heating and surface friction is Phillips' two-level quasi-geostrophic model on the βplane (Phillips 1956). The present paper uses this model to derive a set of ordinary differential equations that are able to provide a minimal yet meaningful model of the above described interaction between mean flow and eddy activity. The novelty of this work is how the stability properties of the derived models are determined. In contrast to the usual approach of normal mode instability analysis (see e.g. Pedlosky (1979); Hoskins and James (2014)), which focuses on studying the stability properties of the zonal flow, we use here methods of dynamical systems theory, which allow us to show the existence of a second attracting steady state instead of growing normal mode baroclinic instabilities. Another novel aspect is that this second steady state exhibits the above described eddy saturation properties in a model with parameterised eddies.
In section 2 we briefly review the two-level model and its reduction which follows work by Phillips (1956) and by Thompson (1987). In section 3 we introduce the nondimensionalisation used to define our reduced model and determine the steady states. The stability and dependence on relevant parameters thereof is analysed and physically explained in section 4 and compared to the normal mode baroclinic instability analysis of Phillips' two-level model in section 5. In section 6 a quality factor is introduced to describe the oscillatory behaviour of the model and section 7 summarizes the results and discusses them in comparison to the current literature.

The model equations
We use the two-level quasi-geostrophic model in pressure coordinates of Phillips (1956), consisting of two vorticity equations, coupled by a thermodynamic energy equation and including surface friction and a β-plane approximation, i.e. f = f 0 +βy. The equations are where ψ is the streamfunction, V = (u, v) the horizontal velocity vector, ω = dp/dt the vertical velocity and ∇ is the horizontal gradient. The pressure is denoted by p, where δp = 500 hPa is the pressure difference between the two vorticity levels as well as the pressure-thickness of each layer. Subscripts 1, 2, 3 and 4 denote the pressure levels 250, 500, 750, and 1000 hPa respectively. The fields are defined in the spatial domain (x, y) ∈ [0, L]× [0, W ], where L and W are the length and width of the β-plane channel. We remark that, as customary, the x coordinate is aligned with the longitude and the y coordinate is aligned with the latitude. The fields above are defined for non-negative times t ≥ 0 starting from smooth initial conditions. Clearly, we have that A(0, y, t) = A(L, y, t) ∀y ∈ [0, W ] and ∀t ≥ 0 for all fields A.
The parameters A and κ describe the eddy diffusion and the surface friction diffusion. Additionally, J is the diabatic heating and Γ = RT /pc p − ∂T /∂p the basic state static stability, where T is the temperature, c p the isobaric specific heat capacity of dry air and R the specific gas constant for dry air.
The diabatic heating J is given by the sum of radiative and diffusive contributions. The radiative contribution represents the net local radiative energy gains and losses and is chosen to vary linearly in y, where H is the mean rate of heating per unit mass for y ∈ [0, W/2] (or cooling for y ∈ [W/2, W ]). The diffusive contribution represents the effect of lateral eddy diffusion of temperature at level 2: We remark that the parameters A and κ control the dissipation processes that remove energy from the system. The parameter H controls the input of energy in the form of an increase in temperature difference between low and high latitudes, hence fuelling the Lorenz energy cycle of the system by production of zonal mean available potential energy. The competition between dissipation and forcing determines the degree of instability and turbulence of the system, which, after transients, reaches a statistically steady state, the model climate.
Since we want to obtain evolution equations for the longitudinally averaged or mean flow (denoted by overbars) and longitudinally varying disturbances or eddies (denoted by primes) separately, we introduce the zonal mean A(y, t) = 1/L L 0 dxA(x, y, t) and define the deviation from such mean as A ′ (x, y, t) = A(x, y, t) − A(y, t) for any field A.
Now the evolution equations for the mean zonal wind and mean thermal wind or shear are obtained by taking half the sum and the difference of equations (1a) and (1b) where we introduced the geostrophic assumption of horizontally non-divergent flow, yielding the geostrophic relation u = −∂ψ/∂y and v = ∂ψ/∂x for the velocities in all terms but the one containing the Coriolis parameter. Substituting equation (1c), taking the derivative with respect to y and applying the definitions and assumptions above, we arrive at where u m is the mean zonal wind and u T the mean thermal wind or mean zonal shear.
Here we additionally assumed a specific meridional shape of the mean wind, namely ∂ 2 u/∂y 2 = −λ 2 y u where λ 2 y is an unspecified wavenumber. Following Thompson (1987), we next restrict the eddy component of the barotropic and baroclinic streamfunction to be of the form ψ ′ m,T (x, y, t) = A m,T (t) sin kx sin ly + B m,T (t) cos kx sin ly, where the wavenumber k remains unspecified and l = π/W . We note that ∇ 2 ψ ′ m,T = −λ 2 ∇ ψ ′ m,T with λ 2 ∇ = k 2 + l 2 = k 2 + π 2 /W 2 . Using this assumption, all eddy product terms in equations (4a) and (4b) besides the net poleward eddy temperature flux v ′ m ψ ′ T cancel. Using the explicit expression (5) for We see that the net poleward heat transport attains its maximum at y = W/2 and vanishes at y = 0 and y = W . Indeed, this qualitatively corresponds to what is observed in more complex versions of the same model (Lucarini, Speranza, and Vitolo 2007) and to what is observed in the actual climate (Peixoto and Oort 1992). Following again Thompson (1987), we assume that the mean thermal wind u T has the same meridional structure, and inspecting equation (4a) we see that the same must be true for the mean zonal wind u m . Additionally, from the assumption about the zonal mean wind, it implies that λ 2 y = 2π 2 /W 2 . Finally, equations (4a) and (4b) can be substantially simplified to which is now a closed system by an evolution equation for the net poleward temperature flux. This equation and evolution equations for three other occurring eddy correlation terms, namely the mean meridional kinetic energy v ′2 m , the temperature variance v ′2 T and the cross-correlation between temperature and geopotential v ′ m v ′ T , are derived following to some extend the technique by Thompson (1987), with the difference that our model includes a surface friction term. The poleward temperature flux equations is The evolution equations for the three other eddy correlation terms can be found in appendix A.

The non-dimensional model
Up to this point the meridional wavenumber k of the eddy component of the streamfunctions remained unspecified. In principle, one has that k = 2pπ/L, where p is a non-vanishing natural number. Now, we let k be a multiple of the longitudinal wavenumber l and we define k = jl = jπ/W where j is a parameter describing the aspect ratio of the eddies, that is -for simplicity -assumed to be a non-negative number. Indeed, j can only assume discrete values, but since L ≫ W , they are closely spaced and so the discrete nature of j is neglected in what follows. Hence, j larger than 1 means that the eddies are elongated in the meridional direction. If j is smaller than 1 but positive, the eddies are elongated in the longitudinal direction.
To simplify the notation, we introduce the dimensionless variables and define the time variable τ = tβ/λ R , so that equations (8a), (8b), (9) and (A1a)-(A1c) take the form where the definitions of the dimensionless constants α, γ, δ, ǫ, ν and η can be found in appendix B. For a non-zero heating rate H this system exhibits three steady states. The first one is a zonal steady state where all eddy components are zero: Here, M 0 and S 0 are always positive because the denominator of the pre-factor and all parameter values including the heating rate are positive within the considered parameter space (see appendix C). The second steady state, denoted by P * = (M * , S * , T * , K * , V * , X * ), is given by where K * , V * , X * and the new parameters B, C and D can be found in appendix D. Note that B, C and D are independent of the heating H and therefore the steady state mean zonal wind and shear are independent of the heating.
The third steady state of the system can be excluded as unphysical since the values of K and V are negative in the given parameter space despite being non-dimensional variables for v ′2 m and v ′2 T , respectively, which being squared real quantities must always be non-negative. Therefore the system exhibits two physical solutions, in the following referred to as the zonal (P 0 ) and the eddy steady state (P * ).

Physical mechanisms of zonal and eddy saturated steady state
The linear stability of these two solutions is determined by numerically calculating the eigenvalues of the Jacobian matrix of system (10a)-(10f) at each steady state. Therefore, the parameters in table C1 were kept fixed, whereas the surface friction diffusion κ, the eddy diffusion A, the mean rate of heating H and the parameter j describing the eddy aspect ratio were varied within the ranges given in table C2.
In the following the notion of an attracting (stable) and repelling (unstable) steady state is used to describe the linear stability of the two steady states. Hence, if the eddy steady state is attracting, the model converges to the steady state with finite eddy contributions. This state then breaks the zonal symmetry of the zonal state: the phase of the eddies is not determined but the eddy correlation statistics are fixed in time. On the other hand, an attracting zonal steady state describes a state of the atmosphere where eddy activity decays until the flow is purely zonal, and once in that state remains zonal.
The considered parameter space is divided into two opposite stability regions where either the zonal or the eddy steady state is attracting and the respective other state is repelling. Keeping κ, A and j fixed but increasing the heating rate H, a transcritical bifurcation occurs (see for example Guckenheimer and Holmes (1983)) and the system switches stability from the zonal to the eddy steady state. From equation (11) it is clear that the mean zonal wind and shear of P 0 grow linearly with the heating H. However, this proportionality is lost for the eddy steady state, where M * and S * are independent of the heating rate. This can be seen directly from equation (12) and from the left hand side of figure 1, where contours of the mean zonal wind at the top level for the respective attracting steady state are shown. Such an insensitivity of the mean zonal wind to the forcing by the heating rate closely resembles the eddy saturation mechanism discussed in the introduction.
In contrast to that, the (non-dimensional) net poleward heat transport T * grows linearly with H (see again equation (12) and figure 1 on the right) and is non-negative in its attracting parameter region only. For values of H where the zonal steady state is attracting, T * is negative and therefore the net eddy heat transport is equatorward for the eddy state. This would imply that the system transports heat from cold towards warm regions, which suggests a condition that is thermodynamically not realisable. This further clarifies that the eddy steady state is physically irrelevant in its repelling region. Furthermore, it is clear from equation (12) that the steady state value of the net poleward heat transport only depends on the shear and incoming heat but not on the three other eddy correlation variables K * , V * and X * . Hence, the dynamics of the system are resembled by the mean zonal shear (proportional to the mean zonal wind) and the net poleward heat transport, and K * , V * and X * do not need to be considered separately.
The above described stability switch is also dependent on the eddy diffusion parameter A. This is illustrated in figure 1, where a higher eddy diffusion yields a larger heating rate at which the zonal steady state loses stability. The relation between eddy Figure 1. Contours of the mean zonal wind at the top level u 1 (left) and the net poleward eddy temperature flux (right) for the respective attracting steady state; the dashed region is the parameter space where the zonal steady state P 0 is attracting, the non-dashed region where the eddy steady state P * is attracting; note that the range of H in the right figure is four times the range in the left; κ = 4 × 10 −6 s −1 , j = 2.5. diffusion and heating rate in P 0 is determined by equation (11). Replacing S 0 by the mean poleward temperature gradient at the middle level T 2 yields (13) Since in the considered parameter space κ is much smaller than one, the function G is approximately linear in A. In physical terms, for a fixed equator to pole mean temperature gradient, a linear increase of the heating rate H can be balanced by a linearly increased eddy diffusion. In this sense, the stabilization of the zonal state against heating is always realized by an eddy heat transport, be it parameterised, as an eddy diffusion, or explicit, as in eddy heat transport.
It is also clear from (13) that the surface friction κ is not able to balance the heating and the stability analysis shows that a change in κ does not influence the critical value of heating for realistic values of κ (see figure 3 on the left).
These opposed roles of the zonal and eddy steady state describe two physically related mechanisms of the atmosphere to compensate incoming heat. In the attracting region of the zonal steady state the relatively low heating can be compensated by the eddy diffusion. However, this balance is no longer achievable for increasing H and decreasing A so that in the eddy steady state the eddies have to reach a finite size to perform poleward transport of heat and compensate the atmospheric radiative energy imbalance. In the actual atmosphere -or in more complex models -this corresponds to the case where the flow is baroclinically unstable, and baroclinic cyclones grow and decay as manifestation of an active Lorenz energy cycle and, at the same time, transport heat poleward. Although it is a similar mechanism, the compensation of the heating does not depend on the friction diffusion κ. This is explained later in this section.
Another parameter that influences the stability regions of the steady states is the aspect ratio j ∈ R + of the eddy shape. For the previous analysis j = 2.5 was chosen, which corresponds to eddies elongated by a factor of 2.5 in the meridional direction. For this aspect ratio, and even more meridionally elongated eddies, a stability switch always occurs and the larger j, the smaller is the heating rate H at this switching point. In contrast to that, for j smaller than one, corresponding to longitudinally elongated eddies, there is no exchange of stability and the zonal steady state is attracting for the whole parameter space under consideration.
For an aspect ratio above one, the value of H at which the zonal state loses stability is additionally dependent on the eddy diffusion parameter. The left hand side of figure  2 shows the values of H at which the switch occurs as a function of j for several orders of magnitude of A. For a vanishing small eddy diffusion (e.g. A on the order of 10 m 2 s −1 ) the eddy steady state is attracting for the whole range of H and an aspect ratio larger than one. Increasing A yields an exchange of stability for higher H, which was already seen in figure 1, and additionally an increasing value of j is required for the eddy steady state to become attracting.
Hence, for meridionally more elongated eddies the eddy steady state becomes attracting at a lower heating rate. It is indeed well known that, in the full quasigeostrophic flow, baroclinic instability is facilitated for these geometrical conditions (Pedlosky 1979). Only for a very large aspect ratio j above five and for a large eddy diffusion the critical value of H increases again. As before, besides this behaviour at very high aspect ratios, a change in the surface friction does not change this critical heating rate.
The right hand side of figure 2 shows the mean zonal wind at the top level of the eddy steady state as a function of the aspect ratio. Eddies close to a circular shape, i.e. j close to one, yield unrealistically high wind speeds just above 215 m s −1 . For more meridionally elongated eddies the wind speed decreases for an aspect ratio of up to 5 and thereafter increases again for very high eddy diffusion values. For j ≥ 2.5 one gets reasonably realistic values for the wind speed. However, before this upward slope, the wind speed neither depends on the eddy diffusion nor on the surface friction diffusion (the latter not shown).
For the net poleward heat transport the opposite holds: the longer the eddies in the meridional direction, the higher the heat transport until it drops again for very large j. Besides for these very high aspect ratios, the surface friction diffusion κ has again no impact on the values, whereas a very high eddy diffusion A yields a decreased heat transport for the whole range of j.
As briefly mentioned earlier, the surface friction κ could perhaps be thought to act similarly to the eddy diffusion A, taking out energy from the system and thus compensating for the heating H. However, the heating forces a vertical zonal shear by Figure 3. Contours of the mean zonal wind at the top level u 1 for realistic (left) and extremely high (right) values of the surface friction parameter κ; the dashed region is the parameter space where the zonal steady state P 0 is attracting, the non-dashed region where the eddy steady state P * is attracting; note that in the left figure the scale for the surface friction is logarithmic; A = 10 5 m 2 s −1 , j = 2.5. thermal wind balance, but would not force a zonal wind directly at the surface. The model can develop very low surface winds but increase its shear by increasing upper level zonal winds. In this sense, the surface friction does not directly compensate for the incoming heating H.
However, insofar that the heating induces surface winds, the surface friction would represent a sink for those surface winds. In that case, one would expect that a higher surface friction is able to dissipate the incoming energy and therefore moves the stability switch to higher values of the heating. This frictional adjustment can not be seen for realistic values of κ (figure 3 on the left) but is observable for extremely high surface friction values (figure 3 on the right).
Additionally, and in contrast to the eddy diffusion A, an extremely high surface friction increases the mean zonal wind to unrealistically high values. This behaviour of the mean flow was also observed by Novak, Ambaum, and Harvey (2018) and Marshall et al. (2017) for the Antarctic Circumpolar Current.

Short and long wavelength cut-off
For sufficiently large eddy diffusion A, the left hand side of figure 2 shows a short (large j) and long wavelength (small j) cut-off of the instability of the zonal solution. In other words, our model only exhibits finite-size eddies within a certain range of wavelengths and the zonal flow stabilizes for very short and long waves.
A similar result follows from the classical stability analysis of the Phillips model presented in, for example, Holton and Hakim (2013). This is done for a simpler version of the two-level model, without surface friction, eddy diffusion or external forcing due to heating, and is obtained by performing a linear perturbation analysis. Assuming wave-like solutions for the barotropic and baroclinic development of the streamfunction, they obtain a dispersion relation for the phase speed of the waves and from this derive conditions for instability of the zonal flow and the growth of perturbations. The zonal flow is stable for very short waves because they are inefficient in converting available into potential energy, and for very long waves, because they are stabilized by the differential rotation of the plane along latitudes (β-effect). This instability is similar to the repelling regime of the zonal steady state in our model and the results of Holton and Hakim (2013) are compared to ours in the following.
The first condition in Holton and Hakim (2013) for instability to occur is This is also Phillips' criterion for baroclinic instability of a two-level model, see Phillips (1954). Using the values in table C1 we obtain a minimal shear of approximately 3.65 m s −1 , meaning that for a shear below this value the zonal steady state must always be attracting. In fact, it can be shown that for values of u T below this level our model is always in the attracting region of the zonal steady state and as the shear is increased this state becomes repelling. According to Holton and Hakim (2013), the minimum value of u T for which P 0 becomes repelling occurs when k 2 = √ 2λ 2 R /2, where k is the meridional wavenumber. Using the relation k = jπ/W , this condition becomes where we used again the values in table C1. Analysing our model yields exactly that value of j at which the mean zonal shear is minimal. The corresponding minimum is u T ≈ 4.1 m s −1 and is similar to the value given by Holton and Hakim (2013). This value of j also yields the wavenumber of maximal instability of the purely zonal flow. The short wavelength cut-off observed for large eddy diffusion A in figure 2 on the left is given in Holton and Hakim (2013) by the critical wavelength For waves shorter than this value, instability can not occur and the zonal steady state remains attracting. In our model this value corresponds to j > 6.6 and matches roughly the value for which the zonal steady state becomes attracting again if the eddy diffusion is very high and the heating rate is approximately the standard value given in table C2.

Q factor
The quality factor Q is a dimensionless parameter that is used in mechanics and electronics to describe the oscillatory behaviour of a damped system. In general, the Q factor is defined as the ratio between the natural undamped frequency of a dynamical system and its damping coefficient. Therefore, a high Q factor indicates a larger number of oscillations before the system reaches the attracting steady state, in contrast to a low Q factor indicating that the damping is more dominant and the attracting steady state is reached after fewer oscillations. A Q factor below 1/2 corresponds to an overdamped system that does not oscillate at all when displaced from its steady state, but returns to it by exponential decay. This concept can be applied to the theory of damped oscillators, which satisfy the linear equationẍ where ω 2 0 is the frequency of free oscillations and λ the damping coefficient. The general solution of equation (17) is where A is a fixed complex amplitude and p = −λ/2 ± i ω 2 0 − λ 2 /4 1/2 . The generally accepted definition of the Q factor for the damped oscillator is We now consider a linearised system obeying the differential equation dx/dt = αx. This is associated with a linear damped oscillator if we take p = α = α r + iα i , where α r < 0 is the real and α i the imaginary part of the eigenvalue. With this identification, we find α r = −λ/2 and α i = ± ω 2 0 − λ 2 /4 1/2 and hence the Q factor of a linearised system can be expressed as By construction, Q ≥ 1/2 for eigenvalues with a non-zero imaginary part, which is associated with the oscillatory behaviour of the system. For real eigenvalues, the Q factor is exactly 1/2 and the system is said to be critically damped. Similar to an over-damped system, the steady state is approached without oscillations. The left hand side of figure 4 shows contours of the Q factor in the heating rate H versus eddy diffusion A plane. As expected, Q is always greater than 1/2 so the system oscillates everywhere in the considered parameter space. Furthermore, it can be seen that an increase of the eddy diffusion parameter A leads to a decrease of the Q factor. This corresponds to the fact that the diffusion acts as a damping and stabilizing factor of the system and leads to a reduction of the available energy. Therefore, the higher the diffusion, the less the initial transient relies on the eddy heat flux to compensate the heat imbalance. In contrast to that, a larger heating rate H increases the Q factor since more heat needs to be transported away from the equator and therefore more oscillatory cycles of the eddies are needed during the initial transient. The contours of Q in the heating rate versus surface friction κ plane are similar (not shown) and can be interpreted in the same way.
On the right hand side of figure 4 the Q factor is shown as a function of the eddy aspect ratio j for different values of A. The highest Q factor occurs for the smallest j, i.e. for the most symmetric eddy shape, and the larger j, the smaller is the Q factor. Hence, for eddies that are more elongated in the meridional direction, the system undergoes fewer oscillatory cycles before it reaches the attracting steady state. In physical terms, meridionally elongated eddies are more efficient in balancing the heat in the atmosphere than eddies close to symmetric. This also coincides with the findings Figure 4. LEFT: Contours of the Q factor in the heating rate H versus eddy diffusion A plane, the dashed region is the parameter space where the zonal steady state P 0 is attracting, the non-dashed region where the eddy steady state P * is attracting, κ = 4 × 10 −6 s −1 , j = 2.5; RIGHT: Q factor as a function of j for various values of the eddy diffusion A in m 2 s −1 , κ = 4 × 10 −6 s −1 , H = 3.5 × 10 −3 KJ ton −1 s −1 .
in the earlier chapters. Apart from that, the dependence on the diffusion parameter A is the same as in the contour plot on the left side.

Discussion and conclusion
We used Phillips' classical two-level quasi-geostrophic model on the β-plane (Phillips 1956) and performed a model reduction to develop a system of six ordinary differential equations describing the interaction between the mean flow and shear, and the eddy activity in the mid-latitude atmosphere. This has been achieved by imposing a specific shape on the eddy streamfunction and thus looking at the evolution of the amplitude of the atmospheric modes considered in our spectral projection.
Instead of taking the usual approach of performing a normal mode instability analysis, we used classical dynamical systems theory and determined two physically possible steady states, namely a zonal solution, where the flow is purely zonal and has no eddy contribution, and an eddy solution, where zonal flow and eddies coexist. The considered parameter space is divided into two separate regions where either the zonal or the eddy solution is attracting, while the alternative solution is unstable and repels nearby initial conditions.
The two competing states are in close correspondence with two possible physical mechanisms that allow the atmosphere to reach a steady state in such a way that the imbalance in the diabatic heating between equator to the poles is compensated. Relatively weak imbalances in the diabatic heating can be balanced by diffusion and no eddies need to develop to flatten the meridional temperature gradient. Indeed, the flow is, in this case, baroclinically stable. In this state the mean flow increases for stronger forcing by diabatic heating, because the latter imposes a stronger temperature difference between low and high latitudes. This monotonicity is lost as soon as the heating reaches a certain threshold, where the eddies start growing to compensate the incoming heat because the damping effect of the eddy diffusion is no longer sufficient. The zonal flow is then baroclinically unstable so that eddies grow and reach a finite steady-state amplitude, in such a way that the resulting poleward eddy heat flux weakens the baroclinicity and therefore the mean flow. This process is referred to as baroclinic adjustment (Stone 1978).
In this baroclinically unstable eddy steady state the zonal flow is independent of the external forcing by the heating. In contrast, the net poleward heat transport by the eddies and the eddy energy increase with incoming heat, so incoming energy in form of diabatic forcing of the zonal mean flow goes directly into the eddies in form of eddy kinetic energy, where it is dissipated. Therefore, the stability switch from zonal to eddy steady state describes a transfer of the location of dissipation of incoming energy. In terms of the Lorenz energy cycle, for sufficiently large heating, the zonal flow can no longer act as the reservoir for the incoming energy, so a new reservoir needs to develop, namely the eddies. Thus, the reduced model is exhibiting complete eddy saturation.
In contrast to Munday, Johnson, and Marshall (2013), stating that eddy-resolving models of the ACC lead to a much lower sensitivity of the volume transport to increased forcing, Marshall et al. (2017) suggested the possibility to capture the physics of eddy saturation in models with parameterised eddies. Thereafter, Mak et al. (2017) used a highly idealised model configuration of the ACC with parameterised eddies to obtain a near total eddy saturation. For the atmosphere, Novak, Ambaum, and Harvey (2018) found a model where the mean baroclinicity is nearly independent of the thermal forcing, whereas the eddy intensity responds more strongly to the forcing and is independent of the eddy dissipation. In this context, the model presented in this work is the first to exhibit complete eddy saturation in a model of the mid-latitude atmosphere with parameterised eddies.
Another aspect of the eddy saturation mechanism described in Marshall et al. (2017) is the increase of the volume transport with the bottom drag, and one would expect a similar effect of the surface friction in our model, dissipating energy from the system and therefore compensating for the incoming heat. This would lead to a stability switch at higher values of the heating and to an increase of the zonal wind speed. However, this effect can only be seen for unrealistically high values of the surface friction, perhaps indicating a limitation of our model. Due to the two-level setup the model is able to develop a high wind shear induced by the incoming heating by increasing the upper level zonal wind but maintaining a very low surface wind. Hence, the surface friction can not directly compensate for the incoming heating in the model, although in the atmosphere the heating actually induces stronger surface winds.
In addition to the zonal wind being independent of the incoming heating and the surface friction (for realistic values) in the baroclinically unstable eddy steady state, the wind speed is independent of the eddy diffusion as well. Instead, an increased diffusion leads to a decrease of eddy energy, indicating the damping effect of eddy dissipation which is known as dissipative control (e.g. Novak, Ambaum, and Harvey (2018)), and is in accordance with Marshall et al. (2017), stating that the equilibrium volume transport is controlled by the ACC requiring sufficiently strong vertical shear induced by external forcing to overcome the stabilising role of the eddy dissipation.
Despite the simplicity of the model and the limiting assumptions made on the shape of the eddies, our model exhibits complete eddy saturation, a regime of the mid-latitude atmosphere that appears to be supported by numerical models. The eddy saturation mechanism is explained as a switch of stability from a zonal mean state to a mean state with finite eddy amplitude. This work is of relevance to climate modeling, where the response of the mid-latitude storm tracks to changes in the diabatic heating is still not fully understood. 571. Straub, David N. 1993. "On the Transport and Angular Momentum Balance of Channel Models of the Antarctic Circumpolar Current." Journal of Physical Oceanography 23 (4): 776-782. Thompson, Philip D. 1987. "Large-Scale Dynamical Response to Differential Heating: Statistical Equilibrium States and Amplitude Vacillation." Journal of the Atmospheric Sciences 44 (8): 1237-1248.
Here, B, C and g are positive for the whole parameter space given by tables C1 and C2, D and f are non-negative for j larger than 1 and for the positivity of D additionally κ ≥ 3 × 10 −6 s −1 is required.