We propose a new reduced fluid model for the study of the drift wave–zonal flow dynamics in magnetically confined plasmas. Our model can be viewed as an extension of the classic Hasegawa-Wakatani (HW) model and is based on an improved treatment of the electron dynamics parallel to the field lines, to guarantee a balanced electron flux on the magnetic surfaces. Our flux-balanced HW (bHW) model contains the same drift-wave instability as previous HW models, but unlike these models, it converges exactly to the modified Hasegawa-Mima model in the collisionless limit. We rely on direct numerical simulations to illustrate some of the key features of the bHW model, such as the enhanced variability in the turbulent fluctuations and the existence of stronger and more turbulent zonal jets than the jets observed in other HW models, especially for high plasma resistivity. Our simulations also highlight the crucial role of the feedback of the third-order statistical moments in achieving a statistical equilibrium with strong zonal structures. Finally, we investigate the changes in the observed dynamics when more general dissipation effects are included and, in particular, when we include the reduced model for ion Landau damping originally proposed by Wakatani and Hasegawa.

In toroidal magnetically confined plasmas, the level of cross field heat and particle transport in the edge region is in large part set by drift-wave turbulence driven by temperature and density gradients.17 This turbulence can itself generate zonal flows,8 which are known to mediate turbulent transport by shearing the turbulent eddies and absorbing some of the drift wave turbulence energy.3,6,8,11,18,28,32 Understanding this process is a critical step towards the goals of minimizing heat transport in magnetically confined plasmas and designing compact, economically viable fusion devices.

A wide variety of models have been used to study the drift wave turbulence–zonal flow dynamics in toroidal plasmas, with a varying degree of physics fidelity;1,11,14,15,21,28,29,32,33 interested readers may find a longer list of relevant references in Ref. 8. In the present work, we propose a new fluid model which continues the rich tradition of studies based on reduced fluid models1,7,13,16,21–24,30,31 initiated by the pioneering works of Hasegawa and Mima14 and Hasegawa and Wakatani.15 The approximations made to derive many of the reduced models, including ours, are often only satisfied for parameter regimes which do not correspond to the values measured in the edge of magnetic confinement fusion plasma experiments. We therefore do not expect numerical results obtained with reduced fluid models to be in quantitative agreement with measurements. However, because of their relative simplicity, reduced models have played an important role in characterizing the mechanisms responsible for the formation and dynamics of zonal flows and continue to be important tools to provide more insights into these important phenomena. In order to put our new model in context, we propose a very brief review of the main reduced models considered in the plasma physics literature.

The Hasegawa-Mima (HM)14 model is a one-field partial differential equation (PDE) which is the simplest known model containing the drift wave turbulence–zonal flow feedback loop mechanism. It was recognized in the 1990s10,11 that in the context of magnetic confinement fusion, it is more physically relevant to consider a corrected form of the HM equations, now referred to as the modified Hasegawa-Mima (mHM),7 in which the magnetic surface averaged electron density response is subtracted from the original HM (oHM) equation.10 The mHM model has the desirable feature of being Galilean invariant under boosts in the poloidal direction7 and is known to lead to a stronger generation of zonal flows. This modification is also at the heart of the new model we propose.

A limitation of the HM models is that they do not contain a natural instability, so that turbulent forcing must be added externally in order to study the turbulence–zonal flow dynamics. In contrast, Hasegawa and Wakatani15 have shown that a generalized version of the HM models which includes electron ion friction in the parallel direction naturally contains a drift instability due to the finite plasma resistivity and thus drift wave turbulence induced transport. This is the original Hasegawa-Wakatani (oHW) model, a system of coupled PDEs for the ion vorticity (which is the Laplacian of the electrostatic potential) and the ion density. However, the oHW model has shortcomings too. Just like the oHM, it is not Galilean invariant. Furthermore, zonal flows are not generated in the oHW model,21 unless one makes the adiabaticity parameter wave number dependent.24 Numata et al. have proposed a modified Hasegawa-Wakatani (mHW) model obtained by subtracting the zonal components from the resistive coupling term. The mHW model addresses the issues of the oHW mentioned above: its drift wave turbulence can lead to the strong generation of zonal flows, and it can be shown that the model has the desired Galilean invariance. Even if so, the mHW has the weakness that it does not converge to either the oHM model or the mHM model in the limit of zero resistivity and zero dissipation, as we will demonstrate numerically in this article.

We propose a generalized Hasegawa-Wakatani model with an improved treatment of the electron response along the magnetic field lines. Our model addresses the aforementioned limitation of the mHW model by solving for the mHM potential vorticity instead of the vorticity and reduces to the mHM model in the adiabatic, nondissipative limit, as desired, while maintaining the Galilean invariance of the mHW model. We call our model the flux-balanced Hasegawa-Wakatani model (bHW) to emphasize the fact that in this model, the surface-averaged electron density is constant in time when the electrons are adiabatic, as one would physically expect.7,10,11 Our simple modification to the mHW model leads to major differences in the observed dynamics. Most significantly, the generation of zonal flows is enhanced and the turbulent fluctuations about the zonal mean state are increased. In our new Hasegawa-Wakatani model, we have also chosen to generalize the form of dissipative effects introduced in the oHW model to represent a wider variety of physical processes. As we acknowledged previously, we do not expect quantitative agreement between our fluid model and measurements from the edge of magnetic confinement fusion plasma experiments, but we hope to improve our understanding of the interplay between competing physical effects and of the implications of making simplifying fluid assumptions to model subtle kinetic effects. As an illustration, in this article we make the observation that the simplified linear Landau damping term of the oHW model15 which we also included in our model, acting mostly on the largest scale modes, may not always act as a pure dissipative effect and can effectively increase the variability in the flow fluctuations. To help the reader throughout the article, we provide a summary of the similarities and differences between the flux-balanced Hasegawa-Wakatani model and the modified Hasegawa-Wakatani model in Table I.

TABLE I.

Summary of the similarities and differences between the mHW model and the bHW model observed in the results of direct numerical simulations.

Flux-balanced Hasegawa-Wakatani (bHW) modelModified Hasegawa-Wakatani (mHW) model
Zonal jet structure (Figs. 5 and 6Stronger and more turbulent zonal jets are generated Jets are less persistent, but more steady when present, with weaker fluctuations 
Hydrodynamic limit α0 (Figs. 2, 3, and 5The bHW model maintains zonal jet structures, with a reduced particle flux ũñ¯ The mHW model reduces to fully homogeneous turbulence with strong particle flux 
 Stronger variability in the entire variance spectrum and especially in zonal modes with wavenumber ky = 0 Weaker variability in the variance spectrum and no anisotropy in the zonal direction 
Adiabatic limit α (Figs. 2 and 4The bHW model converges uniformly to the mHM model with little small-scale fluctuation The mHW model saturates in a final state different from the mHM final state, with intermittent small-scale vortices 
Third-order moment statistics (Fig. 3The mHW and bHW models both have important third-order moment feedbacks to the statistical equations for the mean and variance  
 Stronger third-order moment feedback, especially along zonal modes Weaker zonal feedback in the third-order moments 
Flux-balanced Hasegawa-Wakatani (bHW) modelModified Hasegawa-Wakatani (mHW) model
Zonal jet structure (Figs. 5 and 6Stronger and more turbulent zonal jets are generated Jets are less persistent, but more steady when present, with weaker fluctuations 
Hydrodynamic limit α0 (Figs. 2, 3, and 5The bHW model maintains zonal jet structures, with a reduced particle flux ũñ¯ The mHW model reduces to fully homogeneous turbulence with strong particle flux 
 Stronger variability in the entire variance spectrum and especially in zonal modes with wavenumber ky = 0 Weaker variability in the variance spectrum and no anisotropy in the zonal direction 
Adiabatic limit α (Figs. 2 and 4The bHW model converges uniformly to the mHM model with little small-scale fluctuation The mHW model saturates in a final state different from the mHM final state, with intermittent small-scale vortices 
Third-order moment statistics (Fig. 3The mHW and bHW models both have important third-order moment feedbacks to the statistical equations for the mean and variance  
 Stronger third-order moment feedback, especially along zonal modes Weaker zonal feedback in the third-order moments 

The structure of this paper is as follows. We first briefly review the classical Hasegawa-Mima and Hasegawa-Wakatani models in Sec. II, in order to provide a historical context for our modifications to the Hasegawa-Wakatani model. In Sec. III, we present our new Hasegawa-Wakatani model and derive its conserved quantities. In Secs. IV and V, we study the main differences between the mHW model and our model using numerical simulations; Sec. IV focuses on statistical considerations while Sec. V illustrates our conclusions and observations from Sec. IV with snapshots of results of direct numerical simulations and explores the role of the different dissipation terms. We summarize our work in Sec. VI and present the main linear stability properties of our new balanced model in  Appendix A.

In this section, we review the central features of the classical Hasegawa-Mima (HM) and Hasegawa-Wakatani (HW) models, as well as their modified versions,10,21 known to excite more realistic and stronger drift-wave–zonal flow dynamics. This short review is meant to provide the motivation for our new model. The presentation is purposefully brief as clear, longer, and more detailed presentations are readily available in the literature.7,8,17,21

For simplicity, all the models we will discuss in this article are considered for a shearless slab geometry,2,7 in which the toroidal magnetic surfaces are imagined to be flattened into planes parallel to the y and z axes, as shown in Fig. 1, where (x, y, z) is the Cartesian coordinate system used to describe the geometry, x representing the radial distance, which can be viewed as a flux surface label, and y and z playing the roles of the poloidal and toroidal angles, respectively. The magnetic field is assumed to be solely in the z direction, B=B0z, and B0 is constant and uniform. The equilibrium density depends on the radial variable, n0(x), and we will treat the density profile within the framework of the standard “local approximation,”27 in which n0(x)/n0(x)=constant.21,31 The electron temperature is uniform throughout the plasma, and the ion temperature is assumed to be small compared to the electron temperature: Ti/Te1. In all the models discussed in this article, the physical quantities will be assumed to be uniform in the z-direction, corresponding to the fact that in strongly magnetized plasmas, the dynamics is highly anisotropic, with much slower variations of the physical quantities along the magnetic field than across the magnetic field. The problem is therefore reduced to a two-dimensional problem in which all quantities only depend on x, y, and the time t, and the computational domain Ω is a rectangle whose sides have lengths Lx and Ly. Finally, the boundary conditions on the perturbed quantities (n, φ), which are the quantities the models solve for, are periodic in both x and y.31 

FIG. 1.

Nested toroidal flux surfaces in a tokamak geometry. The subplot (a) in the bottom left corner gives an illustration of the slab approximation for the plasma edge, in which the curved flux surfaces are flattened into planes, and shows the coordinate system used in this article.

FIG. 1.

Nested toroidal flux surfaces in a tokamak geometry. The subplot (a) in the bottom left corner gives an illustration of the slab approximation for the plasma edge, in which the curved flux surfaces are flattened into planes, and shows the coordinate system used in this article.

Close modal

The Hasegawa-Mima models, both original and modified, are equations for the ion vorticity obtained by combining the ion continuity equation and the ion momentum equation, together with the assumption of an adiabatic electron response to close the system of equations.7,14 The only difference between the oHM and mHM models lies in the treatment of that adiabatic electron response. More than a decade after the introduction of the oHM model, it was indeed recognized that one should subtract the zonal mean of the electrostatic potential in the equation relating the adiabatic electron density and the electrostatic potential in order to prevent unphysical net radial transport of electrons.7,9,10,24

The oHM and mHM equations can be formulated under the same framework by defining a switch parameter j with j =0 for the oHM model and j =1 for the mHM model, which accounts for the different treatment of the adiabatic electrons. The unified equation is13,17

qt+J(φ,q)κφ̃y=0,q=2φ(φ̃+δj0φ¯).
(1)

In Eq. (1), J(φ,q)=xφyqyφxq is the Jacobian associated with the advection term vE·q, where vE is the E×B velocity, t is the time normalized to the ion cyclotron frequency ωci=eB0/m, and x and y are normalized in terms of the hybrid ion thermal Larmor radius ρs=ωci1(Te/mi)1/2=miTe/eB0,δj0 is the Kronecker delta, which is equal to 1 if j =0 and 0 otherwise, q=2φ(φ̃+δj0φ¯) is the potential vorticity, with φ=eϕ/Te the normalized electrostatic potential, where e is the charge of the electron and ϕ the electrostatic potential, and κ=dlnn0/dx. A bar over a quantity f represents the zonally averaged mean of that quantity, which only depends on x, and a tilde represents the fluctuation component of f, obtained by removing the mean from f

f¯(x)=1Lyf(x,y)dy,f̃=ff¯,

where Ly is the length of the domain in the y direction.

The mHM modification to the oHM model is simple, only appearing in the definition of the potential vorticity q. Yet it has important physical implications. First, unlike the oHM model, the mHM model is Galilean invariant under boosts in the poloidal direction,7 which is a desired property for our shearless slab geometry. Second, it can be shown, using linear theory, that in the absence of mean flow, the drift wave dispersion relation is identical in both models and given by ω=kyκ/(1+k2), where ky is the y component of the wavevector k, and k the magnitude of k. However, in the presence of a constant and uniform background mean flow in the y direction, vE=v¯ŷ, the dispersion relation is modified in different ways in the two models7 

oHM:ω=kyκ1+k2+k21+k2kyv¯,mHM:ω=kyκ1+k2+kyv¯.

We see that at small scales, k1, the dispersion relations of the oHM and mHM models agree and are given by ω=ω*+kyv¯, where ω*=κky/(1+k2) is the drift wave frequency. However, at larger scales, i.e., scales comparable to ρs corresponding to k1, the dispersion relations differ. In the mHM model, the mean flow leads to a simple Doppler shift in the dispersion relation, as expected, but in the oHM model, the Doppler shift is reduced by the factor k2/(1+k2).

It is clear from the dispersion relations above that in the absence of mean flow or in the presence of a steady, uniform mean flow in the poloidal direction, the Hasegawa-Mima models do not have instabilities. This does not remain true if we assume the presence of a radially varying zonal mean flow. As we will show in an article currently in progress and more focused on the HM models, a drift instability can grow in these conditions and break up the zonal jet structure. Here too, the oHM and mHM models differ, and we find that the zonal jets are more likely to be broken up by this instability in the oHM model than in the mHM model. This provides a hint for the most important physical difference between the oHM and mHM models from the point of view of the present article as well as applications to magnetic confinement fusion, namely, the well-known result that the proper treatment of the electron adiabatic response in the mHM model leads to much stronger zonal jet structures when a random forcing term is added to the equation in order to mimic an instability leading to turbulent behavior.5,7,11 This crucial observation provides the motivation for the mHW model and our new bHW model, which will be the focus of the remainder of this article.

The Hasegawa-Mima models capture essential features of the drift wave–zonal flow dynamics when a forcing term is included in the equation but does not include any internal drift instability to drive the drift wave–zonal flow feedback loop in the absence of forcing. The Hasegawa-Wakatani models21,31 address this limitation by including electron-ion friction, which relaxes the slaving relation between the electron density and the electrostatic potential, leading to a drift wave instability.30 Because the one-to-one correspondence between density and potential is lost, the HW models are two-field models for the ion vorticity and the ion density fluctuation. As for the HM models, one can make the distinction between the original HW model (oHW)31 and the modified HW model (mHW),21 which differ in their treatment of the parallel current. The mHW model accounts for the fact that zonal modes, with wavenumber ky = 0, do not contribute to the parallel current, while the oHW does not. As for the HM models, we can write the two HW models in a unified form, with the switch parameter j such that j =0 for the oHW model and j =1 for the mHW model

(2)
ζt+J(φ,ζ)=α[(φ̃ñ)+δj0(φ¯n¯)]+μΔζ,
(2a)
nt+J(φ,n)+κφ̃y=α[(φ̃ñ)+δj0(φ¯n¯)],
(2b)
where ζ=Δφ is the ion vorticity and n=n1/n0 is the relative density fluctuation, with N=n0+n1 the total ion density, and where all the quantities have been normalized in the same way as for the HM models. The term μΔζ is an approximate model for collisional ion viscosity perpendicular to the magnetic field, where μ=μ̃/(ρs2ωci) with μ̃=3Tiνii/(10miωci2) the kinematic ion viscosity coefficient, Ti the ion temperature, νii the ion-ion collision frequency, and mi the ion mass. We note that strictly speaking, Eqs. (2a) and (2b) are not the mHW model, in the sense that the mHW model has different dissipation terms than the oHW model. Specifically, in the mHW model, dissipation for the ion vorticity is written as DΔ2ζ and dissipation for the ion density is written as DΔ2n, with D the same unspecified constant in both equations. We do not dwell further on this distinction, since these terms are not given any physical justification in Ref. 21, and seem to only be included in order to guarantee numerical stability.

The parameter αTekz2/n0e2ηωci, where η is the resistivity parallel to the field lines, is often referred to as the adiabaticity parameter.4,21,24 In the collisionless limit, α, the electrons have an adiabatic response along the field lines, whereas for small α electron-ion friction decouples the density and electrostatic potential. In order to make the mathematical connection with the HM models in the asymptotic collisionless limit α and μ0, we can derive an equation for the potential vorticity q=2φn valid for both models, by assuming that the oHW and mHW models both only include the viscous dissipation μΔζ as the dissipation mechanism (these dissipative terms will not matter in the end since we will take the limit μ0). Subtracting Eq. (2b) from Eq. (2a), we find

qt+J(φ,q)κφ̃y=μΔζ,q=2φn.
(3)

We see that in the collisionless limit μ0, Eq. (3) has the same form as the equivalent equation (1) for the HM models, with the HM potential vorticities 2φ(φ̃+δj0φ¯) replaced with the HW potential vorticity q=2φn. In the collisionless limit, α, Eqs. (2a) and (2b) give n = φ in the oHW model, so the oHW potential vorticity becomes q=2φφ, and the oHW model coincides with the oHM model as desired.

However, in the mHW model, α leads to ñ=φ̃, leaving n¯ unspecified. The potential vorticity converges to q2φφ̃n¯ and the mHW equations (2a) and (2b) do in general not coincide with the mHM model in the zero resistivity limit. Physically, the mHW model does not guarantee the absence of a net radial transport of electrons when the electrons are adiabatic. Throughout this article, we will return to this important conclusion, which is confirmed by our numerical simulations, as it is a key difference between the mHW model and the new bHW model we introduce in Sec. III. Even if the mHW model does not converge to the mHM model in the collisionless limit, it has been shown that the improved treatment of the parallel current as compared to the oHW model leads to stronger zonal mean structures than in the oHW model,21 just as the mHM model does compared to the oHM model.

 Appendix A contains a detailed discussion of the linear stability of resistive drift waves in the HW models. At this point, it suffices to say that for α0,κ0, and μ = 0, all the modes are unstable and the growth rates are smaller for small-scale modes, k1.

We are now ready to introduce the new bHW model, which addresses a major shortcoming of the mHW model, namely the fact that it does not converge to the mHM model in the collisionless limit. In this new model, we also consider generalized dissipation effects, which allow us to study separately the roles of the different dissipation terms, which are often included for the mere sake of numerical stability, and whose simple forms are rarely justified physically.

Let us start with the dissipation terms. Starting from the mHW model in Eq. (2), we relax the ad hoc assumption that the viscosity coefficients are identical for the ion vorticity and the density fluctuation. Furthermore, we include the model Landau damping term considered in the original Hasegawa-Wakatani article when the drift wave frequency ω* is comparable to kzvTi, where vTi is the ion thermal velocity.31 The equations become

(4)
ζt+J(φ,ζ)=α(φ̃ñ)+μΔζ+Cω*φ,
(4a)
nt+J(φ,n)+κφ̃y=α(φ̃ñ)+DΔn.
(4b)
Here, μΔζ is the same ion collisional viscosity term as in (2a), and DΔn can be viewed as a dissipation term due to the collisional diffusion of the electrons perpendicular to the magnetic field.4 Finally, Cω*φ, with C=Ti/(ωciTe), is the simple model for ion Landau damping introduced in Ref. 31. Observe that Hasegawa and Wakatani suggest a cutoff for this term, which they set to zero when ω* is large. For simplicity, we do not impose such a cutoff, noticing that the absence of cutoff will only be noticeable when kx is small and ky is moderate, since ω* is small for large k, as well as for small ky. Furthermore, even if ω* can be larger than kzvTi when kx is small and ky moderate, the ion viscosity term, μΔζ, is still likely to dominate in these circumstances, since it contains a term proportional to ky4.

We would like to stress the fact that to the best of our knowledge, the forms of our dissipation terms μΔζ,DΔn, and Cφ cannot be derived rigorously from kinetic theory for magnetic confinement plasmas. However, our choice to consider distinct coefficients μ and D and to introduce different forms of dissipation allows us to analyze the relative importance of these effects and determine which effects would need to be modeled more accurately in a higher fidelity model. For example, as we will discuss in more detail in Sec. V in the light of our numerical results, and in our linear stability analysis in  Appendix A, we find that Landau damping as modeled in Eq. (4) has a stabilizing effect on the largest scales and at the same time increases the variability of the small-scale fluctuations.

We now turn to the key modification to the mHW model we make in the bHW model. As we did in Sec. II, if we subtract Eq. (4b) from Eq. (4a), the terms with the adiabaticity parameter cancel, and we obtain an equation for the potential vorticity in the mHW model qmζn=Δφn which reads

qt+J(φ,q)κφ̃y=[(μD)Δ2+Cω*]φ+DΔq.
(5)

As discussed in Sec. II, the issue with the equation above is that in the limit μ,D,C0 and α it violates the balanced electron response on the magnetic surfaces in the mHM sense, due to the inclusion of a non-zero zonal mean state n¯ in the potential vorticity. We propose a simple modification to the mHW model to address this limitation, namely, to replace the original potential vorticity qm with the corrected form qb=ζñ. The resulting flux-balanced equation for the new potential vorticity qb does not depend on the zonal mean density n¯ explicitly. The immediate and desired implication of this modification is that in the collisionless limit, α gives the slaving relation ñφ̃ from Eq. (4b), and the potential vorticity qb converges to the mHM potential vorticity q2φφ̃ unlike the mHW potential vorticity. As a result, in the adiabatic limit and in the absence of the dissipative terms, Eq. (5) when expressed for qb is identical to Eq. (1) of the mHM model.

We now have all the elements to introduce the flux-balanced Hasegawa-Wakatani (bHW) equations for the balanced potential vorticity qb=2φñ and the relative particle density n=n¯+ñ as

(6)
qbt+J(φ,qb)κφ̃y=[(μD)Δ2+Cω*]φ+DΔqb,qb=Δφñ,
(6a)
nt+J(φ,n)+κφ̃y=α(φ̃ñ)+DΔn.
(6b)
The potential vorticity equation (6a) has a form analogous to the vorticity equation in the mHM model, Eq. (1), but unlike the mHM model, it contains a resistive drift instability through the coupling with the density fluctuation ñ. The equation for the relative particle density fluctuation (6b) is the same as in the mHW equation (2b), except for the dissipation term DΔn.

In this section and Sec. III C, we present some elementary properties of our model. We begin with a direct comparison with the mHW model. A good starting point is to write the equations for the balanced potential vorticity qb=Δφñ according to the mHW model and to the bHW model side by side

mHW:qbt+J(φ,qb)+(ũñ¯)x+(n¯xκ)φ̃y=dissip.,
(7a)
bHW:qbt+J(φ,qb)κφ̃y=dissip.,
(7b)

with ũ=φ̃/y the velocity fluctuation along the x direction and

(ũf̃¯)xJ¯(φ̃,f̃)=1Ly(φ̃xf̃yφ̃yf̃x)dy,

with f̃=ñ as in the equation above, or f̃=q̃ as we will have below. The two differences in the mHW model are (1) the zonal density gradient feedback term (xn¯)φ̃/y, which modifies the original background density gradient profile κ=d(lnn0)/dx, and (2) the additional eddy flux feedback, x(ũñ¯), due to the transport of the particle density along the x direction. Both terms are zonal mean density feedbacks which can act on the largest scales (to change the zonal jet structure) and can also modify the structures in the small-scale fluctuations (to change the particle flux). It is important to observe that although the zonal mean density n¯ does not appear in Eq. (7b), qb depends on n¯ in the bHW model through the dependence of ñ on n¯, which we will write explicitly shortly. For the same reason, the vorticity ζ=2φ depends on n¯ in the mHW model.

Let us now turn to the equations for the zonal mean quantities q¯b(x) and n¯(x). If, in order to highlight similarities between the mHW and bHW models, we modify the ad hoc dissipation terms in the mHW model in such a way that the dissipation in the density equation is DΔn in both models, we obtain the following equations for the mean balanced potential vorticity q¯b=2φ¯/x2:

(8)
mHW:q¯bt+(ũq̃b¯)x+(ũñ¯)x=D2q¯bx2,
(8a)
bHW:q¯bt+(ũq̃b¯)x=[(μD)4x4+Cω*]φ¯+D2q¯bx2,
(8b)
and an equation for the zonal mean density n¯ which is identical in both models

n¯t+(ũñ¯)x=D2n¯x2.
(9)

Equations (8) and (9) highlight the feedback of the fluctuations on the zonal mean states, through the nonlinear coupling terms x(ũq̃b¯) and x(ũñ¯). In the mHW model, the density fluctuation feedback x(ũñ¯) is also present in the equation for the mean balanced potential vorticity (8a), just like it is in Eq. (7a). Considering the steady-state version of Eq. (9), we see that a strong particle flux ũñ¯ can lead to a large zonal mean density gradient n¯/x which can then become a significant contribution in the term (xn¯κ)φ̃/y in Eq. (7a) and effectively alter the background density gradient in the mHW model.

A summary of the similarities and differences between the bHW and mHW models is given in  Appendix B, including the equations for the total energy E and enstrophy W which we will discuss in Sec. III C. Once more, in order to draw the appropriate parallels between the two models, we modified the form of the ad hoc dissipation originally proposed in the mHW model21 in such a way as to make the dissipation term agree in the density equation of both models.

For the analysis of the fundamental properties of a new model and the verification of numerical codes, it is always helpful to identify dynamical invariants. We define the total energy E and potential enstrophy W21 for the bHW model as

E=E¯+Ẽ=12Ω(|φ¯x|2+n¯2)dV+12Ω(|φ̃|2+ñ2)dV,W=W¯+W̃=12Ωqb¯2dV+12Ωqb̃2dV,
(10)

where the integrals are over the computational domain Ω, which is a periodic box with sides Lx and Ly, and where we have separated E¯,W¯, the energy and enstrophy of the zonal mean state with qb¯=2φ¯/x2, from Ẽ,W̃, the energy and enstrophy of the fluctuations about the zonal mean state, with qb̃=Δφ̃ñ. Note that the enstrophy W in the bHW model is defined in terms of the balanced potential vorticity qb=2φñ. It is easy to verify that the nonlinear terms J(φ,qb) and J(φ,n) in (6) conserve both the energy and the balanced enstrophy, and we obtain the following dynamical equations for the total energy and potential enstrophy:

dEdt=κΩũñ¯(1+κ1v¯)dVαΩ(ñφ̃)2dVDE,dWdt=κΩũñ¯dVDW,
(11)

where v¯=φ¯x and DE and DW come from the dissipation terms in the bHW equations

DE=Ω(Cω*φ2+μ|Δφ|2+D|n|2)dV,DW=Ω[(Dμ)Δφ̃Δñ+Cω*ñφ̃]dV+Ω[Cω*|φ|2+μ|qb|2+(Dμ)|ñ|2]dV.

Comparing the energy and enstrophy equations (11) with the analogous equations for E and W in the mHW model, as given in  Appendix B, we notice the additional term Ωv¯ũñ¯dV in the energy equation for the bHW model. This contribution to the total energy from the fluctuations originates from the eddy diffusivity term (ũñ¯)x in the equation for the vorticity in the bHW model, which is absent in the mHW model. This additional term, in which we recognize the advection by the mean velocity v¯, represents the zonal flow transport of the particle flux ũñ¯.

As in the mHW model, the resistive coupling of ñ and φ̃ through the adiabaticity α(φ̃ñ) gives a negative-definite term in the energy equation and acts as an energy sink. That term does not enter in the evolution equation for the enstrophy W. The dissipation for the total energy DE is always positive definite. This is not true for the dissipation DW of the potential enstrophy, whose terms may not all always be positive. Specifically, the first integral on the right hand side contains interactions between the fluctuating potential and density Δφ̃Δñ,ñφ̃ due to ion Landau damping and the fact that the damping coefficients μ and D are not equal; the last term in the second integral, (Dμ)|ñ|2, is always negative if D<μ and can thus act as a source of potential enstrophy instead of a sink. In Sec. V B, we illustrate these different dissipation effects on the flow field and energy spectra through direct numerical simulations.

To close this section, we propose a more mathematically rigorous proof of the convergence of the bHW model to the mHM model in the adiabatic limit α. We have already observed that formally, the limit α implies that φ̃=ñ, meaning that the bHW potential vorticity is equal to the mHM potential vorticity and that the models are unified in this limit. Here, we rely on the energy and enstrophy equations (11) to find a more precise description of the relation between the bHW model and the mHM model in the adiabatic limit. For simplicity, we assume that μ=D and C =0, so that the dissipation terms in the energy/enstrophy equations take the simple form DW=D|qb|2 and DE=D(|Δφ|2+|n|2). We then consider states in statistical equilibrium, i.e., such that the time derivatives for the statistical expectations E and W vanish. Here, · represents an ensemble average, which could in principle be estimated with a Monte-Carlo approach, using a large number of numerical simulations. When statistical equilibrium has been reached, the enstrophy equation in (11) gives us the following estimate for the total particle flux:

ΓeqΩũñ¯dVeq=κ1DΩ|qb|2dVeq>0,

where ·eq can be viewed as the time average along the stationary trajectory of the functional f by making the usual ergodic hypothesis. This relation shows that there always exists a statistically positive particle flux towards the boundary of the domain and that the total particle flux Γeq depends on the ratio D/κ.

Turning to the equilibrium statistical equation for the total energy, we obtain the following estimate from the balance between the flux and dissipation terms:

Ω(ñφ̃)2dVeq=α1(κΩũñ¯(1+κ1v¯)dVDE)eqα1(VDWDE)eq.

To derive the inequality above, we used the identity κΓeq=DWeq and assumed that the zonal mean flow can be bounded as V=maxΩ|1+κ1v¯|. The statistical expectation on the right hand side of the inequality is positive and finite in the equilibrium state. Therefore, the expectation (ñφ̃)2eq vanishes as α since the right-hand side of the inequality goes to zero. This demonstrates that in the adiabatic limit, the fluctuations ñ and φ̃ approach the Boltzmann relation ñ=φ̃ in the mean square sense under expectation. Note however that there may still exist a non-zero zonal structure in the density n¯(x) in this limit, in addition to the fluctuation ñ=φ̃.

For the remainder of this article, we turn to direct numerical simulations to highlight characteristic features of the bHW model and compare them with the mHW model. We solve the equations on a doubly periodic square domain of size L along each side, so that the smallest wavenumber is 2π/L, which is also the spacing Δk between any two wavenumbers. We write the quantities (φ,n,ζ) we solve for as the following Fourier series:

φ=kφ̂k(t)eik̃·x,n=kn̂k(t)eik̃·x,ζ=kk̃2φ̂k(t)eik̃·x,

where the spatial variables x=(x,y)[L/2,L/2]×[L/2,L/2] and the corresponding wavenumbers are

k̃=2πLk,k=(kx,ky)[N2+1,N2]×[N2+1,N2].

where N is the number of modes we keep in our simulations. For all the results shown in this article, we used L =40 and N =256. We know from our linear stability analysis (see  Appendix A) that the larger the domain size L, the more unstable modes we will have inside the circle of strong instability. We solve the equations using a pseudo-spectral approach, in which the nonlinear terms are calculated in real space instead of Fourier space. To stabilize the truncated numerical system, hyperviscosities νΔsqb and νΔsn are added in the vorticity and density equations, respectively, with ν=7×1021 and order s =8. We rely on a fourth-order explicit-implicit Runge-Kutta scheme for time stepping, where the stiff hyperviscosity operator νΔ8 is the only term treated implicitly. In agreement with the difference in the formulations of the mHW and bHW models, we numerically solve for the unknowns (ζ,n) when we consider the mHW model and solve for the balanced potential vorticity qb=2φñ and n when we consider the bHW model.

We keep the background density gradient fixed at κ=0.5 and will vary the value of α in the range [0.01,10]. In this section, the dissipation coefficients μ and D are set to μ=5×104,D=5×104, and the ion Landau damping term is turned off by setting C =0. In contrast, in Sec. V we will consider different values for the parameters μ,D, and C to analyze the role of dissipation effects more closely. The parameters we used for all our numerical simulations in this article are summarized in Table II.

TABLE II.

Model parameters and their chosen values for our numerical simulations.

Domain size L 40 
Spatial resolution N 256 
Time step Δt 1×102 
Mean density gradient κ 0.5 
Adiabaticity parameter α 0.01–10 
Collisional ion viscosity parameter μ 5×104,20×104 
Cross-field diffusion D 5,20×104 
Ion Landau damping C 0, 0.01 
Hyperviscosity ν 7×1021 
Order of hyperviscosity s 
Domain size L 40 
Spatial resolution N 256 
Time step Δt 1×102 
Mean density gradient κ 0.5 
Adiabaticity parameter α 0.01–10 
Collisional ion viscosity parameter μ 5×104,20×104 
Cross-field diffusion D 5,20×104 
Ion Landau damping C 0, 0.01 
Hyperviscosity ν 7×1021 
Order of hyperviscosity s 

1. Transition from the turbulence dominated regime to the zonal flow dominated regime

Previous work has shown that the randomness of the turbulent flow depends sensitively on the parameters κ and α.21 This can be readily understood from the linear stability analysis we present in  Appendix A, since κ/α determines the size of the region of instability in k-space. Hence, if κ is held fixed, a smaller value of α, corresponding to higher plasma resistivity, leads to a larger region of instability and a more energetic and turbulent vorticity field. When the value of α is increased, the vorticity field is regularized into dominant anisotropic zonal jets. This is the phenomenon we now study, comparing the bHW and mHW models.

Because of the turbulent nature of the flows, it is more appropriate to adopt a statistical viewpoint, and instead of looking at physical quantities at a given instant in time, we take time-averages of the solution once the stationary state is reached. We will focus on the total kinetic energy and the relative enstrophy for both models, given by

Ω|φ|2dV=kk̃2|φ̂k|2,Ωζ2dV=kk̃4|φ̂k|2.

The relative enstrophy gives more weight to the small scale modes, while the kinetic energy puts more emphasis on the large scale structures. We decompose φ and ζ into their time averages φeq,ζeq and their time fluctuating parts φ,ζ according to φ=φeq+φ,ζ=ζeq+ζ, and compute the energy and enstrophy in the time-averaged mean

Ω|φeq|2dV,Ωζeq2dV,

and the time-averaged kinetic energy and enstrophy of the variances

Ω|φ|2eqdV=Ω|φφeq|2dV,Ωζ2eqdV=Ω(ζζeq)2eqdV.

The results are shown in Fig. 2 for simulations of the bHW equations and the mHW equations, with α varying in the interval [102,10]. Since one of the main reasons to study the reduced HM and HW fluid models is to better understand the drift wave–zonal flow interplay, we separately plot the kinetic energy and enstrophy contained in zonal modes, i.e., with ky = 0, and the total kinetic energy and enstrophy obtained by summing over all the resolved modes. We first focus on the plots for the variances, in the left-hand column, to observe that in both models, the kinetic energy and enstrophy of the variances decrease as α increases. This is a clear signature of the transition from a strongly turbulent regime to a zonal flow dominated regime, as we expected. For α small, the flow is strongly fluctuating with large and small scale modes both containing considerable amount of energy. As α increases, stronger zonal jets develop and the variances in kinetic energy and enstrophy decrease. In the regime α1, the system tends to the HM equations which do not have an internal instability, and the variances reach minimum values.

FIG. 2.

Comparison of the statistical energy of the variance Ω|ffeq|2dV and of the mean Ω|feq|2dV for both the bHW and mHW models and varying values of α, with f=φ (corresponding to the kinetic energy) and f=ζ (corresponding to the enstrophy). The dashed lines show the energy contained in the zonal modes ky = 0. The other parameters are kept fixed, with κ=0.5, μ=D=5×104. (a) Kinetic energy Ω|φ|2dV of the variance and of the mean. (b) Relative enstrophy Ωζ2dV of the variance and of the mean.

FIG. 2.

Comparison of the statistical energy of the variance Ω|ffeq|2dV and of the mean Ω|feq|2dV for both the bHW and mHW models and varying values of α, with f=φ (corresponding to the kinetic energy) and f=ζ (corresponding to the enstrophy). The dashed lines show the energy contained in the zonal modes ky = 0. The other parameters are kept fixed, with κ=0.5, μ=D=5×104. (a) Kinetic energy Ω|φ|2dV of the variance and of the mean. (b) Relative enstrophy Ωζ2dV of the variance and of the mean.

Close modal

Still looking at this first column, we highlight differences between the bHW and mHW models. The bHW model always contains a large variance for the kinetic energy in comparison to the mHW model, which has weaker variability. Furthermore, much of the variance in the bHW model is in zonal modes, unlike the mHW model. Our improved treatment of the electron dynamics parallel to the field lines enhances the zonal flow variability as the system transitions to the turbulence dominated regime when α decreases. The plot for the variance of the relative enstrophy, which is a physical quantity which gives more weight to smaller scale structures, helps to explain the different dynamics in the two models. Indeed, we notice that the mHW model has a much larger variance in enstrophy as α<0.1 than the bHW model. This is due to the generation of many small and strong vortices, as Fig. 5 in Sec. V illustrates explicitly with a contour plot of the vorticity field. Even if so, the zonal mode variance in the relative enstrophy remains larger in the bHW model. This is a central point of this article: even for very small values of α, the variability of the turbulent fluctuations is in zonal modes in the bHW model, while the mHW equations describe fully homogeneous turbulence in the small α limit.

We now consider the kinetic energy and enstrophy in the mean, in the second column of Fig. 2. The kinetic energy and enstrophy in zonal modes almost overlap with the total kinetic energy and enstrophy, proving that the mean state is always dominantly in the zonal direction. We can distinguish three different regimes depending on the value of α. Starting at around α0.1, we notice a significant jump in the value of the kinetic energy and the enstrophy in the mean as compared to situations with α<0.1, indicating the presence of stronger zonal jets. A close comparison between the plots for the kinetic energy and enstrophy shows a large kinetic energy and a moderate enstrophy for α(0.1,0.5), which is a signature of the presence of one or very few zonal jets. For α>0.5, the kinetic energy is a bit lower and fairly constant as a function of α, whereas the enstrophy increases. This indicates the presence of multiple zonal jets. These conclusions are confirmed by the snapshots of the vorticity field for each of the α regimes discussed here, which are shown in Fig. 5 in Sec. V.

2. Statistical energy equations and higher-order moment feedbacks

Besides the kinetic energy and enstrophy in the statistical equilibrium state discussed in Sec. IV A 1, it is also important to investigate the dynamical evolution of the statistical quantities, in particular to better understand the nonlinear energy exchange processes in the model and the model sensitivity to various perturbations. While in many situations we are mostly interested in the statistics of the two lowest order moments, that is, the statistical mean and variance, the higher-order feedbacks to these low order moments play a key role in the dynamics and should be accurately treated, as has been highlighted and illustrated in Refs. 19, 20, and 25 for general systems with quadratic nonlinearity. In this section, we highlight these mechanisms for the bHW and mHW models.

We consider here the statistical average k̃2|φ̂k|2 of the kinetic energy in each spectral mode φ̂k of the electrostatic potential. To derive the dynamical equations for k̃2|φ̂k|2, we first project the equations for the vorticity in the mHW and bHW models onto each individual spectral mode, to obtain equations for each φ̂k. We then multiply by φ̂k* on both sides of the equations and take the statistical ensemble average (see Ref. 20 for details) to obtain the desired equations for k̃2|φ̂k|2. For the bHW model, the equations are

ddt12k̃2|φ̂k|2+QF,kb=α(|φ̂k|2φ̂k*n̂k)Dk̃4|φ̂k|2+c.c.,ky0,ddt12k̃2|φ̂k|2+QF,kb=Dk̃4|φ̂k|2+c.c.,ky=0.
(12)

Here, c.c. stands for the complex conjugate completion. Importantly, the nonlinear interaction term in the vorticity equation produces the higher-order feedback QF,kb=φ̂k*J(φ,Δφ+n¯)k+c.c. for the non-zonal modes ky0; and QF,kb=φ̂k*J(φ,Δφñ)k+c.c. for the zonal modes ky = 0. For the mHW model, the equations are

ddt12k̃2|φ̂k|2+QF,km=α(|φ̂k|2φ̂k*n̂k)Dk̃4|φ̂k|2+c.c.,ky0,ddt12k̃2|φ̂k|2+QF,km=Dk̃4|φ̂k|2+c.c.,ky=0,
(13)

with the higher-order feedback QF,km=φ̂k*J(φ,Δφ)k+c.c. The removal of the zonal mean density n¯ in the term involving the adiabaticity parameter α in both models, as well as in the definition of the potential vorticity in the bHW model, is responsible for the different form the equations take for the zonal modes. We finally note the presence of the cross-covariance term φ̂kn̂k* in both models, representing the interactions between the potential and density.

As emphasized in Refs. 19 and 20, the third-order moments on the left hand side of Eqs. (12) and (13) play the critical role of mediating the growing unstable modes and driving the system to the final equilibrium. As we know from the linear stability analysis, the linear operators on the right hand sides of the statistical equations contain positive eigenvalues, corresponding to the linearly unstable directions. If the third order moments on the left hand sides are not included, the internal instability will lead to fast growth in kinetic energy among the unstable modes and fast decay in the other over-damped modes. When the third-order moments are kept, they transfer the growing energy from the unstable subspace to the stable subspace, where the energy eventually gets dissipated.

We thus see the importance of a rigorous analysis of the contributions from the third-order moments. Unfortunately, it is usually expensive to compute these moments from direct numerical simulations since it requires the inclusion of all the triad modes across the entire spectrum.20 On the other hand, the third-order moment feedbacks in statistical equilibrium are much easier to compute. Indeed, in statistical equilibrium, the time derivatives on the left hand sides of Eqs. (12) and (13) vanish, so the third-order moments can be expressed in terms of the first and second order moments appearing on the right-hand sides. Specifically, for the bHW model, we can write

QF,k,eqb=α(|φ̂k|2eqφ̂k*n̂keq)+Dk̃4|φ̂k|2eq+c.c.,ky0,=Dk̃4|φ̂k|2eq+c.c.,ky=0.

The equilibrium moments, |φ̂k|2eq,φ̂k*n̂keq, can be readily calculated from the equilibrium mean and variance by averaging along the stationary time trajectory. The equilibrium third-order moments for the mHW model can be computed in an analogous way.

To numerically investigate the differences and similarities between the mHW and bHW models from the point of view of the statistical moments, we choose a regime with small α, namely, α=0.01, which corresponds to a physical situation with strong mean flow–drift wave interactions. Figure 3 presents the equilibrium statistics of the bHW and mHW models for both the second-order variances and the third-order moment feedbacks in the spectral domain. We first look at the first column, which shows contours of the logarithm of the second-order variance of each energy mode k2|φ̂k|2. The stronger variability for the zonal modes (ky = 0) in the bHW model is clearly visible. This is consistent with the observation we made in Sec. IV A 1 that even in the small α regime, the bHW model generates energetic zonal flows which are responsible for most of the flow variability. In contrast, we note that the contour plot for the mHW model is radially symmetric in the Fourier domain, indicating the absence of anisotropy in the flow statistics.

FIG. 3.

Equilibrium second-order moments and third-order moment feedbacks in the Fourier domain for the bHW and mHW models with α=0.01,κ=0.5. The first column shows the logarithm of the variance k̃2|φ̂k|2 of each spectral mode. The second column plots the logarithm of the third-order moment feedbacks. Negative values, which correspond to effective damping, are plotted with solid lines and positive values, which correspond to an injection of energy, are plotted with dashed lines. In the third column, the contours of the third-order moment feedbacks are plotted without taking the logarithm, thus emphasizing the modes with the strongest third-order moments. Blue colors are for the largest negative values and red colors for the largest positive values. (a) bHW model and (b) mHW model.

FIG. 3.

Equilibrium second-order moments and third-order moment feedbacks in the Fourier domain for the bHW and mHW models with α=0.01,κ=0.5. The first column shows the logarithm of the variance k̃2|φ̂k|2 of each spectral mode. The second column plots the logarithm of the third-order moment feedbacks. Negative values, which correspond to effective damping, are plotted with solid lines and positive values, which correspond to an injection of energy, are plotted with dashed lines. In the third column, the contours of the third-order moment feedbacks are plotted without taking the logarithm, thus emphasizing the modes with the strongest third-order moments. Blue colors are for the largest negative values and red colors for the largest positive values. (a) bHW model and (b) mHW model.

Close modal

We next turn to the third-order moment feedbacks shown in the second and third columns of Fig. 3. In the middle column, the logarithm of the magnitude of the third-order moments is plotted in order to amplify the small values, and we use dashed lines to represent positive values and solid lines for negative values. In the last column, we plot the contours of the third-order moments themselves (i.e., not their logarithm). The red colors correspond to positive values and the blue colors to negative values. We first observe that both the bHW and the mHW models show strong non-zero contributions from the third-order moments, particularly at the largest scales in the system. Comparing the results of the linear stability analysis we present in Fig. 8 of  Appendix A with these figures, we see that the linearly unstable modes are subject to negative third-order moment feedbacks which act as an effective damping mechanism against the growth in energy. At the same time, the linearly stable modes are subject to positive third-order moment feedbacks. This is characteristic of a nonlinear transfer of energy from the unstable modes, whose energy grows in time, to the stable modes, which then dissipate the energy. Note that the third-order moment feedback conserves energy as a whole, so it cannot act as a source or sink of energy, and instead redistributes energy at different scales. Finally, one observes much stronger third-order moment feedbacks for the zonal modes in the bHW model than in the mHW model, which we would expect given our previous results: once again, even in the regime α1, zonal modes play a central role in the dynamics in the bHW model, while they are essentially absent in the mHW model. These results also highlight the need for accurate approximations of the third-order moment feedback terms when applying model reduction strategies to investigate the statistical properties of a system and for uncertainty quantification. In the present case, the differences between the bHW and mHW models, localized in the vicinity of ky = 0, lead to major differences in the observed dynamics.

Both the bHW and the mHW models converge to an HM-like model without instability in the limit of large α, as we confirmed numerically in Sec. IV A. We showed in Sec. III C that the bHW model converges to the mHM equation (1) (with an additional dissipation term) in the limit α, with the balanced potential vorticity converging to its limiting form qb=2φñ2φφ̃. We also explained in Sec. III A that the mHW equations do not guarantee exact convergence to the mHM model due to the incomplete treatment of the zonal mean density n¯. In this section, we verify our reasoning with numerical simulations. Specifically, we consider a case in which α is relatively large, α = 5, which is a regime where the linear growth rates are small, so the energy in the small-scale modes is mostly dissipated, and a dominant zonal mode with multiple jets sets in, as is clearly visible in the first column of Fig. 4, which shows contour plots of the ion vorticity ζ at the final time of the simulation, in the top left corner for the bHW model, and in the bottom left corner for the mHW model. To obtain a more quantitative idea of the selective energy decay in each spectral mode, we plot the kinetic energy k̃2|φ̂k|2 of each mode at different simulation times in the second column of Fig. 4, with the bHW results at the top and the mHW results at the bottom as before. The energy decay in the small scale modes can be observed for both models at early times, and for both models, the decay saturates at around t =2500. At the final time, only a dominant zonal mode with three zonal jets remains in the flow field, as can be seen in the contour plots in the first column of Fig. 4. However, the decay rate in the mHW is much slower, with a considerable amount of kinetic energy maintained in the small scales, while the small scale kinetic energy in the bHW model decays quickly to negligible amounts after a short transient state. This fundamental difference between the mHW and bHW models is clearly apparent in the contour plots of the vorticity field in the first column. In the bHW simulation in the upper left corner, the vorticity field is largely dominated by the large scale zonal jet structure, with very limited smaller scale features. On the other hand, the mHW simulation shows many persistent small scale vortices in the vorticity field.

FIG. 4.

(Left column) Contour plot of the ion vorticity at the final time of the simulation for the bHW model (top) and the mHW model (bottom). (Middle column) Kinetic energy spectra at different simulation times for the bHW model (top) and the mHW model (bottom). (Third column) Profile of the zonally averaged velocity in the y-direction for the mHM model in solid blue and the bHW model in dashed black (top) and kinetic energy spectra at different simulation times for the mHM model (bottom). The final state of the bHW simulation was chosen as the initial state for the mHM simulation, with random noise added at small scales. For all these results, we used α=0.5 and the values used throughout Sec. IV for the other parameters.

FIG. 4.

(Left column) Contour plot of the ion vorticity at the final time of the simulation for the bHW model (top) and the mHW model (bottom). (Middle column) Kinetic energy spectra at different simulation times for the bHW model (top) and the mHW model (bottom). (Third column) Profile of the zonally averaged velocity in the y-direction for the mHM model in solid blue and the bHW model in dashed black (top) and kinetic energy spectra at different simulation times for the mHM model (bottom). The final state of the bHW simulation was chosen as the initial state for the mHM simulation, with random noise added at small scales. For all these results, we used α=0.5 and the values used throughout Sec. IV for the other parameters.

Close modal

To further verify the convergence to the mHM model, we numerically solve the mHM equations with dissipation and random noise perturbations added at small scales, starting from an initial state given by the final state of the bHW simulation with parameters κ=0.5,α=5. The results are shown in the right column of Fig. 4. Since the mHM model does not contain any instability, the fluctuations get dissipated while the zonal jet structure is maintained in time. The top right figure shows the final and converged zonally averaged velocity v¯=xφ¯ for both the bHW and mHM model simulations. The zonal mean flow clearly converges to the same saturated limit. We have verified that this is true whatever the initial state is chosen for the mHM simulation. Although not shown in the figure, the contour plot for the mHM ion vorticity has strong similarities with the bHW figure, with the same dominant zonal jets, and very few small scale structures, unlike the mHW results.

This concludes our numerical proof of the convergence of the bHW model to the mHM model in the adiabatic limit and of the incomplete convergence of the mHW model in that same limit.

In this section, we illustrate the key results discussed in Sec. IV with plots of the ion vorticity and mean flow obtained with direct numerical simulations for different values of the adiabaticity parameter α. We also study the separate roles of the dissipation terms μΔζ,DΔn, and Cω*φ by considering different values for the parameters μ, D, and C than what we had in Sec. IV.

In Sec. IV, we highlighted the critical role of the adiabaticity parameter α in determining the turbulent flow regime in the HW models, with a transition to strong zonal jets as α increases. This transition is clearly visible in Fig. 5, which shows contours of the ion vorticity ζ for the bHW model in the first row and for the mHW model in the second row. These figures are obtained for the parameters μ=D=5×104 and C =0, and for the final time of the simulation, which is much later than the time when statistical equilibrium has been reached. The three contour plots in each row correspond to three different values of the adiabaticity parameter α. From left to right, the values are α=0.01,α=0.1, and α=0.5. We see that for α=0.01, many small-scale vortices are visible on top of the background mean vorticity. As the parameter value increases to α=0.1, a stronger single jet pattern begins to form, while still coexisting with many smaller scale fluctuations. Finally, for α=0.5, the vorticity field becomes fully dominated by several strong zonal jets. These figures offer an explicit illustration of the statistical results given in Fig. 2.

FIG. 5.

Snapshots of the ion vorticity ζ=Δφ at the final time of the simulation for the bHW model (top row) and the mHW model (bottom row). For each figure κ=0.5, and each column corresponds to a different value of α: α=0.01,0.1,0.5 from left to right. Notice the clear zonally elongated structures obtained in the bHW model for α=0.01 in comparison with the homogeneous field in the mHW model. (a) bHW model and (b) mHW model.

FIG. 5.

Snapshots of the ion vorticity ζ=Δφ at the final time of the simulation for the bHW model (top row) and the mHW model (bottom row). For each figure κ=0.5, and each column corresponds to a different value of α: α=0.01,0.1,0.5 from left to right. Notice the clear zonally elongated structures obtained in the bHW model for α=0.01 in comparison with the homogeneous field in the mHW model. (a) bHW model and (b) mHW model.

Close modal

Figure 5 also confirms another key observation regarding differences between the bHW and mHW models. The bHW model maintains jet structures for a wide range of values of the adiabaticity parameter α, throughout the transition from a regime with dominant zonal jets (large α) to the strong drift wave turbulence regime (small α). For small values of α, the jets are more turbulent and shift in time, but the anisotropic zonal dynamics persists. This is in stark contrast to the mHW model, which loses the jets as α0, as a regime with fully homogeneous turbulence and strong vortices sets in. Because of the persistence of the zonal jets, the particle flux is always smaller in the bHW model than in the mHW model for small values of α. Physically, the unbalanced density flux in the mHW model is responsible for the highly turbulent vorticity and strong particle transport in the limit of small α.

In the large α regime, we highlight another critical difference between the mHW and the bHW models, namely, the fact that in the bHW model, both the zonal mean flow and the fluctuations have a much larger variability than in the mHW model. To demonstrate this, in Fig. 6 we plot the time series of the zonal mean velocity field v¯=φ¯/x for α=0.5, which corresponds to a strong zonal jet regime. The jets generated in the bHW model have large amplitude variations in time, whereas the zonal velocity in the mHW model is mostly steady in time with an almost constant jet amplitude.

FIG. 6.

Time-series of the zonally averaged mean flow v¯=φ¯x for the bHW model (top) and the mHW model (bottom) in the zonal jet dominated regime α=0.5,κ=0.5.

FIG. 6.

Time-series of the zonally averaged mean flow v¯=φ¯x for the bHW model (top) and the mHW model (bottom) in the zonal jet dominated regime α=0.5,κ=0.5.

Close modal

Thus far, we have always turned off ion Landau damping by setting C =0 and set the dissipation coefficients D and μ to be equal, with μ=D=5×104. In this last section, we take a closer look at the right-hand side of Eq. (6a) by dissociating μ and D and considering finite values for C. Specifically, we study the consequences of increasing the values of μ to μ=2×103 while keeping the other parameters to their original values and of increasing D to D=2×103 while keeping the other parameters to their original values. In addition, we test a small value for C, namely, C =0.01. Since Landau damping mostly acts on the largest scales, we focus in this section on the regime corresponding to α=0.5,κ=0.5, where there exist strong but fluctuating large-scale zonal mean modes.

Figure 7 provides a summary of our main results. The top left panel shows time series of the total energy in the system for the different values of μ, D, and C considered in this section. The blue, red, and yellow curves correspond to cases in which there is no Landau damping, and we see that increasing D by a factor of 4 to D=2×103 only decreases the total energy by a small amount. As expected from our linear stability analysis in  Appendix A, increasing μ by a factor of 4 to μ=2×103 has a slightly stronger effect, decreasing the total energy to a noticeably lower value, but it does not significantly change the nature of the dynamics. In contrast, ion Landau damping greatly changes the energy time series, even for small values of C, as shown with the purple curve. In this regime with strong zonal jets, Landau damping increases the total energy and leads to much larger fluctuations. This is because our naive model for Landau damping damps the large scales but increases the energy at the intermediate and small scales, as we now explain with a spectral viewpoint.

FIG. 7.

(Top row) Time-series of the total energy and time-averaged energy spectra for the bHW model, for different values of the coefficients μ, D, and C. (Bottom row) Time-series of the zonally averaged mean flow v¯=φ¯x for the bHW model when ion Landau damping is included, with μ=D=5×104 and C =0.01. Notice how the jets are more turbulent when C0, with stronger small-scale fluctuations. For all these results, α and κ are fixed, with α=0.5,κ=0.5.

FIG. 7.

(Top row) Time-series of the total energy and time-averaged energy spectra for the bHW model, for different values of the coefficients μ, D, and C. (Bottom row) Time-series of the zonally averaged mean flow v¯=φ¯x for the bHW model when ion Landau damping is included, with μ=D=5×104 and C =0.01. Notice how the jets are more turbulent when C0, with stronger small-scale fluctuations. For all these results, α and κ are fixed, with α=0.5,κ=0.5.

Close modal

To better understand the results found for the time series of the total energy, we show the energy spectra of the variance Etot2eq and the statistical mean Etoteq in the top center and top right panels of Fig. 7, respectively. These figures confirm that increasing D has a small impact on small scale modes and does not significantly modify the largest zonal mean modes of the flow either. The impact of a larger value of μ is more noticeable but approximately uniform across all scales. On the other hand, setting the Landau damping coefficient C to a finite value has a strong effect in changing the zonal flow profile v¯. Landau damping effectively removes the energy of the zonal modes at the largest scales and drives the original zonal flow with 3 jets to a configuration with 5 dominant jets in the flow. At the same time, it increases the variance at all scales and particularly strongly at intermediate scales. We illustrate this effect by showing the time-series of the zonal mean flow v¯=xφ¯ with Landau damping in the second row of Fig. 7. Comparing this figure with the analogous Fig. 6 for the same regime in the absence of Landau damping, we see that the jets are more turbulent when Landau damping is included, meandering in time with stronger fluctuations. Moreover, instead of a persistent 3-jet structure as in the C =0 case, the flow shows a dominant 5-jet structure, with sudden merging events leading to a 4-jet flow, before the 5 jets reemerge after a short period of time.

In this article, we presented a new two-field fluid model for the study of the drift wave–zonal flow dynamics in magnetically confined plasmas. Our model is inspired by the rich tradition of Hasegawa-Wakatani (HW) models,21,24,31 which are the simplest known two-field models for magnetized plasmas which contain a drift-wave instability and can generate zonal flows through nonlinear drift wave interactions. The main novelty of our balanced HW (bHW) model as compared to previous HW models is its improved treatment of the electron dynamics parallel to the magnetic field lines, which guarantees a balanced electron density response7,10 in the limit in which electrons are adiabatic. To achieve this, the bHW model does not solve for the ion vorticity, but instead for a well constructed potential vorticity defined by qb=2φñ, which we call the balanced potential vorticity and does crucially not include the zonal mean density n¯.

Our new bHW model inherits the desirable features of former HW models as well as those of the Hasegawa-Mima model with a modified adiabatic electron response, which we have called mHM.7 Like other HW models, the bHW model contains a drift wave instability, and like the mHW model, the bHW model has the desired Galilean invariance along a magnetic flux surface; on the other hand, like the mHM model, the improved electron response in the bHW model leads to enhanced zonal flows and strong fluctuations. Mathematically, the bHW model has the satisfying property of converging to the mHM model in the proper adiabatic limit, which is not the case of previous HW models.

We relied on direct numerical simulations to investigate the main features of the bHW model, with numerical simulations of the recent modified HW model (mHW)21 as a point of comparison. We paid close attention to the transition from a strongly drift wave turbulent regime to a much more organized regime with strong zonal structures, which is associated with a decrease in the plasma resistivity, allowing electrons to become more and more adiabatic. We observed that the bHW model has stronger fluctuations than the mHW model in the turbulent regime and throughout the transition to the regime with strong zonal jets. Remarkably, however, the bHW model maintains mainly zonal dynamics even at large plasma resistivity, while in contrast the drift wave turbulent regime in the mHW model is dominated by small scale dynamics. In other words, in the bHW model, the zonal flows are more robust, in the sense that they are observed for any value of the plasma resistivity, but they have a larger variability than in the mHW model in the regime in which they are seen in both models. We also used direct numerical simulations to confirm the convergence to the mHM model in the zero resistivity limit and to study the role of the dissipation terms introduced in our model. We noticed that the simple model for Landau damping introduced by Wakatani and Hasegawa31 leads to more turbulent dynamics, with zonal structures which are fluctuating more, and more energy and variability at intermediate and small scales. The significant effect of the Landau damping term on the bHW dynamics provides yet again a clear justification for the development and implementation of more accurate fluid approximations for this effect, as was first done in the pioneering work of Hammett, Perkins, and Dorland.10,12

This article is meant to be the first in a series of articles in which we study the interplay between drift waves and zonal flows using the new bHW model presented here. In a second article,26 more focused on the analysis of detailed numerical simulations, we investigate the effect of the shape and size of the computational domain on the properties of the drift wave turbulence and zonal structures. In a subsequent article, we will treat the bHW model and the drift wave–zonal flow dynamics it describes as a test bed to explore the applicability of statistical model reduction strategies20,25 for the efficient description of turbulent magnetized plasmas. The latter is the topic of ongoing research, with progress to be reported in the near future.

A.J.C. would like to thank Jeffrey B. Parker for the suggestions he made during an insightful conversation which led to this work. A.J.M. was partially supported by the Office of Naval Research through MURI N00014-16-1-2161 and DARPA through W911NF-15-1-0636. D.Q. was supported as a postdoctoral fellow on the second grant. A.J.C. was partially supported by the U.S. Department of Energy, Office of Science, Fusion Energy Sciences under Award Nos. DE-FG02-86ER53223 and DE-SC0012398.

In this appendix, we present a detailed linear stability analysis for the bHW model given by Eq. (6). For this analysis, we assume a zero background mean flow, so the linearized equations for the perturbations (φ1,q1b,n1) to the background state are given by

tq1bκyφ̃1=[(μD)Δ2+Cω*]φ1+DΔq1b,tn1+κyφ̃1=α(φ̃1ñ1)+DΔn1.

The dispersion relation for the bHW model is found by considering a single mode solution for the perturbations, i.e.,

φ1=φ̂1exp(i(k·xωt)),n1=n̂1exp(i(k·xωt)),ζ1=k2φ̂1exp(i(k·xωt)),

where we assume ky0. Substituting the above single mode solution in the linearized equations, we find the following equations for φ̂1 and n̂1:

iωk2φ̂1=α(φ̂1n̂1)+μk4φ̂1+Cκky1+k2φ̂1,iωn̂1=α(φ̂1n̂1)iκkyφ̂1Dk2n̂1.
(A1)

We first derive and analyze the dispersion relation for the situation without dissipation D=μ=C=0. We will then discuss the modifications due to the various dissipation effects.

1. Linear stability analysis without dissipation

From the coupled equations (A1) in the absence of dissipation effects, μ=D=C=0, it is straightforward to derive the dispersion relation

ω2+ib(ωω*)=0,
(A2)

where the drift wave frequency ω* and the wavenumber dependent coefficient b are given by

ω*(k)=κky/(1+k2),b(k)=α(1+k2).

We immediately see that if α = 0, b =0 and ω = 0, so all the modes are marginally stable. Furthermore, when α0 but κ = 0, ω = 0 and ω=ib are solutions, so there are only damped and marginally stable modes. For the general case with α0,κ0, we write the solution to equation (A2) as ω=ωr+iωi, with

ωr=±b2(1+16γ2)1/4cosθ2,ωi=b2±b2(1+16γ2)1/4sinθ2,θ=Arg(1+i4γ),γ=ω*b.
(A3)

Here, Arg stands for the principal value argument. Thus, for the case α0,κ0, linear stability is determined by comparing the term (1+16γ2)sin4θ2 with 1. θ=Arg(1+i4γ) implies that cos2θ=(1+16γ2)1. Hence

(1+16γ2)sin4θ2=14(1+(1+16γ2)1/2)2(1+16γ2)=14(1+16γ2+1)2>1foranyγ0.
(A4)

We conclude that for all modes with ky0, there is a solution with ωi > 0. In other words, in the absence of dissipation, all modes with ky0 are unstable.

Keeping in mind the complete equations including dissipation, it is valuable to construct a boundary curve outside of which the growth rate can be regarded as small compared to dissipation mechanisms. We know that the system is marginally stable if γ=ω*/b=0. The growth rate of linear instability is small if ω*(k)b(k), which we write as

ω*b=ϵ,ϵ1.

The above relation can be rewritten as

(k+k1)2=kyϵκα.

In the asymptotic regime ϵ1, the above equation has one solution for k1 and one solution for k1. We discard the former since for our finite size periodic box, the smallest wavenumber has magnitude 2π/L. We thus have, for ϵ1

k+1k=(kyϵ)1/2(κα)1/2.

The curve separating slowly growing modes from faster growing modes for our drift wave instability is therefore given by the following Cartesian equation in wavenumber space:

kx2+(ky12ϵκα)2=14ϵ2(κα)2.
(A5)

It describes a circle in the spectral domain, centered in (0,κ/2ϵα) and with radius κ/2ϵα. Inside the circle, the modes have large linear growth rates, and outside the circle, the modes are quasi-stable. We note that our analysis shows that the range of wavenumbers which correspond to very unstable modes only depends on the ratio κ/α.

Now, returning to the general situation with arbitrary values for ω*/b, we can combine our expression for the growth rate in Eq. (A3) and the relation (A4) to write

ωi=b2[1+12(1+16γ2+1)1/2].

The magnitude of the term in the square parentheses is determined, for a given wave vector, by the magnitude of the ratio r = κ/α. We can derive insightful explicit formulae for the wavenumber dependence of ωi for the two asymptotic regimes r1 and r1.

For r1,γ1 and we find

ωiα(k2+1)k2γ2=κ2αky2k2(1+k2)3;

For fixed values of κ and α, we can look for the mode with the largest growth rate. Given the symmetry of the linearized equations with respect to the wavenumber kx, we look for a mode with kx = 0 and ky finite. In that case

ωiκ2αky4(1+ky2)3;

The maximum growth rate ωimax=427(κ2/α) is reached for kymax=2, corresponding to γmax=229(κ/α) and bmax=3α/2.

For r1,γ1 and we find

ωiα(k2+1)2k22γ=ακ2ky.

In this regime, the largest growth rate is obtained for the smallest finite ky.

Our linear stability analysis is confirmed by our numerical results. The first two plots on the left in Fig. 8 are contour plots of the growth rate ωi in Fourier space, calculated from Eq. (A3). For these plots, we chose ϵ=2×102 and did not plot contours corresponding to slow growing modes, outside of the circle of meta-stability given by Eq. (A5) and marked by the dashed black lines. The size of the region with strong growing modes increases with the ratio κ/α, as predicted. The axial symmetry with respect to the line kx = 0 is clearly visible, and we verify that the maximum growth rate is reached for kx = 0 and ky2 for κ/α1. The rightmost plot in Fig. 8 shows contours of the maximum growth rate in κ-α space. The dependence on κ and α in the regimes r1 and r1 agrees with our analytical expressions derived for these asymptotic regimes.

FIG. 8.

Linear growth rates from the solution of (A3) with different ratios κ/α. The black dashed line corresponds to the circle of meta-stability given by (A5), with ϵ=2×102. The fast growing modes are located inside the circle. The figure on the right shows the maximum linear growth ωimax rate as a function of the model parameters α and κ.

FIG. 8.

Linear growth rates from the solution of (A3) with different ratios κ/α. The black dashed line corresponds to the circle of meta-stability given by (A5), with ϵ=2×102. The fast growing modes are located inside the circle. The figure on the right shows the maximum linear growth ωimax rate as a function of the model parameters α and κ.

Close modal
2. Modifications due to dissipative effects

We now consider the effects of the dissipation terms, namely, μΔζ+Cω*φ for the vorticity equation and DΔn for the density equation. The linearized equation (A1) can be written as a 2 × 2 system of equations for each wavenumber, given by

[iωαk2(μk2+Cω*(k)k2)αk2α+iκkyiω+α+Dk2][φ̂n̂]=0.
(A6)

Nontrivial solutions to the system exist, provided that the following dispersion relation has solutions:

ω2+i(b+L1)ωibω*(k)αL2D(μk4+Cω*(k))=0,
(A7)

where

L1(k)=Dk2+μk2+Cω*(k)k2,L2(k)=(D+μk2+Cω*(k)k2).

Because of the complexity of the dispersion relation (A7), we will only solve it numerically. To do so, we choose the same numerical values for the parameters as we used in the numerical simulations discussed in the main text. Specifically, if μ0, then μ = 5 × 10–4, if D0, then D =5 × 10–4, and if C0, then C =0.01. Our results comparing linear growth rates with and without dissipation terms are shown in Fig. 9. The first two plots on the left show the dependence of the linear growth rates on ky along the line kx = 0 in k-space, which is the line along which the maximum growth rate is reached. The situation without dissipation is considered in the first plot on the left, for different values of κ and α. The curves confirm our analytical results in Subsection 1 of the  Appendix A, with a maximum reached for ky2 for r1 and moving to the left as r increases. In the second plot, we analyze the effect of including dissipation terms, for α = 0.1 and κ = 0.5. We see that as expected, our simple model for Landau damping stabilizes the large scale modes but has a negligible effect on small scale modes. In contrast, the terms μΔζ and DΔn only stabilize the small scale modes and have little effect on the large scale modes. In order to determine the relative importance of these two terms, we turn to the plot on the right of Fig. 9, which shows the contour of meta-stability, set here at ωi=1×103, for different values of the parameters μ, D, and C, and α = 0.1 and κ = 0.5. We see that the term μΔζ has a much stronger damping effect, since setting μ to zero leads to a much increased region of strong instability. We also have the confirmation that the model Landau damping term stabilizes the large scale modes.

FIG. 9.

Comparison of growth rates with changing parameter values for α,κ, μ, D, and C. The left and middle plots show the linear growth rates as a function of ky, for kx = 0 and different values of α,κ, μ, D, and C. The plot on the right shows the contours of meta-stability with small growth rate ωi=1×103 as the coefficients of the dissipation terms are changed.

FIG. 9.

Comparison of growth rates with changing parameter values for α,κ, μ, D, and C. The left and middle plots show the linear growth rates as a function of ky, for kx = 0 and different values of α,κ, μ, D, and C. The plot on the right shows the contours of meta-stability with small growth rate ωi=1×103 as the coefficients of the dissipation terms are changed.

Close modal

In this appendix, we present a summarized comparison between the modified Hasegawa-Wakatani model and the new flux-balanced Hasegawa-Wakatani model we introduced in this article.

1. Modified Hasegawa-Wakatani (mHW) Model

Governing equations for the vorticity ζ=2φ and the density n

ζt+J(φ,ζ)=α(φ̃ñ)+DΔζ,nt+J(φ,n)+κφ̃y=α(φ̃ñ)+DΔn.

Auxiliary equation for the balanced potential vorticity qb=2φñ, with ũ=φ̃y

qbt+J(φ,qb)+(ũñ¯)x(κn¯x)φ̃y=DΔqb.

Zonal mean and fluctuation equations for the vorticity ζ=ζ¯+ζ̃

ζ¯t+(ũζ̃¯)x=D2ζ¯x2,ζ̃t+J(φ,ζ̃)+[ζ¯xũx(ũζ̃¯)]=α(φ̃ñ)+DΔζ̃.

Conservation laws for the energy E and the enstrophy W

E=12Ω(|φ|2+n2)dV,W=12Ω(ζn)2dVdEdt=κΩñũdVαΩ(φ̃ñ)2dVDE,dWdt=κΩñũdVDW.
2. Flux-balanced Hasegawa-Wakatani (bHW) Model

Governing equations for the balanced potential vorticity qb=2φñ and the density n

qbt+J(φ,qb)κφ̃y=[(μD)Δ2+Cω*]φ+DΔqb,nt+J(φ,n)+κφ̃y=α(φ̃ñ)+DΔn.

Auxiliary equation for the vorticity ζ=2φ, with ũ=φ̃y

ζt+J(φ,ζ)+[n¯xũ(ũñ¯)x]=α(φ̃ñ)+μΔζ+Cω*φ.

Zonal mean and fluctuation equations for the vorticity ζ=ζ¯+ζ̃

ζ¯t+(ũζ̃¯)x(ũñ¯)x=μΔζ¯+Cω*φ¯ζ̃t+J(φ,ζ̃)+[ζxũx(ũζ̃¯)]+n¯xũ=α(φ̃ñ)+μΔζ̃+Cω*φ̃.

Conservation laws for the energy E and the enstrophy W, with v¯=φx¯

E=12Ω(|φ|2+n2)dV,W=12Ω(ζñ)2dVdEdt=κΩñũ(1+κ1v¯)dVαΩ(ñφ̃)2dVDE,dWdt=κΩñũdVDW.
1.
A.
Ashourvan
,
P. H.
Diamond
, and
Ö. D.
Gürcan
, “
Transport matrix for particles and momentum in collisional drift waves turbulence in linear plasma devices
,”
Phys. Plasmas
23
(
2
),
022309
(
2016
).
2.
R.
Balescu
,
Aspects of Anomalous Transport in Plasmas
(
CRC Press
,
Boca Raton
,
2005
).
3.
M. A.
Beer
, “
A gyrofluid models of turbulent transport in tokamaks
,” Ph.D. thesis (
Princeton University
,
1994
).
4.
S. J.
Camargo
,
D.
Biskamp
, and
B. D.
Scott
, “
Resistive drift–wave turbulence
,”
Phys. Plasmas
2
(
1
),
48
62
(
1995
).
5.
C.
Chandre
,
P. J.
Morrison
, and
E.
Tassi
, “
Hamiltonian formulation of the modified Hasegawa-Mima equation
,”
Phys. Lett. A
378
(
13
),
956
959
(
2014
).
6.
C.
Connaughton
,
S.
Nazarenko
, and
B.
Quinn
, “
Feedback of zonal flows on wave turbulence driven by small-scale instability in the Charney-Hasegawa-Mima model
,”
Europhys. Lett.
96
(
2
),
25001
(
2011
).
7.
R. L.
Dewar
and
R. F.
Abdullatif
, “
Zonal flow generation by modulational instability
,” in
Frontiers in Turbulence and Coherent Structures
(
World Scientific
,
2007
), pp.
415
430
.
8.
P. H.
Diamond
,
S. I.
Itoh
,
K.
Itoh
, and
T. S.
Hahm
, “
Zonal flows in plasma—A review
,”
Plasma Phys. Controlled Fusion
47
(
5
),
R35
(
2005
).
9.
W.
Dorland
, “
Gyrofluid models of plasma turbulence
,” Ph.D. thesis (
Princeton University
,
1993
).
10.
W.
Dorland
and
G. W.
Hammett
, “
Gyrofluid turbulence models with kinetic effects
,”
Phys. Fluids B
5
(
3
),
812
835
(
1993
).
11.
G. W.
Hammett
,
M. A.
Beer
,
W.
Dorland
,
S. C.
Cowley
, and
S. A.
Smith
, “
Developments in the gyrofluid approach to tokamak turbulence simulations
,”
Plasma Phys. Controlled Fusion
35
(
8
),
973
(
1993
).
12.
G. W.
Hammett
and
F. W.
Perkins
, “
Fluid moment models for landau damping with application to the ion-temperature-gradient instability
,”
Phys. Rev. Lett.
64
,
3019
3022
(
1990
).
13.
A.
Hasegawa
,
C. G.
Maclennan
, and
Y.
Kodama
, “
Nonlinear behavior and turbulence spectra of drift waves and rossby waves
,”
Phys. Fluids
22
(
11
),
2122
2129
(
1979
).
14.
A.
Hasegawa
and
K.
Mima
, “
Pseudo-three-dimensional turbulence in magnetized nonuniform plasma
,”
Phys. Fluids
21
(
1
),
87
92
(
1978
).
15.
A.
Hasegawa
and
M.
Wakatani
, “
Plasma edge turbulence
,”
Phys. Rev. Lett.
50
(
9
),
682
(
1983
).
16.
A.
Hasegawa
and
M.
Wakatani
, “
Self-organization of electrostatic turbulence in a cylindrical plasma
,”
Phys. Rev. Lett.
59
,
1581
1584
(
1987
).
17.
W.
Horton
, “
Drift waves and transport
,”
Rev. Mod. Phys.
71
,
735
778
(
1999
).
18.
K.
Itoh
,
S. I.
Itoh
,
P. H.
Diamond
,
T. S.
Hahm
,
A.
Fujisawa
,
G. R.
Tynan
,
M.
Yagi
, and
Y.
Nagashima
, “
Physics of zonal flows
,”
Phys. Plasmas
13
(
5
),
055502
(
2006
).
19.
A. J.
Majda
,
Introduction to Turbulent Dynamical Systems in Complex Systems
(
Springer
,
2016
).
20.
A. J.
Majda
and
D.
Qi
, “
Strategies for reduced-order models for predicting the statistical responses and uncertainty quantification in complex turbulent dynamical systems
,”
SIAM Rev.
60
(
3
),
491
549
(
2018
).
21.
R.
Numata
,
R.
Ball
, and
R. L.
Dewar
, “
Bifurcation in electrostatic resistive drift wave turbulence
,”
Phys. Plasmas
14
(
10
),
102312
(
2007
).
22.
J. B.
Parker
, “
Dynamics of zonal flows: Failure of wave-kinetic theory, and new geometrical optics approximations
,”
J. Plasma Phys.
82
(
6
),
595820602
(
2016
).
23.
J. B.
Parker
and
J. A.
Krommes
, “
Generation of zonal flows through symmetry breaking of statistical homogeneity
,”
New J. Phys.
16
(
3
),
035006
(
2014
).
24.
A. V.
Pushkarev
,
W. J. T.
Bos
, and
S. V.
Nazarenko
, “
Zonal flow generation and its feedback on turbulence production in drift wave turbulence
,”
Phys. Plasmas
20
(
4
),
042304
(
2013
).
25.
D.
Qi
and
A. J.
Majda
, “
Low-dimensional reduced-order models for statistical response and uncertainty quantification: Two-layer baroclinic turbulence
,”
J. Atmos. Sci.
73
(
12
),
4609
4639
(
2016
).
26.
D.
Qi
,
A. J.
Majda
, and
A.
Cerfon
, “
A flux-balanced model for collisional plasma edge turbulence: Numerical simulations with different aspect ratios
,” J. Plasma Phys. (to be submitted).
27.
P.
Ricci
,
B. N.
Rogers
,
W.
Dorland
, and
M.
Barnes
, “
Gyrokinetic linear theory of the entropy mode in a Z pinch
,”
Phys. Plasmas
13
(
6
),
062102
(
2006
).
28.
B. N.
Rogers
,
W.
Dorland
, and
M.
Kotschenreuther
, “
Generation and stability of zonal flows in ion-temperature-gradient mode turbulence
,”
Phys. Rev. Lett.
85
,
5336
5339
(
2000
).
29.
G. M.
Staebler
,
J.
Candy
,
N. T.
Howard
, and
C.
Holland
, “
The role of zonal flows in the saturation of multi-scale gyrokinetic turbulence
,”
Phys. Plasmas
23
(
6
),
062518
(
2016
).
30.
T.
Stoltzfus-Dueck
,
B. D.
Scott
, and
J. A.
Krommes
, “
Nonadiabatic electron response in the Hasegawa-Wakatani equations
,”
Phys. Plasmas
20
(
8
),
082314
(
2013
).
31.
M.
Wakatani
and
A.
Hasegawa
, “
A collisional drift wave description of plasma edge turbulence
,”
Phys. Fluids
27
(
3
),
611
618
(
1984
).
32.
R. E.
Waltz
,
G. D.
Kerbel
, and
J.
Milovich
, “
Toroidal gyro-landau fluid model turbulence simulations in a nonlinear ballooning mode representation with radial modes
,”
Phys. Plasmas
1
(
7
),
2229
2244
(
1994
).
33.
T. H.
Watanabe
,
H.
Sugama
, and
S.
Ferrando-Margalet
, “
Gyrokinetic simulation of zonal flows and ion temperature gradient turbulence in helical systems
,”
Nucl. Fusion
47
(
9
),
1383
(
2007
).