Global ITG Eigenmodes: From Ballooning Angle Andradial Shift To .

Physics of Plasmaspressure configurations. We have been developing global simulationsof ITG turbulence using an electrostatic adaptation of the JOREKcode, which is a global finite element fluid code originally designed forsimulations of ELMs in complete X-point geometry. The implementation of this electrostatic fluid description of ITG modes into JOREKtakes advantage of the previous extensive development of the codewith respect to equilibrium, geometry, and numerical methods, whichwere originally applied to MHD simulations. In particular, the realisticequilibrium obtained from the Grad–Shafranov solver, flux-alignedgrid needed for modeling, and all numerical developments of thesparse matrix solver were used here. These features allow us to consider the excitation and interaction of ITG modes across many rationalsurfaces, along with including macroscopic tokamak phenomena. Inthis paper, we first demonstrate successful benchmarking to the lineargrowth rate data in Refs. 31 and 32. Although our equilibria are poloidally symmetric, these simulations do indeed develop the up-downasymmetry in the mode structure, enumerated by the finite ballooning6 0. In varying the magnetic shear, while maintaining theangle, h0 ¼same temperature gradients, we have shown a transition from the typical global ballooning mode structure, where h0 0, to modes whichare radially localized between adjacent rational surfaces, which retaintheir up-down symmetry, h0 ¼ 0, as is typical for slab-like modes. Thetransition to localized modes occurs for low values of magnetic shear, s ¼ dðln qÞ dðln rÞ ⱗ 0:7 (where q is the safety factor and r is theminor radius), at which the elongated ballooning structure breaksdown due to weaker coupling between poloidal harmonics.In this paper, we also study the structure of the quasilinearReynolds stress, which arises due to global geometric effects. Reynoldsstress drives plasma rotation, and can lead to the generation of poloidal shear flows, which have been shown to suppress turbulent transport by stretching and breaking apart turbulent vortices.33 This limitsthe spatial scale over which the vortices can transfer energy, resultingin turbulent self-stabilization. If the flows become strong enough, theycan lead to the formation of transport barriers, which are most likelyresponsible for the transition between low and high confinementregimes (L-mode and H-mode, respectively). Transport barrier formation is commonly studied in the context of zonal flows, which arepoloidally and toroidally symmetric flows, n ¼ m ¼ 0, however thereis growing interest in poloidally asymmetric flows (convective cells),which have been detected experimentally.34 In particular, the m ¼ 1component of the turbulent Reynolds stress is of interest for the generation of Geodesic Acoustic Modes (GAMs).It has been pointed out that the eigenmode asymmetries (i.e., ballooning angle) predicted within the generalized 2D ballooning theorycan also lead to significant residual Reynolds stress.1,20,21 Generalizedballooning theory, which includes global profile effects, also predicts acombined monopole m ¼ 0 and m ¼ 1 dipole structure of theReynolds stress, even in poloidally symmetric equilibria.1,3,6,35 In thispaper, we investigate these mechanisms by studying the poloidal structure of the Reynolds stress depending on the value of the magneticshear and demonstrate significant generation of Reynolds stress due tomagnetic shear, diamagnetic, and global profile effects.Mode saturation in multi-harmonic simulations of multiplen 0 is also demonstrated and analyzed. This occurs for all shear configurations, however, the resultant energy level distribution of eachconstituent n harmonic can be markedly different. Furthermore, forlow shear configurations, we find the possibility of nonlinearPhys. Plasmas 27, 072507 (2020); doi: 10.1063/5.0006765Published under license by AIP on in single-n simulations (without coupling to n ¼ 0), occurring by the coupling of unstable poloidal harmonics to stable ones.This occurs since we use a non-Boussinesq form of the vorticity36which leads to cubic interactions between three variables, allowing asingle toroidal mode n to couple to itself.In Sec. II, we describe the two constituent components of ourITG model, and describe their linear analytic behavior. Following this,in Sec. III we briefly overview JOREK, the code we use for the simulation, outline our calculation of equilibrium, and describe the combinednonlinear model used in simulations. Section IV, describes our benchmarking comparisons to existing simulations, and Sec. V offers a thorough analysis of the two distinct mode types which arise insimulations of a single toroidal mode number. Section VI then investigates the Reynolds stress of each of these mode types. Section VII covers the nonlinear behavior in our simulations, beginning by discussingsaturation occurring via coupling between m modes, then moving onto turbulent states including multiple n 0 modes. Finally, we conclude in Sec. VIII by offering a summary of our results.II. LINEAR ITG THEORYThere are two main types of ITG modes, commonly referred toas the toroidal37–40 and slab41–44 ITGs. Both modes are driven by negative effective compressibility, where an increase in density occurs witha decrease in pressure (@pi @ni 0, where pi and ni are the ion pressure and density, respectively). The toroidal ITG (tITG) mode is aninterchange mode, and is thus only unstable in a particular configuration—namely, when a force is directed against a density gradient. Theprototypical example of an interchange mode is the Rayleigh Taylorinstability, where for a heavy fluid above a light one, the downwardforce of gravity is directed opposite the upward density gradient.Conceptually, this mechanism is the fluid equivalent of a pendulum inan inverted state. For the tITG, the force is provided by the magneticgradient and curvature, and this force opposes the density gradient onthe low field side, in the so-called “unfavorable curvature” region.For the slab ITG (sITG) mode, different contributions of the parallel and perpendicular fluxes can result in negative effective compressibility. When the gradient in temperature surpasses the gradient indensity, a region can experience a net inflow of pressure, being carriedmost rapidly perpendicular to the magnetic field, while the net flux ofdensity is negative, being carried mostly rapidly parallel to the field.This allows local changes in density and pressure to be of oppositesign, and thus @pi @ni 0.In our model, stabilization of short wavelengths is provided byion inertia. This effect results in an ion sound Larmor radius term,which stabilizes the mode at a high k? qi (where k? is the perpendicular wavenumber and qi the ion gyroradius). In the present work, weignore the consideration of ion Landau damping, which can beincluded later using an appropriate closure relation.45,46A. Toroidal ITGWe begin our derivation of the tITG by considering the electrondynamics, which is governed by the electrostatic parallel electronmomentum balance in which inertia is ignored rk p e :0 ¼ ene rk /(1)Here, ne (later ni) is the electron (ion) number density, / is the potential, pe (later pi) is the electron (ion) pressure, and in the local model,27, 072507-2

Physics of PlasmasARTICLErk is the derivative along z. This approximation is valid for the lowfrequency regime, x kz vTe (where kz is the wave number along themagnetic field direction and vTe (later vTi ) the electron (ion) thermalvelocity), and treats the electrons as adiabatic, resulting in theBoltzmann relation e ¼n 1 e/n0 ;s Ti0(2)where s ¼ Te0 Ti0 is a constant based on centerline temperatures, n0is the equilibrium density, and the over-tilde is used to highlight someof the fluctuating quantities (mainly variables).To describe the ion dynamics, we utilize the ion continuity andenergy equations@niþ r ½ni ðvE þ v pi Þ þ ni0 r v Ii ¼ 0;@t@pi52þ vE rpi þ pi r ðvE þ vpi Þ þ r q ¼ 0;33@t(3)(4)where ni0 ¼ n0 is the equilibrium ion density, and q is the heat flux.The E B, diamagnetic, and inertial drifts arevE ¼ b r/;Bv pi ¼ b rpi;eni BvIi ¼1 d0b vE ;xcidtrespectively, and have been obtained from the perpendicular momentum equation assuming x xci (xci the ion cyclotron frequency)and d0 dt ¼ @ @t þ vE r is the convective derivative includingonly the E B drift. Note that we have used the Boussinesq approximation in (3), which omits density convection by the inertial drift.Furthermore, for the meantime, the diamagnetic drift is neglectedwithin the inertial drift, however, it will be included later in the nonlinear JOREK system via the gyroviscous cancelation.47 This system ofequations is supplemented using quasineutrality ne ¼ ni, and the diamagnetic heat flux,48 q ¼ qð0Þqð0Þ ¼5 pi b rTi :2 eB0(5)It is important to note that this closure includes only fluid effects, andthus the effect of ion Landau damping is not investigated at this time.This model for tITGs is quite similar to the Weiland model,39,40 whichhas been used extensively in ITG investigations.Using local slab coordinates (x radial, y poloidal, and z along B), i0 ; P ¼ p pi0 , the continuityand the normalized variables U ¼ e/ Tiand energy equations become1 @U@@U sq2s r2? U ¼ v i vDi rðU þ PÞ;s @t@t@y @P@U 51¼ ð1 þ gi Þv i vDi r U 1 þ 2P ;@t@y 3s(6)(7)where the magnetic curvature and gradient drift, and the equilibriumdiamagnetic drift are2Ti Ti @ðln ni0 Þ;ðb r ln BÞ; v i ¼eBeB @xrespectively, gi ¼ @ ln Ti0 @ ln ni0 , and q2s ¼ Te mi e2 B2 , is the ionsound Larmor radius. Here, the terms proportional to v i originatev Di ¼Phys. Plasmas 27, 072507 (2020); doi: 10.1063/5.0006765Published under license by AIP Publishingscitation.org/journal/phpfrom the E B advection of density and pressure, whereas those proportional to v Di originate from the divergence of the E B and diamagnetic drifts, along with r qð0Þ . Note that due to quasineutralityand our use of adiabatic electrons, the continuity equation nowappears as an equation for potential.Neglecting the inertial term, proportional to q2s , and consideringplane waves, this system can be reduced to the dispersion ��ffiffiffi5sx ¼ xDi þ ðxDi x i Þ6 sxDi x i ðgcr gi Þ;32(8)where xDi ¼ vDi ky ; x i ¼ v i ky (ky is the wavenumber in the y direction), and 2 s xDi x i10 xDiþ 2 þ:(9)gcr ¼ þ3 4 x i xDi9s x iThus, the tITG system is unstable for gi gcr .B. Slab ITGAs with the tITG, we consider adiabatic electrons for the sITG, asdefined in (2). We neglect the toroidal effects, provided by the diamagnetic drifts, and instability is now provided via parallel ion motion,which is included by incorporating the parallel ion velocity, v ik , andthe parallel ion momentum equationhi@ni þ ni r vIi ¼ 0;(10)þ r ni ðv E þ v ik bÞ@t@pi rpi þ 5 pi r ðv E þ v ik bÞ ¼ 0;þ ðv E þ v ik bÞ(11)3@t@ v ik rk p i : r (12)þ ðv E þ v ik bÞv ik ¼ eni rk /ni mi@tThe same simplifications for v E occur as in the tITG, and v ik isretained as a variable. Utilizing the normalized variable Vik ¼ v ik vTi ,along with U and P as before, results in the linearized system@Vik1 @U q2s @ 2@U r U ¼ v i vTi;s @t@ys @t ?@z@P@U 5 @Vik¼ ð1 þ gi Þv i vTi;@t@y 3@z@Vik@U@P vTi:¼ vTi@z@z@t(13)(14)(15)Although the tITG mechanism is easily visualized through theRayleigh–Taylor instability, the sITG is not so simple. Nevertheless,insight into how instability arises within the sITG can be found byexpanding the spatial components of these equations onto a Fourierbasis. This yields@U¼ iv i ky U ivTi kz Vik ;@t@P5¼ ið1 þ gi Þv i ky U i vTi kz Vik ;@t3@Vik¼ ivTi kz ðU þ PÞ;@t(16)(17)(18)where we have ignored inertia and consider s ¼ 1 for simplicity. For apositive perturbation in potential and pressure, i.e., a locally increased27, 072507-3

Physics of Plasmasion density, (18) shows that the ions will spread away along the fieldlines, departing from the considered region. The correspondingdecreases in density and pressure are encapsulated in the rightmostterms in (16) and (17). Note that these first two equations are nearlyidentical, except for the factors of ð1 þ gi Þ, and 5/3 in the pressureequation (17). For ð1 þ gi Þ ¼ 5 3 these two equations will yield differing rates of change in density and pressure, however, they can certainly not take on different signs (what is needed for a negativecompressibility instability). Considering, now, ð1 þ gi Þ 5 3, and allterms on the right hand side being positive, we see that the time rate ofchange in pressure can, in fact, take on a different sign than the timerate of change in density, which is proportional to U. This is the sourceof the instability.To see how this occurs practically, consider Fig. 1, in whichð1 þ gi Þ 5 3, and there is a positive perturbation in ion density(and subsequently potential and pressure). This perturbation results inthe formation of a E B vortex, owing to the first term on the righthand side of (16) and (17). On the lower part of this vortex, since thepressure gradient is larger than the density gradient, more pressurewill be convected into the area than density—to make this clear, wecan imagine the density gradient as vanishingly small, thus eliminatingthe density convection entirely. At the same time, density leaves theregion along the field line. Thus, the lower part of the vortex experiences increasing pressure (from the E B convection) and decreasingdensity (from the parallel motion)—negative compressibility. In thecase when the density gradient is finite, the situation simply amountsto the net change in pressure being positive, due to its E B convection into the region outpacing its parallel motion out of it, while thenet change in density is negative, due to its parallel motion outpacingits E B convection.The system of equations defined by (13)–(15) can be reduced tothe following dispersion relation:ARTICLEscitation.org/journal/php"# xx5x2þ x i 1 þ gi x i ¼ 0:x 2 2kk vTi s3s3(19)In the limit of kk vTi 1, the cubic term dominates, and we are leftwithx3 þ k2k v2Ti sðgi 2 3Þx i ¼ 0:(20)This is always unstable, provided gi ¼6 2 3, as is the case for commonfluid models of the sITG,41,42 unlike the critical thresholds for thesITG which arise in kinetic models.43,44C. Linear analysis of the general modelThe general ITG model is composed of both the tITG and sITGmechanisms, and is given byhi@ni þ ni0 r v Ii ¼ 0;(21)þ r ni ðvE þ v pi þ v ik bÞ@t@pi rpiþ ðvE þ v ik bÞ@t5 ¼ 2 r q;(22)þ pi r ðvE þ v pi þ v ik bÞ33@ v ik rk p i ; r (23)ni miþ ðvE þ v ik bÞv ik ¼ eni r/@talong with (1) and quasineutrality. The linearized system is given by@Vik1 @U q2s @ 2@U r U ¼ v i v Di rðU þ PÞ vTi;s @t@ys @t ?@z@P@U¼ ð1 þ gi Þv i@t@y 515 @Vik vDi r U 1 þ 2P vTi;3s3@z@Vik@U@P vTi:¼ vTi@z@z@t(24)(25)(26)For local analysis, we selectkz ¼ FIG. 1. Schematic of the sITG mechanism (see the text for description). Note,N ¼ n i ni0 U is only used in this image—U is the main variable.Phys. Plasmas 27, 072507 (2020); doi: 10.1063/5.0006765Published under license by AIP Publishingxky ;Ls(27) as is geometrically appropriate for slab coordinates (ikz ¼ b r),wherex is the distance from the rational surface, Ls ¼ Rq s is the magneticshear length, and s ¼ @ ln q @ ln r is the magnetic shear. We estimatean appropriate distance x, using a relevant simulation presented laterin Fig. 5(a). It will be shown that in this case (which is the most radially localized global mode), the mode appears directly between tworational surfaces, at a distance of x ¼ 1 cm from each of them.Since kx only appears within k? , and subsequently only leads tomode stabilization through the ion sound Larmour radius term, weutilize the value kx ¼ 0. With these values, along withq ¼ 1:4; R ¼ 171 cm; and s ¼ 1 2, the growth rate and frequency ofthe local tITG, sITG, and ITG models are shown in Fig. 2. We can seethat the general ITG behavior is dominated by the tITG, with somecompetition in mode growth49 arising from the slab mechanism. Thefrequency of the general mode is, however, larger than that of bothtITG and sITG, and increases for increasing magnetic shear. The magnetic shear only effects the sITG components (via the determination of27, 072507-4

Physics of PlasmasFIG. 2. Analytical growth rate and frequency of tITG, sITG, and ITG models.Normalized with vTi ¼ 4:48 105 m/s and Ln ¼ 76:6 cm.kz), and the sITG is stable as s ! 0, where the growth rate curve isflattened to larger ky qi , and for s 1:2, where the region of instabilityis pushed toward the origin. This behavior affects the ITG mode differently, for which the overall behavior approaches that of the tITG inthe limit s ! 0, whereas for s 1:3, the ITG is stabilized.III. THE JOREK CODEJOREK was designed to simulate large-scale MHD phenomena,such as edge-localized modes and disruptions in realistic tokamakgeometry.50–52 JOREK simulates the complete tokamak domain(including the scrape-off-layer) by using Fourier decomposition in thetoroidal direction, and the finite element method in the poloidal plane.Each element is decomposed on the basis of cubic Hermite polynomials, which allow for accurate representations of the variables andtheir second order derivatives on the scale of the grid size. For the purpose of this investigation, we have adapted the code to simulate ourtwo-fluid model at scales relevant to ITG troduced to match with the original benchmarking paper, Ref. 31.This paper used a circular modified equilibrium without theShafranov shift, so for fair comparison, we were obliged to adapt tothis situation, calculating first our equilibrium and the correspondingmesh by utilizing realistic and comparable values for the current profile, but imposing an artificially reduced pressure profile. It would bepossible to achieve the same result through careful and significantmodification of the current profile, however, this methodology wasnot pursued in depth, although several results which employed thismethod will be discussed below.For the duration of this paper, we use the same magnetic boundary conditions, and input profiles, and we vary only the value of Pf.This leaves the gi profile as a function of w or the minor radius, r, aswell as the gradients as a function of w unchanged. The gradients dochange slightly with respect to r due to variation of the Shafranov shiftwith different values of Pf, and the subsequent differences between theLow and High Field Sides (LFS and HFS, respectively). For this reason,all results will be reported with respect to w. The input profiles arebased on the well-known CYCLONE base case,31 where we matchedtheir values of LT ¼ 25 cm and s ¼ 0:78 where our modes are formed.A plot of the gi parameter, along with the magnetic shear for threeselect values of Pf ¼ 1; 1 2; 1 100, is shown in Fig. 3. These threevalues of shear will be used throughout the paper, and will be referredto as the low, medium, and high shear cases, respectively. The value ofq where the modes are typically formed for these three cases areq ¼ 1:1; 1:25; and 1:6, respectively. The profiles of temperature, density, and FF 0 are taken as w 0:321þ 0:03;(29)TðwÞ ¼ 0:0906tanh0:3kB l0 n0 w 0:32qðwÞ ¼ 0:5tanhþ 0:5 2mi n0 ;(30)0:5A. Equilibrium in JOREKBefore the numerical simulation of our ITG equations can begin,an equilibrium state must be defined. JOREK uses a built-in equilibrium solver which calculates a full MHD equilibrium, utilizing theGrad–Shafranov equation @ 1 @w@2wþ 2 ¼ Pf l0 R2 p0 ðwÞ l20 FðwÞF 0 ðwÞ:(28)R@R R @R@ZThis equation is represented in the poloidal plane with coordinates Rand Z, where the angular coordinate, u, represents the toroidal direction. The boundary conditions for the magnetic flux coordinate,w ¼ wðR; ZÞ, define the poloidal boundaries, and the pressure, pðwÞ,and toroidal field function, FðwÞ, are defined as inputs. For our purposes, the “pressure factor,” Pf 1 is utilized to modify the calculatedequilibrium by artificially reducing the pressure gradient. This reducesthe Shafranov shift and increases the safety factor and shear profiles byincreasing the relative significance of the FF 0 term. This factor wasPhys. Plasmas 27, 072507 (2020); doi: 10.1063/5.0006765Published under license by AIP PublishingFIG. 3. Shear profile for the low (Pf ¼ 1), medium (Pf ¼ 1 2), and high(Pf ¼ 1 100) shear cases, alongside the gi profile used for all simulations. Alsoshown are typical global profiles of / identifying where the mode is commonlylocated.27, 072507-5

Physics of PlasmasARTICLE 122FF ðwÞ ¼ 3:6ð1 1:8w þ w Þ þ sech ðw 10Þ2 1 1w 51; tanh2 20:03l20(dropping the 0 ), and using f ; g ¼ ð@ @RÞð@ @ZÞ ð@ @ZÞð@ @RÞ as the Poisson bracket, the system of equationsdefined by (1), and (32)–(35) becomes0(31) 19 3m .where w is normalized between 0 and 1, and n0 ¼ 6:1 10The input pðwÞ is simply the product of (29) and (30). Normalizationsare discussed in Sec. III B, and for the present investigation the equilibrium electrostatic potential is set to zero. In this paper, all our simulations are conducted in circular geometry with a minor radiusa ¼ 62:5 cm and a major radius R0 ¼ 170 cm.B. ITG model implemented in the JOREK codeFor the simulation in JOREK, we adapt the system of equationsfor the local ITG system, defined by (2) and (21)–(23). The JOREKmodel includes the tITG and sITG, and also includes parallel ionmotion in the convective derivatives. In addition, we retain the densityconvection by the ion inertial drift, v Ii rni , therefore forgoing theBoussinesq approximation.36 This is achieved by retaining the total niwithin the divergence of the (poloidal) vorticity term ni mir/:(32)X¼r ?eB2Utilizing the gyroviscous cancelation,47 with the convective derivative r in v Ii ¼ x 1 b defined as d dt ¼ @ @t þ ðvE þ vpi þ vik bÞci dvE dt, the inertial term becomes r ðni vIi Þ ¼ dX dt. Thus, thedensity, pressure, and parallel velocity equations used for the simulation are53d ¼ r ðD? rni Þ;ðni XÞ þ ni r ðvE þ v pi þ vik bÞdtd1 pi 5 þ pi r ðv E þ vpi þ vik bÞ3dt2 k rk Ti Þ;¼ r qð0Þ þ r ðK? r? Ti þ bK3d1 vik¼ eni rk / rk pi þ lk r2 vik ;nidt(33)Phys. Plasmas 27, 072507 (2020); doi: 10.1063/5.0006765sIC F0 @/ sIC nini f/; wg@uR2RF0 @niTe0þ sICþsIC Te0 2fni ; wg;R @uR@ ðni XÞR¼ sIC Rf/; ni Xg þ sIC fpi ; Xg@tniF0 @ðni XÞ 1 vik fni X; wg 2 vik@uRR@/@pi F0 @vikþ 2sIC niþ2sIC ni@Z@Z R2 @u1 ni vik ; w þ r ðD? rni Þ lX r2? X;R@piF0 @pi 1¼ sIC Rf/; pi g 2 vik vik fpi ; wg@tR@u R@/@piF0 @vikþ 2csIC Ti c 2 piþ 2csIC pi@Z@ZR@u1@Ti c pi vik ; w þ 2csIC piR@Z þ r K? r? Ti þ bKk rk Ti ;0¼ Bni@vikF0 @vik¼ Bni sIC R /; vik Bni 2 vik@tR@uBniF0 @/ 1 ni f/; wgvik vik ; w 2 niRR @u RF0 @pi 1 fpi ; wg þ Blk r2 vik ; 2R @u R(36)(37)(38)(39)respectively, where(34) X ¼ s2IC r ni R2 r? / ;(35)and pi ¼ ni Ti represents the pressure.For the duration of our analysis, we use F0 ¼ 3:247 T m, andunless otherwise specified, we use the value of sIC ¼ 0:01 T 1 . Our diffusivities are maintained at or near D? ¼ K? ¼ lk ¼ 2:5 10 5 T mand Kk ¼ 0:1 T m.JOREK uses Fourier decomposition (in the toroidal harmonic, n)to solve these equations in the toroidal direction. In all of our currentsimulations, we omit the evolution of the n ¼ 0 harmonic, whichmeans that the equilibrium is fixed. A large portion of the analysis presented in this paper is concerned with the behavior of single-n modes.We will refer to these simulations as linear, although a range of m harmonics will still be stimulated, and will indeed interact with each other(through cubic coupling discussed in Sec. VII A) when the perturbations reach significant amplitudes.where D? ; K? , and lk are small diffusion coefficients used to helpmaintain numerical stability, and Kk , is used to encapsulate paralleltemperature flux. The convective derivative in the latter two equations r. As before, the electron dynamicsis d1 dt ¼ @ @t þ ðv E þ vik bÞis adiabatic, and the governing equation is left in the form of (1).As with the Grad–Shafranov equilibrium, equations simulated inJOREK must be represented in cylindrical coordinates, with R and Zbeing the radial and vertical coordinates in the poloidal plane, and ubeing the toroidal coordinate. Normalizations are chosen to leave thepffiffiffiffiffiffiffiffiffiffispatial coordinates in meters, and the time, t ! q0 l0 t (q0 ¼ 2mi n0 ,with n0 ¼ 6:1 10 19 m 3 being the centerline density), in unitsrelated to the Alfv en time. The magnetic field is represented asB ¼ F0 ru þ rw ru, where F0 ¼ RBu , and w is the fixed fluxfunction, calculated via the Grad–S

We present global linear and nonlinear simulations of ion temperature gradient instabilities based on a fluid formulation, with an adapted . ITG model, and describe their linear analytic behavior. Following this, in Sec. III we briefly overview JOREK, the code we use for the simula-tion, outline our calculation ofequilibrium,anddescribe the .