Eddy memory as an explanation of intra-seasonal periodic behavior in baroclinic eddies

The baroclinic annular mode (BAM) is a leading-order mode of the eddy-kinetic energy in the Southern Hemisphere exhibiting. oscillatory behavior at intra-seasonal time scales. The oscillation mechanism has been linked to transient eddy-mean flow interactions that remain poorly understood. Here we demonstrate that the finite memory effect in eddy-heat flux dependence on the large-scale flow can explain the origin of the BAM's oscillatory behavior. We represent the eddy memory effect by a delayed integral kernel that leads to a generalized Langevin equation for the planetary-scale heat equation. Using a mathematical framework for the interactions between planetary and synoptic-scale motions, we derive a reduced dynamical model of the BAM - a stochastically-forced oscillator with a period proportional to the geometric mean between the eddy-memory time scale and the diffusive eddy equilibration timescale. Our model provides a formal justification for the previously proposed phenomenological model of the BAM and could be used to explicitly diagnose the memory kernel and improve our understanding of transient eddy-mean flow interactions in the atmosphere.


I. INTRODUCTION
Large-scale atmospheric dynamics at mid-latitudes is significantly influenced by the baroclinic wave life cycle initiated by shear instability [43]. Constrained by geostrophic and hydrostatic balance, a meridional temperature gradient imposed by the heat imbalance between tropics and high latitudes is proportional to the vertical shear of the mid-latitude jet, which acts as a key parameter controlling the growth rate of baroclinic waves [8,11,41]. During the baroclinic growth of synoptic waves, poleward heat flux increases and the waves propagate upward, after which the waves break near critical latitudes while propagating equatorward. The whole life cycle of unstable baroclinic waves could be represented by energy exchange between zonal mean zonal wind and synoptic eddies [12,28,33,51]. This cycle is reflected in daily weather at mid-latitudes and understood as an engine to transfer residual energy from the tropics to high latitudes. It determines many aspects of large-scale circulations in extra-tropical areas [19].
The life cycle of baroclinic waves is an essential part of variability in mid-latitudes with its time-scale being around 3-4 days. However, the low-frequency intraseasonal variability with a time scale longer than weather but shorter than seasonal [7,16] is not entirely explained by the traditional baroclinic wave life cycle [2,44]. For example, the North Atlantic Oscillation (NAO) [22], which represents the fluctuation of the difference of atmospheric pressure between the Icelandic Low and the Azores High, contains the intra-seasonal time-scale of around 10-days. Similarly, the Arctic Oscillation (AO) * woosok.moon@gmail.com [56] or the Southern Annular Mode (SAM) [25], representing the zonally symmetric seesaw between sea level pressures in polar and temperate latitudes in Northern and Southern Hemisphere, respectively, also vary on intra-seasonal timescales.
The variability of the mid-latitude westeries is commonly defined using the zonal index [35] that was originally defined by the pressure difference between 35N and 55N, but there are many variants [24,46]. The zonal index, especially in the Southern Hemisphere, shows strong intraseasonal and interannual variability based on several observational data [17]. On intra-seasonal time-scales, it is hard to identify a major external forcing and hence the mid-latitude jet variability could be controlled by internal dynamics [7]. Indeed, the low-frequency variability of the zonal index has been identified in simple numerical models with zonally-symmetric thermal forcing [45], implying that variability can be caused by the interaction between synoptic eddies and zonal mean field.
The eddy-mean flow interaction is a highly nonlinear process involving turbulent mixing of synoptic waves leading to energy transfer among different scales that are difficult to accurately parameterize. However, phenomenological models of eddy-mean flow interactions exhibiting oscillatory behavior have been put forward. Thompson and Barnes [53] introduced a stochastic model to explain the quasi-oscillatory behavior of the poleward heat flux in the Southern Hemisphere, also referred to as the Southern Hemisphere baroclinic annular mode (BAM) introuduced by [57]. Specifically, they suggest a two-dimensional stochastic model representing interactions between the poleward heat flux and the meridional temperature gradient to capture the quasi-oscillatory variability with a dominant time-period around 25 days. In this model, the increase of the poleward heat flux is proportional to the growth rate of eddies generated by baroclinic instability, which is proportional to the meridional temperature gradient. At the same time, the time evolution of the anomalous meridional temperature gradient is controlled by the poleward heat flux with a damping. This phenomenology results in a two-dimensional stochastic dynamical system that contains oscillatory solutions depending on the choice of parameters. With this idealized description of the variability, Thompson and Barnes [53] suggest that the quasi-oscillatory lowfrequency variability can be caused by nonlinear eddymean flow interaction. Even though the model reflects the basic characteristics of the baroclinic wave life cycle, its direct connection to the primitive equations is not described, which obscures the physical interpretation of model parameters.
The BAM is a companion index of the Southern Annular mode (SAM) in the Southern Hemisphere representing the characteristics of large-scale low frequencies in the atmosphere. The SAM is defined by zonal mean kinetic energy describing the variability of zonal mean wind mainly influenced by the momentum flux. The BAM is constructed by eddy kinetic energy mainly controlled by the meridional heat flux in lower levels [57]. The two processes, meridional heat flux in lower levels and momentum flux in upper ones, show different characteristics in their variabilities [3,4,40,60,61]. The SAM is well approximated by an autoregressive model of order 1 (AR(1) process), having a red noise spectrum, but the BAM shows a distinct peak around 25 days in its power spectrum. The stochastic oscillatory behavior shown in the BAM is worthwhile to be investigated in detail due to the expectation that it leads to significant progress in predictability in mid-latitudes on sub-seasonal timescales. Furthermore, [55] argue that the BAM in the Northern Hemisphere shows a dominant peak located at around 20 days in its power spectrum. This implies that the oscillatory behavior in meridional heat flux or eddy kinetic energy is an intrinsic feature of eddy-mean flow interactions in large-scale atmospheric dynamics. Therefore, the major question is how the meridional heat flux in lower levels can induce the oscillatory behavior in the interaction with mean field.
Generation of low-frequency oscillations in a rotating fluid does not seem to be limited to the large-scale atmosphere as similar internal oscillations were also found in large-scale ocean currents. An eddy-resolving ocean model simulation with a repeated annual cycle forcing reveals an intrinsic mode in the Southern Ocean with a period of 40-50 year [23,59]. An idealized model of the surface-stress-driven Beaufort Gyre in the Arctic Ocean with mesoscale eddies as a key equilibration process [29] generates an oscillation with a period of 50 years [31]. To explain the oscillation Manucharyan et al. [31] introduce eddy memory into the relation between the eddy buoyancy flux and the mean buoyancy gradient, which leads to a modification of the commonly used Gent-McWilliams (GM) parameterization [15]. The inclusion of eddy memory leads to a stochastic oscillation if the ratio between the memory timescale and the diffusive equilibration timescale reaches a certain threshold [31].
Generally, according to the Mori-Zwanzig formalism [62], the dynamic interactions between slow-evolving modes and fast ones are reflected as a delayed integral on the time evolution of the slowly-varying modes. The memory effect, normally represented by a delayed term or integral, is not new in climate science. One of the simplest models for the El-Ninõ and Southern Oscillation (ENSO) is a delayed ordinary differential equation [52]. Even with such one-dimensional model, it was found that the delayed model contains chaotic dynamics under seasonal forcing [58]. Recently, such delay equation climate models were derived using the Mori-Zwanzig formalism [13]. Considering the complex interactions in fluid-dynamical systems, the memory effect is an intrinsic characteristic leading to various internal variability.
Our research focuses on the memory effect in the interaction between synoptic eddies and a zonal mean zonal wind and its relation to low frequency modes on timescales longer than the weather. More specifically, we propose a formalized explanation of the quasi-oscillatory behavior of the BAM [53] by incorporating the eddy memory [31] into multi-scale equations of atmospheric dynamics [32]. We outline the multi-scale atmospheric model in Section II and we introduce the eddy memory effect onto the planetary -scale equations in Section 3., where we formally derive the reduced model for the zonal mean flow, the stochastic oscillator. We conclude in Section 4.

II. MULTISCALE MODEL OF ATMOSPHERIC MOTIONS
The basic governing equation for Earth's climate starts from the simplest heat flux balance between short-wave radiative flux and outgoing longwave radiative flux, which gives us a global average temperature as an equilibrium. This is possible when we consider the Earth as a point object in the universe maintaining a thermal equilibrium state. If we magnify the size of the Earth from a point to a sphere, we can see that there is an imbalance of heat flux from tropics to higher latitudes, which requires transferring heat from the tropics to higher latitudes. A simple approximation of the poleward heat flux is a turbulent diffusion with a constant eddy diffusivity. In an even more magnified view of the tropics and polar areas, the dominant physics shifts from energy flux balance to potential vorticity conservation, which acts as a theoretical foundation for the generation of mid-latitude weather systems. There are unclearly defined timescales between the weather dominated by the vorticity dynamics and the season controlled by external heat flux. These timescales lie on around 10 to 30 days, where numerous variability has been detected in climate data. There should be an intermediate framework lying between the energy flux balance and the vorticity dynamics to investigate causes of the variability.
The large-scale atmosphere is governed by the primitive equations which are comprised of three momentum balances, the continuity equation and the heat budget. In Cartesian coordinates, the three momentum equations are where P and ρ are atmospheric pressure and density, respectively, and the velocities in the zonal, meridional, and vertical directions are u, v, and w. The Coriolis parameter f is approximated by f 0 + βy, where f 0 is a Coriolis parameter calculated at a reference latitude φ 0 and βy is the latitudinal variation of the Coriolis parameter where β = ∂f /∂y. The material derivative D/Dt is defined as D/Dt ≡ ∂/∂t + u∂/∂x + v∂/∂y + w∂/∂z. Along with the momentum equations, the continuity equation is Finally, to complete the equations, we need the ideal gas law P = ρRT and thermodynamic equation where Θ is the potential temperature, i.e., where P and T represent the pressure and the temperature, respectively, P 00 is a reference pressure, typically 1000mb, R is the gas constant for air, C P the specific heat at the constant pressure, and P 00 is a reference pressure. The full equations are highly nonlinear and complicated and to explain the BAM, we need to simplify the equations using an appropriate scaling. The scale of atmospheric motions is distinguished by the magnitude of the Rossby number, U/f L, where U is a horizontal velocity scale, and L is a horizontal length scale. When the Rossby number is asymptotically small, the geostrophic balance becomes a dominant force balance in the horizontal momentum equations. The largescale atmospheric flows satisfying the geostrophic and hydrostatic balances are called geostrophic motions.
There are two kinds of geostrophic motions depending on the horizontal length scale [42]. If the length scale is similar to the internal Rossby deformation radius (L ∼ 1000km), the motion is called quasi-geostrophic and is described by the conservation of potential vorticity, D H Dt where ψ is a leading-order pressure field acting as a stream function, ρ s is a mean vertical density profile, S represents the average vertical stability, and ∂ 2 ψ ∂x 2 + ∂ 2 ψ ∂y 2 + 1 ρs ∂ ∂z ρs S ∂ψ ∂z + βy is the potential vorticity. All the variables are nondimensionalized [39]. On the other hand, if the horizontal length-scale is close to the external Rossby deformation radius (L D ∼ 3000km), the governing equations become where P L is the planetary-scale pressure, u L , v L and w L are the zonal, the meridional and the vertical velocity, respectively, Θ s is the hemispheric average potential temperature, Θ L is an anomalous potential temperature, β L is the planetary-scale beta effect, and Q L represents radiative processes. The thermodynamic forcing Q L is a residual of local radiative fluxes driving a local temperature change. The subscript L implies planetary-scale variables. Eqs. (7)- (11) are usually referred to as planetary geostrophy, where the heat equation with advection is constrained by the geostrophic and hydrostatic balances, together with the Sverdrup relation [32]. Similar results can be found in [9,10], where the same leadingorder equations are derived. They used a different small parameter = (aΩ 2 /g) 3 instead of the Rossby number U/f L, where a is the earth's radius and Ω is the earth's rotation frequency, thus the detailed derivations toward the leading-order equations are different. On the other hand, [32] derive the same results based on the Rossby number as a small parameter in the asymptotic analysis following the historical development of theories of quasigeostrophic motions in large-scale atmospheric dynamics [39]. The two geostrophic motions coexist in the large-scale atmosphere. Hence, the large-scale atmospheric dynamics should be represented by interactions between the two scales. Because these two scales are asymptotically separate, we can use a multi-scale analysis in spatial and temporal domains. Based on the Rossby number in the planetary scale = U/(f 0 L D ), we can introduce the two scales, (X, Y,t) for the planetary scale and (x, y, t) for the quasi-geostrophic (QG) one, where (x, y, t) = 1/2 (X, Y,t). Thus, the time and spatial derivatives are scaled as The two scales are separated by 1/2 which comes from the estimation that L/L D ∼ 1/2 in the large-scale atmosphere in the earth. The internal Rossby deformation radius L is around 1000km in mid-latitudes and the external (barotropic) one L D around 3000km.
A regular perturbation analysis of the primitive equations (1-4) [32] leads to dz . The Burger number F is defined as L 2 /L 2 D where L is an internal Rossby number and L D the external Rossby number. These L and the numbers 0 and 1 are used as the subscripts to represent planetary-scale and synopticscale variables, respectively. The number 0 (the number 1) represents the leading-order (the next order) in the synoptic-scale. It is assumed that planetaryscale (synoptic-scale) variables are only dependent upon planetary-scale (synoptic-scale) coordinates, X, Y , and t (x, y, and t). The vertical coordinate z is used for the both scales. The equations (12)(13)(14)(15)(16) show the dynamics of the two scales and their mutual interactions. The detailed derivation and discussions are found in [32]. The equations (12)(13)(14) describe that the two scales satisfy the basic balances, geostrophic and hydrostatic balances. The continuity equation (14) in the planetary scale includes the O(1) beta effect, which is known as the Sverdrup relation. The two scales contribute to the energy flux balance in the heat equation (15) with their own temporal and spatial scales. The last equation (16) is the vorticity equation governing the quasi-geostrophic dynamics. Here the planetary geostrophic motion provides a mean field for the development of vorticity on the Rossby deformation scale. [10] also considered a similar multi-scale analysis to represent the interactions between planetary and synoptic scales, which leads to equivalent results. However, their main focus lies on vorticity dynamics at the two scales, and hence a planetary vorticity equation is derived from the combination of the heat equation and the Sverdrup relation along with the quasigeostrophic potential vorticity equation. Our focus is in the planetary-scale heat equation with thermal forcing Q L . Lying between the simple energy flux balance model [37] and the quasi-geostrophic dynamics, the planetary heat equation together with basic dynamic balances connects the planetary thermal heat flux to fluid dynamics on planetary scales.
The above equations become simpler when the planetary-scale thermodynamic forcing Q L is zonally homogeneous. The planetary scale preserves that symmetry, in which case all X-derivative terms vanish. Hence, the zonal means of the planetary meridional and vertical velocities v L and w L become zero. This yields We can introduce a planetary-scale time-and spatialaverage to the QG variables under the assumption that the time-and spatial-average of a QG variable is close to zero. This implies that the overall effect of synoptic scales on the planetary-scale motions is represented by the temporal and spatial average of QG-scale fluxes. In particular, it is important to consider a planetary-scale spatial-average over the terms containing the QG-scale spatial derivative such as ∂/∂x and ∂/∂y. Let's consider a large horizontal length l for the spatial-average in the synoptic-scale and then the meridional average of ∂F ∂y Y0 , where Y 0 represents a position in the planetary-scale coordinate, is then interpreted as Here l and ∆Y are non-dimensional lengths in synoptic-and planetary-scale, respectively. The same length is expressed in two length scales, i.e., l * = ∆Y * where * is used to represent dimensional quantities, which is same as Ll = L D ∆Y . Thus, l = ∆Y / 1/2 where L/L D = 1/2 is used. Due to the asymptotic difference between the two scales, the length l in the synoptic-scale is approximated by the limit of the planetary-scale length ∆Y toward zero. The heat equation (18) and the QG vorticity equation (19), after taking the planetary-scale time-average ( · ), become where the time average of linear terms in synoptic scales such as ∂Θ 0 /∂t and ∂Θ 0 /∂x is assumed to be zero and the nonlinear terms representing heat and momentum fluxes are considered as non-zero average terms. Using ∇ 2 P 0 = ∂v 0 /∂x − ∂u 0 /∂y, we find that the equation (22) is equivalent to Now, we can take the planetary-scale spatial average over the equations (21,22), and combine the two by replacing w 1 by the horizontal vorticity convergence, which leads to The equation (24), after taking the zonal-average ( · ) , becomes where the dominant balance after ignoring the momentum flux contribution in the O( ) order is In the equation (26), we can consider an asymptotic solution Θ L [Θ L ] + 1/2 η, in which case the radiative process Q L is represented approximately as This is understood as a local radiative energy flux balance, and [Θ L ] represents the seasonal mean of the potential temperature mainly determined by the seasonallyvarying radiative flux balance [Q L ]. The simplest rep- where S 0 is the shortwave radiative flux, α is the local albedo, and σ[Θ L ] 4 is the outgoing longwave radiative flux, but we have to consider other heat fluxes such as incoming longwave radiative flux and turbulent sensible and latent heat flux. In the planetary-scale heat equation, the leadingorder governing physics is the local energy flux balance contributing mainly to a seasonal cycle of the potential temperature.
The O( 1/2 ) order balance is The fluctuation η around the seasonal cycle [Θ L ] is controlled by the seasonal sensitivity of the radiative processes a(Y,t) and the meridional heat flux convergence induced by synoptic eddies. [34] focus on the monthlyaverage variability of surface air temperature based on a periodic non-autonomous stochastic differential equationη = a(t)η + N (t)ξ(t), where a(t) is equivalent to the a(Y,t) in (28), ξ(t) a white noise mimicking the overall effect of weather-related processes and N (t) is a monthlyvarying amplitude of noise. The noise forcing N (t)ξ(t) can be understood as an approximation of the residual forcing R(Y,t) which can be considered as a contribution of short-time processes in the equation (28). The monthly statistics including variance and lagged correlations in a given monthly-averaged data such as surface sea temperature or tropical climate indexes is regenerated by the periodic non-autonomous stochastic model with an appropriate choice of the two periodic functions a(t) and N (t). In particular, the positive sign of the a(t) implies existence of positive feedbacks which magnifies the magnitude of a given perturbation. For example, when the stochastic model is applied to the Nino 3.4 index, the a(t) is positive from July to November, which shows the seasonality of the Bjerknes feedback. The construction of the monthly sensitivity a(t) enables us to detect how a positive feedback shapes the monthly statistics of a climate variable. It is essential to understand how phase locking and seasonal predictability barrier of climate phenomena are associated with the magnitude and timing of positive feedback. The equation (28) is an extension of the onedimensional stochastic model considering meridional variation. Especially, it represents a strong influence from synoptic-scale eddies on the planetary-scale temperature. It is in contrast with [3] which provide a similar budget equation on the planetary scale based on the multi-scale analysis suggested by [10]. In their budget equations, the planetary scale is independent from the synoptic scale, thus the two scales interact indirectly through source terms. On the other hand, the equation (28) tells that the variability of planetary-scale temperature is mainly determined by the turbulent heat flux induced by synoptic eddies. One remaining step to close the planetary-scale equations is to parameterize the meridional heat flux v 0 Θ 0 using the planetary-scale temperature Θ L .

III. EMERGENCE OF A GENERALIZED LANGEVIN EQUATION
A. Fickian diffusion model The meridional heat flux at mid-latitudes under the growth of synoptic waves is the major consequence of baroclinic instability. The baroclinic instability is initiated from an unstable mean state measured by the vertical shear of zonal mean wind ∂u L /∂z. The growth of baroclinic waves by a baroclinic instability starts from near surface, which induces a meridional heat flux. Hence, the meridional gradient of the planetary temperature ∂Θ L /∂Y is strongly related to the meridional heat The simplest representation of the mutual relationship between ∂Θ L /∂Y and v 0 Θ 0 is a turbulent flux gradient parameterization, v 0 Θ 0 −K∂Θ L /∂Y , where K is a constant [37]. If we use this parameterization in the equation (26), this leads to This simple energy flux balance model was first introduced by [50] and [36] to include the effect of large-scale atmospheric dynamics upon the local energy flux balance. The consequence of complicated large-scale atmospheric dynamics especially in mid-latitudes is to transfer the surplus of energy in low latitudes to high latitudes, which is simply approximated by a turbulent meridional diffusion. The temporal and spatial evolution of the perturbation η is represented by where we include the R(Y,t) for the contribution of shorttime processes. We can consider two boundary conditions in meridional coordinate, ∂η ∂Y (·, Y = 0) = ∂η ∂Y (·, Y = 1) = 0 implying there is no meridional heat flux near the pole and the tropical areas.
For the Y -direction diffusion operator, consider the eigenvalue problem, with ∂Hn ∂Y (Y = 0) = ∂Hn ∂Y (Y = 1) = 0. Here, H n = A n cos(nπY ) with λ n = Kn 2 π 2 and n = 0, 1, 2, · · · . Hence, the temporal and spatial perturbation of the potential temperature anomaly η can be represented by an infinite series of the eigenfunctions with time-varying coefficients, Consider first the case that a(Y,t) is a constant −γ and that the sum of short-time processes R(Y,t) is represented as a sum of the same eigenfunctions, R = ∞ n=0 R n (t)H n (Y ). The time-varying coefficient x n (t) then satisfies If we consider the last term as a stochastic noise, in particular, a Gaussian white noise, this becomes a Langevin equation with a decay time scale 1/(λ n + γ). The time scale becomes shorter as n increases, which indicates that the first several modes representing large-scale motions dominate in the overall fluctuations. When n is equal to 0, the Langevin equation becomesẋ 0 = −γx 0 + R 0 (t), where γ represents the seasonal sensitivity of local energy flux balance mainly dominated by radiative-convective equilibrium in land and ocean heat flux in ocean boundary layer. The formalism used in this argument was first introduced by [30] to explain the internal variability caused by interactions between meso-scale eddies and mean fields under the framework of the Gent-McWilliam parameterization. The simple Langevin equation was introduced to climate science by [18] to interpret climate variability in terms of stochastic dynamics. The deterministic part is understood as a stabilizing process around a mean state and the stochastic noise as the effect of weather. Climate is understood as a combination of the stability of slowlyevolving backgrounds and short-time processes approximated as a noise. The interpretation was useful to rationalize the ubiquitous emergence of red noise spectra from climate variables including sea surface temperature (SST) [14]. In a symmetric hemisphere, the meridional variation of mean near-surface temperature could be constructed by the equation (29) and the fluctuation around the mean be represented by the equation (33), which might be consistent with the usage of an AR1 process explaining the zonal index variability [27].
B. Quasi-oscillatory behavior of the BAM [53] (TB) studied quasi-periodic variability of the Southern Hemisphere baroclinic annular mode (BAM). Instead of the red noise spectrum originated from a simple linear Langevin equation, the BAM shows a clear sign of quasi-oscillation in its power spectrum. This indicates that the simple Langevin equation leading to a red-noise spectrum is not adequate to describe the quasi-periodic fluctuation.
To explain a periodic behavior of meridional temperature gradient and poleward heat flux, TB introduce mutual interaction and feedback between the meridional temperature gradient b ≡ ∂η/∂Y and the poleward heat flux v 0 Θ 0 . Baroclinic instability theory tells us that the growth rate of the baroclinic waves is proportional to the meridional temperature gradient [26], which could be represented by where α is a constant representing the amplitude of the feedback between the baroclinicity and the meridional heat flux [20]. Here, the noise forcing (t) is added to include the effect of chaotic weather. At the same time, TB consider a feedback of the poleward heat flux on the meridional temperature gradient, which simply assumes that the meridional temperature gradient b increases linearly with respect to the poleward heat flux v 0 Θ 0 . Thus, the second equation in TB is where µ is the degree of the feedback and r is a recovering time-scale of the meridional temperature gradient due to diabatic processes and vertical motions. The combination of (34-35) could generate an oscillation whose period is determined by the relevant coefficients. While ignoring the noise forcing in the equation (34), combining the two equations [53] leads to where the G is the time-derivative of the short-time scale forcing (t). Oscillatory behavior comes out when 1 4r 2 − αµ < 0. TB used r 4 days and αµ 0.05days −2 to generate a dominant peak of the power spectrum of poleward heat flux around 20 days. They suggest a mutual interaction between v 0 Θ 0 and ∂η/∂Y with the damping represented by a Newtonian cooling. In the model, it is hard to guess what determines the damping time-scale r and the feedback parameter αµ. By a simple comparison, the equation (35) is equivalent to the anomalous heat equation (28) and then the equation (34) provides a relationship between the poleward heat flux and the meridional temperature gradient. Even though the equation (35) is likely to be derived from the anomalous heat equation, the damping term −b/r is not clearly understood physically. To understand the relationship between b and v 0 Θ 0 on intra-seasonal timescales, it seems necessary what determines the damping time scale r.
We can incorporate the equation (34) into the equation (28) after taking time-derivative on the equation (28). This results in where −γ is used instead of the a(Y,t) for simplicity. The constant γ implies the damping time scale of seasonal energy flux balance. In land, the outgoing longwave radiative flux dominantly controls the time scale. Similarly, we can consider the eigenvalue problem for the diffusion operator, where η = ∞ n=0 x n (t)H n (Y ). We obtain the timeevolution equation for x n (t), where R = ∞ n=0 R n (t)H n (Y ) is used. The equation for x n is similar to the equation (36), which suggests that the damping time scale is equivalent to 1/γ. Physically, the γ is introduced as a sensitivity of the radiative energy flux balance. In terms of time-scale, it could be interpreted as a time-scale to return to a climatological seasonal cycle. [34] developed a time-series method to construct the γ from monthly-averaged surface air temperature. The time scale of the surface energy flux balance in the Southern Hemisphere is around 1.5 months, which is much larger than 4 days. Therefore, the damping time scale r does not come from diabatic processes related to radiative fluxes.

C. Eddy memory and generalized Langevin equation
The damping time scale introduced to explain the oscillatory behavior of the baroclinic mode in Southern Hemisphere is much shorter than that of the mean seasonal cycle of radiative processes. It is plausible that the time scale is related to synoptic eddies, rather than any external influences which have longer time-scales. Synoptic eddies generated from the baroclinic instability due to an unstable background undergoes a specific energy cycle with a zonal mean steady state. The baroclinic instability enables the synoptic eddies to extract energy from the zonal mean state, from which the poleward heat flux increases near surface. The growing waves propagate upward and equatorward and meet critical latitudes where phase speeds of waves are same as mean wind. The waves break and give their energy back to the mean state by momentum flux. This baroclinic wave life cycle is complete in a few days.
The overall effect of the baroclinic wave life cycle on the meridional heat flux is represented by the planetary potential temperature anomaly η as v 0 Θ 0 = −K ∂η ∂Y . This Fickian diffusion approximation is well applied when the time evolution of the planetary variable is much slower than the turn-over time-scale of synoptic eddies. On intraseasonal time scales, however, the time scale for the evolution of planetary-scale variables is not clearly separated from that for the baroclinic wave life cycle. The advective timescale for the synoptic waves is around 1.2 days (where we use L = 1000km and U = 10m/s) and the same timescale for the planetary waves around 3.5 days (where we use L D = 3000km and U = 10m/s). A complete cycle of baroclinic wave lifecycle takes around 3-4 days, hence the timescale for the planetary-scale motion is comparable with that for the baroclinic wave life cycles. It is questionable to apply the Fickian diffusion as a parameterization of the poleward heat flux in these time scales. Non-Fickian approximation for turbulent transport or diffusion has been a central topic in turbulent closure problems [38]. For a transient and chaotic turbulence, the Fickian approximation is not enough to capture nonlocal and non-Gaussian nature of turbulence. One of the simplest approach is called the minimal τ approximation, where the third-order momentum represented as a forcing for the time-evolution of second-order moments is approximated as a damping term with the timescale τ [5,6]. This is equivalent to apply an integral kernel instead of a constant diffusivity with a finite memory [21].
The main idea is applied to the time-evolution of the poleward heat flux of synoptic-scale waves, i.e., The time-evolution equations for the synoptic-scale meridional velocity v 0 and potential temperature Θ 0 [32] are Hence, where D contains the terms representing the contribution from synoptic eddies. After taking synoptic-time average on both sides, we obtain The main assumption of the minimal τ approximation is that the contribution of the synoptic eddies D is approximated by a damping of the v 0 Θ 0 with the designated time-scale r. Therefore, where K ≡ v 2 0 suggesting how the eddy diffusivity K is related to the second-order statistics of synoptic eddies. The minimal τ approximation for the v 0 Θ 0 seems to be equivalent to the equation (34) in the TB. The feedback parameter α could be understood as K. TB suggests the relationship based on the result from the baroclinic instability. On the other hand, we derived a similar one by the simplest non-Fickian approximation as a closure of turbulent eddies. TB included a random forcing in the relationship to include unresolved processes, but the randomness coming from short-time processes is considered in the heat equation in our derivation. The above parameterization can be represented by an integral form, i.e., The integral form originates from integrating (44) for the v 0 Θ 0 with respect to the timet. From this integral form, we can see that the poleward heat flux is the result of accumulating the baroclinicity during a certain time until present. It represents the memory effect caused by synoptic eddies.
The minimal τ approximation was tested in direct simulations of isotropic 3D turbulence [5]. The statistics of a passive scalar are compared between direct simulations and the parameterized equations, where decent matches are obtained. In particular, the parameterization transforms the main parabolic tracer equation to a wave equation leading to an oscillatory behavior. The simulations show decayed oscillation with a certain choice of the damping time-scale r. In an idealized Beaufort Gyre numerical simulation, [31] introduced a finite memory kernel instead of a constant diffusivity in the Gent-McWilliam parameterization, which generates the quasiperiodic variability in the eddy-mean interaction. The minimal τ approximation can be understood as a finite memory effect represented by an integral kernel.
Taking account of memory effect of eddies, we can introduce an integral kernel κ(t − t ) on the parameterization of the meridional heat flux, which is a generalization of the minimal τ approximation in the equation (45) where η * is an effective temperature anomaly defined by Thus, the equation (30) becomes where we assume that a(Y,t) = −γ for simplicity. Furthermore, if we expand the η * using the eigenfunctions of the diffusion operator, η * = ∞ n=0 x * n (t)Q n (Y ), and the residual forcing, R = ∞ n=0 R n (t)Q n (Y ), each timedependent coefficient x n satisfies which is a generalized Langevin equation [63]. The generalized Langevin equation contains memory terms representing the effect of past states, which is generally originating from interactions with short-time scale components.
The complicated chaotic processes to generate memory are captured by the memory kernel κ(t −t ). The dynamics represented by a generalized Langevin equation is dependent upon the choice of integral kernel. We can follow the theme of the τ approximation, using an integral kernel with a finite-time memory. It has been suggested that baroclinic eddies generated by baroclinic instability are eventually organized to maintain a midlatitude jet characterized by meridional temperature gradient [27,47]. This implies that the effect of baroclinic eddies upon the meridional heat flux is limited within a finite time-scale, which leads to as an approximation Here r represent the finite memory of the baroclinic eddies, which is equivalent to dx * n /dt = −x * n /r + x n /r. Based on the specific memory kernel, the equation (49) becomes There are two-damping time-scales in the above equation defined by r and 1/γ. The 1/γ comes from the seasonal heat flux balance such that it is approximately 1.5 months in Southern Hemisphere. Thus, it follows 1/r γ. The eigenvalues λ n come again from (31) with Q n (0) = Q n (1) = 0 implying no heat flux at both ends, which gives us λ n = Kn 2 π 2 and Q n (Y ) = cos(nπY ). The baroclinic annular mode (BAM) is defined by a dominant mode of the meridional heat flux maximized around the center of mid-latitudes. This is similar to the mode with n = 1, Q 1 (Y ), whose derivative Q 1 (Y ), has its maximum at the center. The time-dependent coefficient of the mode, x 1 , satisfies which should be understood as an equivalent form to the equation (36). The damping time-scale r = 4 days comes from the memory effect of synoptic eddies on meridional heat flux, which seems to be strongly related to the baroclinic wave life cycle. Moreover, αµ in the equation (36) is same as (γ + Kπ 2 )/r, which leads to K 1/5π 2 . Depending on the choice of the parameters r and K, the equation (52) can show decayed oscillatory behavior or exponential decay. Under the assumption that R 1 /r + dR 1 /dt is approximately a Gaussian white noise, we can simulate the equation to generate a stochastic realization. Figure 1 shows several power spectra of stochastic realizations from the equation (52) with different parameters r and K along with the power spectrum of the BAM index [57]. With a fixed K, four different memory time-scales from 1 day to 4 days are used for stochastic realization. In figure 1(a), we see that longer memory (4 days) leads to oscillatory behavior, which has the similar power spectrum with that of the BAM index, and shorter memory (1 day) to a red noise spectrum (Fig. 1a). Similarly, we fix r = 4 days and vary K, which also shows the characteristics from stochastic oscillation to red noise process (Fig.1b). Similarly, when Kπ 2 /r = 1 days −2 , the power spectrum is close to that of the BAM index.
The baroclinic wave life cycle results in recovering the vertical shear of the mean jet. This could be an origin of the finite memory effect in meridional heat transfer in planetary-scale, which leads to oscillatory behavior of zonal mean potential temperature. Following the above formalism, the fluctuation of potential temperature anomaly η is represented as η = ∞ n=0 x n (t)Q n (Y ), hence ∂η/∂Y = ∞ n=0 x n (t)dQ n /dY . By the thermal wind balance, ∂U/∂z ∝ ∞ n=0 x n (t)dQ n /dY , which implies that the vertical shear of the jet in lower levels could show quasi-oscillatory behavior depending on the timescale of the eddy memory and relevant eigenvalues defining the decay time-scale of a specific normal mode.
The main focus lies on how to parameterize the meridional heat flux induced by synoptic eddies v 0 Θ 0 using the planetary-scale potential temperature Θ L . Considering that the time scale for synoptic eddies is not much smaller than that of planetary-scale motions, the meridional heat flux is not entirely determined by an instantaneous potential temperature gradient; instead it is dependent upon the temporal history of the potential temperature gradient. Hence, the Fickian diffusion approximation is not appropriate to parameterize the meridional heat flux induced by synoptic eddies. The simplest way considering the temporal history of the planetary-scale potential temperature is to introduce a time integration of the potential temperature with an appropriate memory kernel, which enables to include the past influence of planetary-scale mean fields on the turbulent heat flux. As a result, fluctuations of potential temperature can be described approximately by a generalized Langevin equation containing a memory term. The time-delayed effect represented as an integral kernel in the generalized Langevin equation can generate various types of variability ranging from simple red-noise process to quasioscillations and chaos. Hence, an appropriate choice of integral kernel is crucial to define the statistical characteristics of planetary-scale variability.

IV. CONCLUSION
The main goal of the study was to rationalize the quasioscillatory variability of the BAM [53] on intra-seasonal time scales using the multi-scale representation of primitive equations of the large-scale atmosphere flow [32]. Conceptually, the multi-scale approximation is a combination of the energy flux balance in the heat equation and the potential vorticity conservation. We found that the critical issue under the framework is the specification of an appropriate parameterization of the horizontal heat flux convergence from the synoptic and planetary-scale eddy dynamics. If the heat fluxes are represented using the simplest diffusive (Fickian) parameterization with a constant eddy diffusivity, the planetary-scale equations simplify to an energy flux balance with a meridional turbulent diffusion, similarly to [36]. Other types of eddy parameterizations are also possible, of those most notable is the residual mean theory of Andrews and Mcintyre [1] with a diffusive parameterization of isopycnal eddy fluxes of potential vorticity. While the residual mean formulation gives a different perspective on the mean flow dynamics compared to the direct diffusive closure of the heat fluxes [49], our study questions their equilibrium-type Fickian assumption for the relation between the eddy fluxes and gradients of the considered tracer fields. Under the equilibrium-type eddy parameterizations, the main equation in the planetary-scales atmosphere is the heat equation with the basic balances resulting in a parabolic partial differential equation that does not support natural oscillations. The resulting temporal variability of observables obeys the Langevine equation leading to a red-noise spectrum that does not ex-plain the quasi-oscillatory behavior of the BAM. Thus the equilibrium-type eddy parameterizations are likely not appropriate at the intraseasonal timescales and can only be used for representing the relatively long-term average of the planetary-scale dynamics.
We find that using a non-equilibrium eddy parameterization with the eddy memory effect represented by the delayed integral in the flux-gradient relation as was proposed for mesoscale ocean eddies [31] can explain the oscillatory behavior of planetary-scale atmospheric flows. The rationalization for the eddy memory effect comes from the notion that the overall effect of synoptic eddies is to maintain the vertical shear of a mid-latitude jet [47], implying that the eddy energy cycle starting from baroclinic instability near the surface could be considered as a finite timescale process. Thus, we hypothesized that the mid-latitude eddies have a finite memory and introduce a simple exponentially decaying memory kernel into the relation between the poleward eddy heat flux and the meridional potential temperature gradient. The inclusion of the finite memory kernel directly converts the heat equation from parabolic to wave equation that allows dispersive planetary-scale waves. As a result, instead of the red-noise model, the temporal variability of the anomalous potential temperature near the surface now obeys a stochastic oscillator model, with spectral characteristics that are similar to the observables associated with the BAM.
The key value of our theoretical study is that it highlights the processes and clarifies the physical interpretation of the crucial parameters that may be governing the spectral characteristics of the BAM. Specifically, our study implies that the period of natural oscillations of the BAM is proportional to the geometric mean between the eddy memory and the diffusive equilibration timescales. The ratio of these timescales controls the existence of a pronounced spectral peak while the variance is directly proportional to the external noise variance in the heat equation. Validating our key hypothesis about the exponential memory kernel and the finite eddy memory timescale in the relation between the mean flow and eddy heat fluxes using atmospheric models would be the next crucial step towards validating our theoretical arguments about the BAM. Furthermore, understanding the nature of the eddy memory and how its timescale depends on the mean flow or eddy characteristics would be necessary to understand the climate conditions under which the BAM could exhibit oscillatory behaviour. Finally, our results emphasize the importance of the external noise in driving the variance of the BAM and the necessity to understand if it is dependent on the mean flow or acts as an external/independent source of energy.