Comparisons and Physics Basis of Tokamak Transport Models and Turbulence Simulations

A. M. Dimits,^{1} G. Bateman,^{2} M. A. Beer,^{3} B.I. Cohen,^{1} W. Dorland,^{4} G. W. Hammett,^{3} C. Kim,^{5} J. E. Kinsey,^{6} M. Kotschenreuther,^{7} A. H. Kritz,^{2} L. L. Lao,^{6} J. Mandrekas,^{8} W. M. Nevins,^{1} S. E. Parker,^{5} A. J. Redd,^{9} G. Rewoldt,^{2} D. E. Shumaker,^{1} R. Sydora,^{10} J. Weiland^{11}

1Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, CA 94550 2Lehigh University, Lehigh, PA 18015 3Princeton Plasma Physics Laboratory, Princeton University, Princeton, NJ 08540 4University of Maryland, College Park, MD 20742 5University of Colorado, Boulder, CO 80309 6General Atomics, Inc., San Diego, CA 92186-5608 7Institute for Fusion Studies, Univ. of Texas, Austin, TX 78712 8Georgia Institute of Technology, Atlanta, GA 30332-0225 9University of Washington, Seattle, WA 98195 10University of Alberta, Edmonton, Alberta, AB T6G2J1 Canada 11Chalmers University of Technology, S-412 96 Goteborg, Sweden

Abstract

The predictions of gyrokinetic and gyrofluid simulations of ion-temperature-gradient (ITG) instability and turbulence in tokamak plasmas as well as some tokamak plasma thermal transport models, which have been widely used for predicting the performance of the proposed ITER tokamak, are compared. These comparisons provide information on effects of differences in the physics content of the various models and on the fusion-relevant figures of merit of plasma performance predicted by the models. Many of the comparisons are undertaken for a simplified plasma model and geometry which is an idealization of the plasma conditions and geometry in a DIII-D H-mode experiment. Most of the models show good agreements in their predictions and assumptions for the linear growth rates and frequencies. There are some differences associated with different equilibria. However, there are significant differences in the transport levels between the models. The causes of some of the differences are examined in some detail, with particular attention to numerical convergence in the turbulence simulations (with respect to simulation mesh size, system size and, for particle-based simulations, the particle number). The implications for predictions of fusion plasma performance are also discussed.

PACS numbers: 52.35.Ra, 52.65.Tt,52.25.Fi 1. INTRODUCTION

We examine the physics basis and predictions of some tokamak plasma thermal transport models which have been widely used for predicting the performance of the proposed ITER tokamak [1]. This topic is of considerable importance and current interest since different models in use give conflicting predictions on whether ITER will or will not achieve thermonuclear ignition [2, 3, 4, 5, 6, 7].

The confinement in current tokamak experiments is generally believed to be degraded primarily by turbulence driven by “low-frequency” microinstabilities (instabilities occurring on drift timescales). The disagreements in Refs. [2]-[7] about predictions of ITER performance arise from using different models of drift-instability-driven transport, which can be distinguished in terms of their treatment of the detailed physics of microinstabilities.

In order to build confidence in predictions of tokamak plasma performance, the models used need to be tested. Here we rely on tests involving cross checks between a number of complementary models. A minimal condition for confidence in this case is agreement between the models where they can be compared under the same physical conditions. Where models disagree, a clear understanding of the reasons for the disagreements is important. This may take the form of knowledge that one of the models is being applied outside of the range of validity of a derivation that underlies it. Even when this agreement is achieved, there is generally no guarantee that the set of models will have predictive power, if applied outside the range of physical conditions for which they have been tested.

This evaluation process presents a great scientific challenge. These models are complex and contain many physical effects. The predictions of the models should be compared at various levels, not only at their final prediction, which is often a radial profile (e.g., density, flow velocity, temperature), but also at the level of comparison of the various submodels, theories and calculations.

In this paper we compare core transport predictions from the following: the IFS-PPPL [8, 9] and Multi-Mode (MMM) [10] transport models; flux-tube gyrofluid [11], flux-tube gyrokinetic [12], and global gyrokinetic [13] turbulence simulations; linear one-dimensional high-n eigenmode and initial-value calculations [14] and their use in transport model calculations. The MMM and IFS-PPPL models are two widely-used transport models for ITER predictions. Also represented here are some of the most advanced large-scale three-dimensional toroidal turbulence simulations, and the most widely used high-n linear and quasilinear calculations.

We first define some terms used here. The term “gyrokinetic” refers to a kinetic model (i.e., one which evolves functions of position-velocity phase-space variables and time), appropriate to charged particles in a strong magnetic field, in which a multipletimescale perturbation expansion in the ratio of the gyroperiod to the timescales of the phenomena of interest is made [15]. When appropriate, such a model is much more efficient than one that tracks the full particle dynamics (including the gyromotion). Gyrokinetic models retain “finite-gyroradius” effects (effects that arise when the scale of the gyro-orbit is comparable to the spatial scale of the phenomena of interest) nonperturbatively. This aspect distinguishes them from “drift-kinetic” models which either ignore or treat finite-gyroradius effects perturbatively. A “fluid” model of a gas or plasma is one that evolves “fluid variables” which are functions of position and time (and not particle velocity). “Gyrofluid” [11] models are a special class of fluid models derived from the gyrokinetic equations, which similarly have non-perturbative aspects to their treatments of finite-gyroradius effects. Those discussed here can be called “gyro-Landau-fluid” models in that they also contain (fluid) models of Landau damping and related processes.

We present detailed comparisons of the various models using a simplified physics problem which still contains the essence of core transport in an ITER-like discharge. A key point of comparison is the predictions of the various models for the χ_{i}, the ion thermal diffusivity (ion thermal flux divided by the ion temperature gradient). Disagreements are found between the predictions of different models for the ion thermal transport. We examine various possible reasons for the disagreement. Particular attention is paid to linear mode growth rates and frequencies, “zonal” flux-surface-averaged flow damping, and noise due to particle discreteness in the nonlinear gyrokinetic codes. Finally, the implications for ITER performance are discussed.

2. TEST PROBLEMS

We focus on the “Cyclone DIII-D base case parameter set” which represents local parameters from an ITER-relevant DIII-D shot (shot #81499) [16], at time t = 4000ms., and minor radius r = 0.5a, where a is the minor radius of the last closed flux surface. A concentric-circular-cross-section model equilibrium is used, with n_{i} = n_{e} and T_{e} = T_{i}, where n_{i} and n_{e} are the ion and electron densities and T_{e} and T_{i} are the electron and ion temperatures. The parameter values in dimensionless form are η_{i} ≡ L_{n}/L_{T} = 3.114, where L_{n} and L_{T} are respectively the density and temperature gradient scale lengths, magnetic “safety factor” q ≡ rB_{t}/RB_{p} = 1.4, where R is the major radius and B_{t} and B_{p} are the toroidal and poloidal magnetic field components, ŝ ≡ (r/q)dq/dr = 0.776–0.796 (some minor variations due to constraints in some codes), R/L_{T} = 6.92, and ϵ ≡ r/R = 0.18.

Additional simplifications which are made in all models discussed here are (1) electrostatic fluctuations, (2) the electrons are taken to be adiabatic, and (3) a single dynamical ion species (which represents the “bulk” ions) is used. Linear stability results from the comprehensive linear gyrokinetic code of Kotschenreuther and from the FULL code which include multiple-ion species and electromagnetic effects are also discussed [14]. These simplifications match those in the nonlinear gyrofluid simulations that underlie the IFS-PPPL model. Scans have been made varying the temperature gradient scale length while keeping other physical parameters fixed. Additionally, cases have been compared in which each one of ŝ and ϵ were set to zero, with the other parameters held to the DIII-D base case values. The first helps isolate differences in the way magnetic shear is treated. The various spatial representations used become very similar in the limit of zero magnetic shear. The second helps isolate the effect of linear damping of flux-surface-averaged poloidal flows. The shear associated with these flows is an important saturation mechanism and the physics of these modes is under continuing study (see further discussion below on radial mode damping).

We also included a case using parameters from TFTR L-mode shot 41309. This was a case used in 1994 in a code comparison within the Numerical Tokamak Project (NTP) [17]. These parameter values are η_{i} = 4.0, q = 2.4, ŝ = 1.5–1.6, R/L_{T} = 10., n_{i}/n_{e} = T_{i}/T_{e} = 1.0, and ϵ = 0.2057. Comparisons are also discussed in which the purely radial modes in this “TFTR L-mode NTP test case” are suppressed. The purpose of running these cases was to verify previous NTP comparisons that showed gyrofluid and gyrokinetic values of χ_{i} which differed by factors of 2 or less (the gyrofluid χ_{i} was higher and the experimentally determined value was higher still by a factor slightly less than 2). Another earlier gyrofluid/gyrokinetic comparison[18], which looked primarily at slab geometry, found that there was good agreement in the slab simulation χ_{i} though the saturation level of the RMS Φ differed by 40%. We also included a case using parameters from the published result of Dimits et al. [12], which had parameters similar to the DIII-D base case.

3. DESCRIPTION OF MODELS

Simulations of a flux-tube sub-domain of the torus which neglect global profile scale effects give a finite ion thermal diffusivity χ_{i} and turbulent correlation lengths much shorter than the scale of profile variation anticipated for ITER. A reasonable conclusion from this is that χ_{i} depends primarily on local plasma parameters, if ITERrelevant values of the normalized gyroradius ρ_{i}/L_{n,T} are used[11, 12, 19]. Motivated by this, all of the models discussed below with the exception of the global gyrokinetic simulations assume that the simulation domain covers only a thin radial extent compared to the plasma minor radius so that quantities such as the density scale length, (dn/dr/n)^{−1}, ω_{∗}, etc. are taken to be independent of minor radius.

The IFS-PPPL [8, 9] model is based on nonlinear gyrofluid simulations [11], which predict the fluctuation and thermal transport characteristics of toroidal iontemperature-gradient-driven (ITG) turbulence, along with comprehensive linear gyrokinetic ballooning calculations [14], which provide accurate growth rates, critical temperature gradients, and a quasilinear estimate of χ_{e}/χ_{i}. A key aspect of the IFS-PPPL model is an interpolation formula which parameterizes both the gyrofluid χ_{i}’s, and calculations of the critical temperature gradients and mixing-length predictions of χ_{i} from the more comprehensive linear gyrokinetic ballooning code. This linear code [14] has full velocity-space dynamics including resonances, trapped particles, Coulomb collisional pitch-angle diffusion, etc. This corrects the somewhat inaccurate critical temperature gradient and the neglect of non-adiabatic electron physics in the gyrofluid simulations that were used as a basis for the IFS-PPPL model.

The nonlinear gyrofluid simulations [11] that underlie the IFS-PPPL transport model use a “gyrofluid” reduction of the gyrokinetic equations. The resulting gyrofluid equations are evolved in toroidal field-line-following coordinates.[20] The gyrofluid equations include toroidal effects (e.g., magnetic curvature drive) and kinetic effects such as toroidal drift resonances, linear and nonlinear FLR orbit-averaging, and parallel wave-particle resonances, as well as nonlinearly generated, fine-scale (k_{r}ρ_{i} ∼ 0.1), sheared poloidal flows, which play a major role in determining the saturation level for the turbulence. In this paper, we focus on gyrofluid simulations with an adiabatic electron response, since this was used in the simulations on which the IFS-PPPL model is based. Gyrofluid simulation codes with bounce-averaged nonadiabatic electrons have since been developed and exercised [21].

The philosophy underlying the Multi-mode transport model (MMM) [10] has been to utilize a collection of theoretically-derived transport models to predict temperature and density profiles in tokamak plasmas and adjust the models as necessary. As the models have improved, less adjustment has been needed to fit the experimental data. The 1995 Multi-mode model combines the fluxes predicted by the Weiland ITG/TEM model [22, 23] with those predicted by the Guzdar-Drake resistive ballooning model [24], with smaller contributions from kinetic ballooning modes and neoclassical transport [25]. The 1995 Multi-mode model was calibrated against a small number of experimental discharges, and then held fixed for all subsequent transport studies [26].

In the simulation of experimental plasmas, the transport predicted by the 1995 Multi-mode model is typically dominated by the contribution from the Weiland ITG/TEM model. This is the only contribution that is kept in the comparisons for the idealized Cyclone parameters in Figs.(1) and (3). The Weiland ITG model is based on a fluid description in which all moments that are driven by sources (i.e. fueling, heating) are included self-consistently. The fluid moments that are not driven by sources generally decay to zero. The model allows free energy exchange between different transport channels, leading to pinch fluxes. The transport coefficients are derived by using quasilinear theory and a mixing-length rule for saturation, which takes k_{θ}ρ_{s} = 0.316 [10]. The transport coefficients therefore have gyroBohm scaling. However, they have been found to agree well with some non-gyro-Bohm L-mode and H-mode experimental data [10, 25, 26]. The Weiland model also includes effects from the impurity profiles[10], fast ions, and T_{e} ≠ T_{i}, and has been extended to include parallel ion motion and electromagnetic effects [23].

The IFS-PPPL and MMM models are both basically gyroBohm-scaling models, though non-gyroBohm scalings can enter in several ways, for example, marginal stability connections to edge boundary conditions, change in particle fueling profiles or density profiles. Recent versions of the IFS-PPPL model (and the related GLF23 model[7]) also add stabilizing E×B shear, which can introduce additional non-gyroBohm scaling effects [4, 6, 7, 27, 28, 29].

The gyrokinetic simulation codes, both flux-tube [12] and global [13], solve the gyrokinetic Vlasov-Poisson system of equations [15, 30] (in the electrostatic limit) using particle-based δf methods [31, 32, 33] and “four-point gyroaveraging” [30]. A single fully toroidal nonlinear gyrokinetic ion species with equilibrium temperature, density, and velocity gradients is used. Adiabatic electrons with a zero response to the flux-surface-averaged potential [34, 35, 36] are used in the present comparisons, both in the flux-tube gyrokinetic and gyrofluid codes. The lack of response of the electrons to the flux-surface-averaged potentials is a key factor in the amplitudes to which these zonal flow modes are driven and therefore to the levels of turbulence and transport seen in the gyrofluid [35, 36] and gyrokinetic [34, 37] simulations. A low-β concentric-circular-cross-section model equilibrium is used here.

The flux-tube gyrokinetic simulations (like the gyrofluid simulations) use a fluxtube domain (bounded by four magnetic field lines) of small perpendicular extent, which spans one or more poloidal circuits in the parallel direction. The flux tube is taken to be periodic in the toroidal direction, and periodic in the radial direction with a toroidal offset such that a magnetic flux sheet is continuous across the radial boundary. This prevents saturation of the turbulence by profile relaxation. The field quantities are defined on a quasiballooning-coordinate grid [38]. The radial differences, interpolation, deposition, and smoothing are formed using shapes in configuration space (not ballooning-coordinate space) that are independent of poloidal location [38]. This combination of coordinates and shapes prevents grid stretching and resolution loss in the presence of magnetic and velocity shear, and allows a smooth implementation of the toroidal periodicity condition across the parallel boundary for arbitrary profiles of the magnetic safety factor q(r). It contrasts with the flux-tube gyrofluid code which uses direct discretization in ballooning coordinates.

The global gyrokinetic simulations [13] typically use a domain which spans the whole tokamak volume, although simulations using only an inner subvolume or a toroidal annulus have been reported [39]. The field quantities are represented on a radial-poloidal-toroidal mesh, and a fully nonlinear form of the gyrokinetic equations is solved [40] instead of the partially linearized form of ref. [15]. The important physics that the global codes allow for is the full radial variation of gradient quantities (e.g., temperature and density gradients, magnetic shear, etc.). These are generally a stabilizing effect, but get weaker in larger tokamaks with larger a/ρ[39, 41, 42]. In global gyrokinetic simulations to date the simulation domain is bounded and no explicit modeling of particle or thermal sources sinks is used. Thus profile relaxation, including kinetic profile relaxation, in which spatial gradients of ions in subregions of velocity space relax [34], may occur.

4. LINEAR COMPARISONS

Figure (1) shows linear frequencies and growth rates as a function of k_{θ}ρ_{i} obtained from several independent linear and nonlinear codes (our convention is that the thermal gyroradius and thermal speed are defined as ρ_{i} = v_{ti}/Ω_{ci} and v_{ti} = \sqrt{T_{i}/m_{i}}). Represented are the linear gyrokinetic codes of Kotschenreuther and of Rewoldt (Full code), the nonlinear gyrofluid code of Beer and coworkers, and the nonlinear gyrokinetic codes of Dimits (flux tube) and Sydora (global), as well as the fluid code of Weiland. Very good agreement between the various codes is observed. Agreement at this level is an important cross check of the codes since all of the gyrokinetic codes should have the same linear physics (with the exception that the global gyrokinetic code has radial variation of profile quantities), and since the closure in the gyrofluid code is designed to reproduce the gyrokinetic linear responses accurately.

Note that the “GK (Dimits)” gyrokinetic results are from a single flux-tube simulation containing many growing modes. The growth rate for the k_{θ}ρ_{i} = 0.1 mode has quite a large uncertainty due to several possible known effects. Since particles in the particle-based simulations constitute a structure that does not have the same periodicity as the mesh, there is linear coupling between different modes (both different toroidal mode numbers and different radial wavenumbers or ballooning angles).

[This linear coupling is relatively small and gets weaker as more particles are added.] This mode has slow time variation, i.e., both the expected ω and γ are small. It is therefore slow to reach its time-asymptotic linear behavior and is also susceptible to the above effects. For example, at a given k_{θ} there are many different θ_{0} modes (θ_{0} is the ballooning parameter value where k_{r} = 0) in the simulation growing simultaneously, and it can take time for the fastest growing mode to sufficiently dominate to get an accurate growth rate.

Whether the small differences in the growth rates shown in Fig. (2) are important to the transport depends on the dynamics of the saturated state. If the longer wavelength modes are primarily driven nonlinearly, then these differences are likely to be unimportant. If the correlation time of the long-wavelength modes is set by their linear growth rates (as is implicit in γ/k^{2}-type mixing laws) which assumes that the nonlinear driving does not set the correlation times, then these differences may be significant since the longer wavelengths dominate if one maximizes γ/k^{2}-type mixing formulas over wavenumber.

Additional points of agreement between the flux-tube nonlinear gyrokinetic code and the linear gyrokinetic code include the linear critical temperature gradients for both the DIII-D base case parameters and for ϵ = 0 but with the other parameters as for the DIII-D base case.

A second linear test is based on the linear damping of purely radial modes of the electrostatic potential, i.e., modes which have no variation within a flux surface. A linear theory for the residual levels of these modes in the collisionless limit has been given by Rosenbluth and Hinton [43]. In this benchmark, the gyrokinetic code is initialized with zero particle weights. A radially sinusoidal potential with no variation within the flux surfaces, which represents a near poloidal E×B flow, is imposed. The particle weights evolve, resulting first in a period of geodesic acoustic oscillations which eventually damp. At late time, the net electrostatic potential is less than the imposed potential. For a circular cross-section equilibrium in the large aspect-ratio limit, the theoretical prediction for the ratio of the late-time net potential to the initial residual level is given as a function 0.6h/(1.0 + 0.6h) of the single parameter h ≡ \sqrt{ϵ}/q^{2} [43].

Figure (2) shows the fractional residual E×B flows for two scans done with the flux-tube gyrokinetic code. In one scan, q is varied, while in the other scan ϵ is varied. Very good agreement is observed, lending confidence to both the Rosenbluth-Hinton theory and the simulation. The gyrofluid code, with the “thesis closure,” [21] which was used in the simulations that underlie the IFS/PPPL formulas used to make the ITER projections of Refs. [2], appears to give E×B flows consistent with the poloidal flow damping to zero [21], i.e., in the case of the setup here, the E×B flow damps essentially to zero, although the short-time response, on the ion-transit timescale agrees reasonably well with gyrokinetic calculations. In more recent work[44], improved gyrofluid closures have been developed for which the long-time response agrees better with Ref. [43], and are discussed in more detail in the next section.

5. NONLINEAR COMPARISONS

Figure (3) shows predictions from the various simulations and models for χ_{i} vs. R/L_{T} for a scan about the DIII-D base case parameters.

Two sets of flux-tube gyrokinetic results are shown. The agreement between these two sets is very good. The “U. Col. Boulder” results were run with 4 × 10^{6} particles, while the LLNL results were obtained with 8–34 × 10^{6} particles. This accounts for the difference at lower values of R/L_{T}. Specifically, the LLNL code run with 4×10^{6} particles also shows nonzero values of χ_{i} which are in agreement with the U. Col. Boulder code for R/L_{T} < 6.9. A remarkably good fit to LLNL gyrokinetic results is given by [45]

χ_{i}L_{n}/(ρ_{i}^{2}v_{ti}) ≃ 15.4[1.0 − 6.0(L_{T}/R)],  (1)

which is shown in Fig. (3). This fit corresponds to an offset linear dependence of the thermal flux on the temperature gradient Q ∝ (R/L_{T}−R/L_{T eff}). Note that R/L_{T eff} ≃ 6.0 > R/L_{T crit} ≃ 4.0, where L_{T crit} is the temperature gradient scale length at which the simulation is linearly marginally stable. This is a newly observed phenomenon and is likely associated (see the discussion of ϵ scans below) with undamped Rosenbluth-Hinton zonal E×B flows. [43]

The linear critical temperature gradient R/L_{T crit} has been checked by several different codes, and is found to be around 4.0 for this model s − α equilibrium in both the flux-tube gyrokinetic particle code and in Kotschenreuther’s gyrokinetic code. It is known that the quantitative details of linear growth rates and critical gradients can be fairly different in β = 0 s − α model equilibrium (which makes large aspect ratio approximations and has concentric circular flux surfaces) versus a realistic numerically-calculated equilibrium with finite aspect ratio (where even at zero β there is still a Shafranov-shift of the flux-surfaces due to the plasma current)[46, 47]. For the β = 0 s − α equilibrium used in these comparisons, the Full code gives a critical temperature gradient of about 3.7, while for a realistic numerical equilibrium the critical gradients drops to 2.5[47].

To be consistent, the 1994 IAEA version of the IFS-PPPL model[8] is used to compare with the various simulations in Fig.(3) because it also used adiabatic electrons in its parameterization of χ_{i} and R/L_{T crit} (the 1995 version of the IFS-PPPL model in Ref.[9] included linear gyrokinetic estimates of the effects of trapped-electrons on χ_{i} and R/L_{T crit} as well as on χ_{e}). The 1994 IFS-PPPL model predicts R/L_{T crit} = 3.1 for this case, but that is low because the IFS-PPPL model was based on Kotschenreuther’s linear gyrokinetic code using a more realistic equilibrium than the s − α model used for the other codes in Fig (3). The IFS-PPPL model was constructed to have a form χ = W_{NL}^{GF} D_{mixing}^{GK} ∝ G(R/L_{T} − R/L_{T crit}), and so by construction should go to zero at R/L_{T crit}. Thus for a more consistent comparison with simulations that assume a simpler s − α equilibrium, the IFS-PPPL curve in Fig.(3) was shifted to the right from R/L_{T crit} = 3.1 to R/L_{T crit} = 4. Part of the reason the IFS-PPPL model was constructed in this way was to attempt to correct for some known inaccuracies in the gyrofluid prediction of the critical gradient and the growth rates near marginal stability, where a slow growing residual mode can sometimes persist in gyrofluid simulations below the gyrokinetic critical gradient[48, 49]. [Improvements to the gyrofluid closures to include frequency dependence or nonlinear effects, along the lines suggested by Chang and Callen[50] or Mattor and Parker[51], may help with this problem.] Thus the IFS-PPPL curve is actually somewhat below the gyrofluid simulations in Fig.(3) using the 1994 gyrofluid closures[21]. The 1994-closure gyrofluid simulations also do not show the roll-over at high R/L_{T crit} that is in the IFS-PPPL model, which would presumably be more apparent at other parameters[28].

For the DIII-D base case parameters (at R/L_{T} = 6.9), the IFS-PPPL model χ_{i} is a factor of 2.7 larger than the gyrokinetic flux-tube results, and the 1994 gyrofluid simulations are a factor of 3.3 higher. Also in Fig.(3) are simulations labeled “98 PPPL GFL” that employ a recent neoclassical improvement of the gyrofluid closure[44]. This improved closure reduces the gyrofluid/gyrokinetic difference to a factor of 2 for the base case, and it is able to reproduce some of the nonlinear upshift in the effective critical gradient seen by the gyrokinetic code. This improved gyrofluid closure, and the possibility of further improved closures, will be discussed more in the next section. The differences between the χ_{i} vs. R/L_{T} curves for the flux-tube gyrokinetic code and the IFS/PPPL model can be characterized partly as a shift in R/L_{T eff}, which is a strictly linear value in the case of the IFS-PPPL model, and partly that χ_{i} shows a more gradual increase as the critical gradient is exceeded in the case of the gyrokinetic model. The flux-tube gyrofluid/gyrokinetic differences in Fig.(3) correspond to a 20-33% difference in the local temperature gradient predicted at fixed heating power, as discussed in more detail in Sec.(7).

The MMM model result (labeled “Weiland QL-ITG”) agrees with the flux-tube gyrokinetic result quite closely for the base case, although comparisons have not yet been carried out for a wider range of parameters. The MMM model also gives a reasonable prediction for the linear R/L_{T crit}, though both the MMM and the IFS-PPPL model miss the nonlinear increase in the effective R/L_{T crit} observed in the gyrokinetic simulations.

The MMM model and the fit to the gyrofluid simulations have a linear scaling for the transport with R/L_{T} − R/L_{T crit}, whereas the gyrokinetic simulations are reasonably well fit by a square-root dependence on R/L_{T} − R/L_{T eff} which is of the same form as the IFS/PPPL model. However, the offset-linear Q vs. R/L_{T} fit to the LLNL gyrokinetic results is a better fit than the best power-law fit. The MMM model result agrees with the gyrokinetic results for the base case parameters. The MMM model also gives a reasonable prediction for the linear R/L_{T crit}.

The global gyrokinetic results are 2.4 times lower than the flux-tube gyrokinetic results for the base case. The global code used the same local dimensionless parameters, a realistic temperature profile, and a value of ρ/a = 1/160 that is comparable to the actual DIII-D experiment but is somewhat large compared to values achieved in larger tokamaks such as JET or TFTR, and even larger compared to proposed designs such as ITER. Variations of the simulated tokamak size in global gyrokinetic simulations [39, 41, 42] show that χ_{i}/χ_{gyroBohm} increases as the simulation is made larger. This is consistent with a theoretical picture that radial variation in various profile and profile-gradient quantities (e.g., the diamagnetic velocity) introduces stabilizing effects [7, 29] that get weaker in larger tokamaks. For very large a/ρ, the global gyrokinetic simulations are expected to asymptote to the value of χ_{i} given by the flux-tube gyrokinetic simulations. Another difference in the global simulations is that the purely radial mode is coherent, fairly stationary, and has a radial scale comparable to the minor radius, whereas in the flux-tube simulations, the radial modes are at shorter wavelengths (smaller than the box size). Recent studies indicate that this difference will also go away as the global simulation domain is made larger and the profile variation weaker [39, 52].

In the DIII-D experiment that the parameters for Fig. (3) were based on, the measured χ_{i} = 0.16 (in the units of Fig. (3)) and the measured R/L_{T} = 6.9. This is very low compared to all of the simulations in Fig. (3). Although this is consistent with the general picture of strong ITG turbulence forcing the plasma to be near marginal stability,[2] one cannot draw this conclusion based solely on the simulations presented here. This is because the primary purpose of Fig. (3) is to compare different simulations with as similar a set of assumptions as possible, so a number of factors that are important in experiments are not included (such as non-adiabatic electrons, equilibrium rotation, and realistic geometry). Several of our simulation methods have been used to study equilibrium-scale sheared rotation[19, 37], which can be particularly important in DIII-D because of its unidirectional beam injection and resulting high toroidal rotation speeds. For example, while the IFS-PPPL model looks pessimistic compared to the experimental measurement in Fig. (3), applying the full IFS-PPPL model[28] including a model of the stabilizing influence of equilibrium-scale E×B flows[19, 53] gives a predicted central ion temperature which is actually somewhat above the measured temperature. However, there are some quantitative uncertainties in the standard models of stabilization due to E×B flow that lead to uncertainties in the predicted temperature profiles of order 10-30%, comparable in magnitude to the 20-33% differences in the gyrofluid/gyrokinetic temperature gradients described above. Some of these effects are hard to distinguish: comparable levels of agreement with experiments can be obtained with modified transport models where the magnitude of the E×B flows is reduced while the the coefficient of χ_{i} is simultaneously reduced[6] (the two effects offset each other somewhat). Nevertheless, there are a wide range of experiments indicating the general importance of equilibrium-scale E×B flows on the transport[27, 54, 55, 56].

6. POSSIBLE REASONS FOR THE DISAGREEMENTS

We now address possible reasons for the disagreements between the various nonlinear results discussed above. We focus first on the differences between the GK/GF flux-tube simulations. This is important since the gyrofluid models are an attempt to approximate the gyrokinetic equations for the problems of interest here, and the physical and numerical parameters have been otherwise matched in the comparisons. The differences between global and flux-tube simulations have already been discussed.

There are many things that might, in the absence of concrete data, be viewed as possible causes for the differences found between the GK/GF flux-tube codes.

These include differences in linear growth rates and critical gradients, differences in the linear damping rates or residual levels of the radial modes, nonconvergence with respect to system size and grid size, nonconvergence with respect to particle number (which may lead to excessive particle noise either in the radial modes or in the modes that have finite poloidal or toroidal variation), and nonlinear wave-particle effects.

The differences in linear growth rates and critical gradients have been addressed above for the Cyclone DIII-D parameters and are probably not important.

6.1 Causes of gyrofluid discrepancies, and improved gyrofluid closures.

As noted above, the original closures used in the gyrofluid simulations which underlie the predictions of Ref. [2] damp the poloidal flow to zero in most cases, and therefore do not properly reproduce the long-time residual levels predicted by Rosenbluth and Hinton [43]. This is probably the main cause of the gyrofluid/gyrokinetic differences near marginal stability, and may account for about half of the differences in stronger turbulence regimes. But there are cases where the gyrofluid/gyrokinetic differences do not appear to be attributable to the differences in the radial mode linear dynamics. The dependence of χ_{i} on r/R from gyrofluid and gyrokinetic simulations is similar, as shown in Fig. (4). The residual Rosenbluth-Hinton (RH) zonal flows vanish in the limit ϵ = r/R = 0, so any differences in that limit can not be attributed to those flows. However, the trapped ion drive of the turbulence also varies with r/R, and so this does not completely isolate just the effect of undamped flows. [The simulations in Fig. (4) are for the TFTR-based NTP test case described in Sec. (2), and show a smaller discrepancy, of about a factor of 2, than the Cyclone base case.] However there is evidence that the Rosenbluth-Hinton undamped component of the zonal flows is a significant part of gyrofluid errors, particularly near marginal stability. The R/L_{T} flux-tube gyrokinetic simulation scan at ϵ = 0 in Ref. [45] showed that χ_{i} becomes non-zero once R/L_{T} becomes slightly larger than the linear ϵ = 0 critical value. This contrasts with the scan done for ϵ = 0.18 (and therefore finite h) in Fig. (3), which indicates that the undamped Rosenbluth-Hinton flows play a role in the departure of R/L_{T eff} from R/L_{T crit}. This is further supported by the observation that in the cases where R/L_{T crit} < R/L_{T} < R/L_{T eff}, radial mode potentials develop stationary structures in which the peak shearing rates are significantly greater than the growth rate of the fastest growing ITG modes. Thus, the undamped Rosenbluth-Hinton flows significantly affect the behavior of χ_{i} near marginal stability.

The Rosenbluth-Hinton (RH) component of the zonal flows are linearly undamped except by collisions. The fact that a non-zero χ_{i} is observed in these collisionless gyrokinetic simulations for R/L_{T} > R/L_{T eff} is an indication that nonlinear damping of the RH zonal flows by turbulent viscosity is able to balance the nonlinear drive of these flows. One might expect that the turbulent viscosity would increase as R/L_{T} increases, so that the RH zonal flows would become unimportant relative to the other components of zonal flows when the turbulence is sufficiently strong that the turbulent damping rate of the RH component of the zonal flows becomes comparable to the damping rate due to collisionless transit-time magnetic pumping that affects the other components of zonal flows. However, certain types of turbulence can exhibit inverse cascades, and these issues warrant more study.

As was noted above, more recent nonlinear gyrofluid simulations have been completed using improved closures[44] that do allow for levels of RH undamped zonal flows in rough agreement with Ref. [43]. As seen in Fig. (3), these simulations show a nonlinear upshift in the effective critical temperature gradient, though not yet as large as the upshift in the gyrokinetic simulations. The improved closure used at present is able to match the RH residual flow to within 20% at k_{r}ρ ∼ 0.2, where the dominant contribution to the effective shearing-rate[57] is usually made, but the residual flow is about a factor of 2 low at very low k_{r}ρ. Future work will investigate further improvements of the neoclassical treatment of the gyrofluid closures to better match the RH residual flows, which should bring the gyrofluid simulations into better agreement with the gyrokinetic simulations in Fig. (3). Other improvements to the gyrofluid model to be investigated include frequency-dependent closures[50, 51]. This may be particularly helpful in improving the approximation of the branch cut in the toroidal kinetic response function[48, 49] and in improving the calculation of the linear critical gradient and growth or damping rates of various modes. An improved frequency-dependent gyrofluid closure may account for most of the remaining difference between the gyrofluid and gyrokinetic simulations in these collisionless-ion adiabatic-electron comparisons.

6.2 Nonconvergence with respect to system size and grid size

We have investigated and demonstrated convergence with respect to system size and grid size for the flux-tube gyrokinetic and gyrofluid simulations. For the NTP test-case parameters, this issue has been addressed in the flux-tube gyrokinetic simulations [12]; and the parameters for the DIII-D base case are similar. It was found that an increase in system size in the parallel direction made essentially no change in the simulation flux, while there was only a very weak dependence on perpendicular system size. Similarly, it has been verified, both for the NTP test-case and Cyclone DIII-D base-case parameters, that the parallel and perpendicular grid sizes are adequate. All of the simulations for the R/L_{T} scans in Fig.(3) were done at fixed system size. Some theories suggest that the turbulence shifts to longer wavelengths at higher R/L_{T}, and so it may be useful to eventually redo the high R/L_{T} simulations in a larger box size, which might cause the χ_{i} to increase some. The gyrofluid simulations have also been tested for and appear to be converged with respect to system size and grid size for the Cyclone DIII-D base-case parameters, though more studies could perhaps be done. A recent port of the gyrofluid code to the massively parallel T3E computer will allow convergence checks at significantly higher resolution. The comparisons made here between the gyrokinetic and gyrofluid simulations were done at similar system sizes, so this is very unlikely to be the cause of the difference. The use of direct discretization in ballooning coordinates in the gyrofluid simulations requires finer grid cells in the radial direction than in the other perpendicular direction if magnetic shear is present. Because of this and the fact that the mesh is explicitly involved in representing the advection of the fluid fields, establishing convergence with respect to grid size in the gyrofluid simulations has been found to be more subtle and to impose more stringent limits on the grid sizes in some cases than previously thought. Dimits has proposed an algorithm based on the periodicity of discrete Fourier transforms that could help reduce the resolution requirements in the gyrofluid code.

Next, we examine noise in the flux-tube gyrokinetic simulations, and address the possibility of nonconvergence with respect to particle number. Fig. (5) shows χ_{i} vs. time from a particle number scan for the base case. The simulations are for different numbers of particles ranging from 5 × 10^{5} to 1.34 × 10^{8}, corresponding to 1 to 128 particles per grid cell. [These simulations used a 128 × 128 × 32 (radial, poloidal, and parallel) grid, with finite-size particle filtering to smooth out fields with k_{i} >∼ 1/Δ_{i}, where Δ_{i} is the grid spacing in the i’th direction.] For 10^{6} or more particles, χ_{i} at late time does not appear to change significantly with particle number. The primary conclusion is therefore that χ_{i} appears to be converged with respect to particle number for 2–4 or more particles per grid cell for the base case parameters. There is some random variation in the late time-averaged χ_{i} for the different cases. When these random variations in χ_{i} (and somewhat larger random variations in the volume-averaged ϕ^{2}) were observed in initial convergence studies over a smaller range of particle number, questions were raised that motivated convergence studies to the very large particle number shown here, and motivated the additional scrambling tests described below. The particle number has now been varied over such a wide range, and the random variations in χ_{i} are sufficiently small and show no systematic dependence on particle number, that particle convergence does not appear to be a problem. This conclusion is made even more convincing by the scrambling tests described below. These random variations are presumably just due to the sensitive dependence on initial conditions of a chaotic system with long time-scale dynamics (for example, interactions with low k modes or zonal flows). These small random variations should average out over longer times or multiple realizations. [The somewhat larger variations in volume-averaged ϕ^{2} are presumably due to similar effects. The m = n = 0 zonal component of ϕ can have quite large amplitudes at small k_{r}, but have little physical consequence because their resulting shearing rate ∝ k_{r}^{2}ϕ_{k} is very small.] There is also some increase in the level of the initial peak in χ_{i} which persists even if the scan is done by increasing the initial weights (as the square root of the particle number) so as to keep the initial mean noise level fixed. The 5 × 10^{5}-particle case, which corresponds to one particle per grid cell, shows secular growth in χ_{i} beyond Time = 700. This is probably due to a noise-driven runaway process in which the RMS average particle weight, related to the detailed δf-particle entropy, increases with the time integral of χ_{i}. The noise causes thermal transport (χ_{i}), both of which increase together.

In order to further assess the impact of particle discreteness, the following scrambling test [58] of the noise level was performed. The gyrokinetic code was run saving restart files at selected times. New restart files were formed from these by scrambling the particle weight list. The gyrokinetic code was restarted from these scrambled restart files. After the restart, the temperature gradient was reduced to slightly below the linear marginally stable value in order to eliminate unstable ITG modes. The test was done using 8 × 10^{6} particles in the simulation.

Once the gyrokinetic code has run in the nonlinear phase longer than a characteristic eddy turnover time, a typical simulation particle has moved from its initial position farther than a characteristic eddy radial scale. Eventually, simulation particles nearby in the (x, v_{∥}, μ) phase space have weights that have large uncorrelated components from which a physical density field must be estimated. Physically this represents the failure of the particle representation to resolve δf. The “worst-case” interpretation is that this uncorrelated component is “noise,” although further work is needed to to establish a clear interpretation. The scrambling of the particle weight list eliminates the physical signal, leaves a state with a similar level of the uncorrelated component to the pre-scrambling state, and therefore provides a measure of this uncorrelated component. The post-scrambling restarted simulations were run long enough for stable geodesic acoustic fluctuations, which are present immediately after the scrambling, to damp. The resulting electrostatic potentials (or the shearing rates derived from them) provide measures of the uncorrelated component of the signal in the gyrokinetic simulation just prior to the scrambling.

Shown in Fig. (6) are the time histories of χ_{i} and the mean squared E×B shearing rate (L_{T i}S/v_{ti})^{2} associated with the flux-surface-averaged electrostatic potential, both in the absence of scrambling and when the scrambling and gradient reduction is done at three times during the run. Both of these quantities decrease after the scrambling. The relative reduction is less the later the scrambling is done, indicating a gradual buildup of noise. However, even at the latest time, the post scrambling values are down by an order of magnitude. This indicates that the relative impact of noise is small (or at most moderate at the latest time), and supports the conclusion that the simulations are converged with respect to particle number. The noise shearing rate for 2 million particles (4 particles per cell) would be a factor of 4 larger than shown in Fig. 6, and thus at late time would be comparable to the pre-scrambling signal. But much of that noise shearing is at high k_{r} and fluctuates rapidly in time, so it is less effective than shearing by low k_{r} modes[57] and can be ignored.

The numerical convergence of the particle codes with respect to particles does seem to be sensitive to how close the system is to marginal stability, however. For values of R/L_{T} somewhat lower than the Cyclone DIII-D base case value, but well above the linear marginal value of 4.0 (e.g., R/L_{T} = 5.3 and 6.0), after the linear growth and nonlinear saturation phases, the system evolves to stable states which have radially dependent flux-surface-averaged temperature gradients and E×B flows. In these situations the radial thermal flux asymptotes to zero. 64 particles per cell are needed to converge to this result for R/L_{T} = 6.0. Evidently, the stable nonlinear states become quite delicate as a threshold value (larger than the linear critical value) of the volume averaged temperature gradient is approached from below.

7. SENSITIVITY OF PREDICTED TEMPERATURE PROFILES AND FUSION GAIN TO TRANSPORT MODEL VARIATIONS

The differences between the χ_{i} vs. R/L_{T} curves for the flux-tube gyrokinetic code and the IFS/PPPL model in Fig.(3) can be characterized partly as a shift in R/L_{T eff}, which is a strictly linear value in the case of the IFS-PPPL model, and partly that χ_{i} shows a more gradual increase as the critical gradient is exceeded in the case of the gyrokinetic model. It is important to note that these differences cannot be characterized by a simple ratio, or multiplication factor. For example, lowering R/L_{T} just a bit from 6.9 to 6 (where the gyrokinetic simulations vanish) causes the relative error to become a factor of infinity. Rather than compare the χ_{i}’s at a fixed R/L_{T}, one can instead turn Fig.(3) around and compare the predicted temperature gradient at a fixed amount of heating power. This way of characterizing the gyrofluid/gyrokinetic differences is more relevant to experiments. When heating power is added to a plasma, the temperature gradient on every flux surface will rise until the resulting χ_{i} is large enough to balance the rate at which the plasma is being heated. The heat flux (or power flow) across a given magnetic surface is given by P = −nAχ∇T in circular geometry, where n is the density and A is the surface area. This can be written as

P = nA\frac{ρ_{i}^{2}v_{ti}}{L_{n}} \frac{T}{R} χ̂ (\frac{R}{L_{T}} ) \frac{R}{L_{T}} = P_{0} g (\frac{R}{L_{T}} ) .  (2)

where χ̂(R/L_{T}) is the normalized χ_{i} in the units of Fig.(3) as a function of R/L_{T}, and g(x) = χ̂(x)x. For a given amount of normalized heating power P/P_{0}, one can then solve this equation to find the resulting temperature gradient R/L_{T} = g^{−1}(P/P_{0}). Carrying this out for the IFS-PPPL curve in Fig.(3), and for the gyrokinetic fluxtube results in Fig.3 (and Eq.(1)), we then take the ratio of these two predicted temperature gradients to measure the relative error. This ratio is plotted vs. normalized heating power P/P_{0} in Fig.(7). At low heating power, the temperature gradient will be close to marginal stability. In this limit in Fig.(3), the IFS-PPPL linear critical gradient R/L_{T crit} = 4, which is 33% low compared to the gyrokinetic nonlinear critical gradient of 6, explaining the result in Fig.(7) at low power. At high heating power, the temperature gradients can pull away from marginal stability, a regime where the differences in the IFS-PPPL and gyrokinetic χ_{i}’s are less and the predicted temperature gradients differ by ∼ 20%. [The range of normalized heating power P/P_{0} of 0 to 250 in Fig.(7) will cause the gyrokinetic predicted temperature gradient to vary from R/L_{T} = 6 to 22.5, i.e., from close to marginal stability to far above marginal stability.] Thus we see that the temperature gradient predicted by the IFS-PPPL model is only 20-33% lower than the temperature gradient predicted by the gyrokinetic simulations over a wide range of heating powers. The fact that the predicted temperature gradient at fixed heating power is less sensitive to model variations than the predicted χ_{i} at fixed temperature gradient is a consequence of the critical gradient feature of ITG turbulence. Given the difficulties of the plasma turbulence problem, a turbulence theory that predicts temperature gradients to within 20-30% can be considered a significant achievement in many ways. But the fusion reaction cross-section scales as ∼ T^{2}, and the resulting fusion power feeds back to give more heating, so the performance of a fusion device at high gain (near ignition) becomes fairly sensitive to the transport. Thus one would like to have even higher accuracy in the transport model.

Next we consider the sensitivity of predictions of fusion power performance to variations in the assumed transport model. Both the 94 and 95 versions of the IFS-PPPL model are of the form

χ_{IFS−PPPL} = G (\frac{R}{L_{T}} − \frac{R}{L_{T crit}} ) W  (3)

where W and R/L_{T crit} are functions of various plasma parameters, G(x) = min(x, x^{1/2}) ×H(x), H(x) is the Heaviside function, and the argument of G is x = R/L_{T} −R/L_{T crit} as shown. In order to fit the gyrokinetic flux-tube simulations shown in Fig.(3) and described in Sec.(5), this was modified to be of the form of Eq.(1),

χ_{GK} = [\frac{R}{L_{T}} − \frac{R}{L_{T eff}} ] \frac{4.0}{R/L_{T}} W,  (4)

where R/L_{T eff} = R/L_{T crit} + 2 is assumed, and the parameterizations of R/L_{T crit} and W given by the IFS-PPPL model are used. While this fits Fig.(3) and matches Eq.(1) for the Cyclone DIII-D parameters, more detailed studies need to be carried out to develop a gyrokinetic-based model which has been tested over a wide range of parameters. For example, the nonlinear upshift in the effective R/L_{T crit} is probably not a constant value of 2 in reality and should depend on various parameters (such as collisionality, as demonstrated in recent gyrokinetic simulations[64]). But for now we will use these two equations, Eq.(3) and Eq.(4), to show the sensitivity of the predicted fusion power gain to variations in the transport model that are roughly of the magnitude represented by these two equations.

Another possible fit to the gyrokinetic flux-tube simulations in Fig.(3), which is not quite as good as Eq.(4) but is more easily compared with the IFS-PPPL model, is χ_{GK2} = 0.8G(R/L_{T} − R/L_{T eff})W. This corresponds to a rescaling and a shifting of the IFS-PPPL model to an effective critical gradient R/L_{T eff} = R/L_{T crit} + 2. The resulting predictions are fairly similar to the predictions of Eq.(4) in Fig.(8). Using R/L_{T eff} = 1.5R/L_{T crit} would lead to slightly (∼10%) lower Q predictions.

Fig.(8) shows the predicted fusion gain Q = P_{fusion}/P_{aux} for a particular reactor design, versus the boundary condition assumed for the temperature at the top the pedestal caused by the H-mode transport barrier. [In H-mode experiments the pedestal temperature can be much higher than the separatrix temperature at the last-closed flux surface, so the distinction between the two can be important.] The assumptions in Fig.(8) are described in more detail below. Predictions are shown for the standard IFS-PPPL(95) transport model, and the gyrokinetic-based version of this model (given by Eq.(4). At a fixed pedestal temperature of 3 keV, Q rises significantly, from Q=5.6 for the original IFS-PPPL model to Q=13.5 for the gyrokineticbased version. However, the results are still fairly sensitive to the assumed pedestal temperature boundary condition, so there is a risk of Q significantly below 10 even with the gyrokinetic fit. The gyrokinetic fit demonstrates primarily the sensitivity of the results to a shift in R/L_{T eff} = R/L_{T crit} + 2. There are other possible sources of uncertainty, so we also show the Q curve predicted if the gyrokinetic χ_{i} of Eq.(4) is further reduced by a factor of 0.5. In this case high Q operation, Q > 10, can be achieved even at fairly low pedestal temperatures. These results are sensitive to the achievable density. The 3 cases just described use the standard assumption in the 1996 ITER baseline scenario of a density 1.5 times the Greenwald density. Because of uncertainties about the achievable density, later ITER designs considered lower density operating points, and the effect of lowering the density to 1.15 times the Greenwald density is shown in Fig.(8) for the most optimistic of these 3 cases (χ = χ_{GK}/2).

To be clear, the 94 version of the IFS-PPPL model[8] was used in comparing with the adiabatic electron results of Fig.(3), while the 95 version of the IFS-PPPL model[9, 28] is used in predicting the fusion reactor gain Q in Fig.(8) because it includes quasilinear/mixing-length estimates of the destabilizing effects of trapped electrons that are thought to be important in real experiments. The trapped electrons cause an increase in W and a drop in R/L_{T crit}. Well above marginal stability, trapped electrons cause roughly a factor of two increase in the 95 version of χ_{i} relative to the 94 version, though this depends on parameters. There is linear and nonlinear evidence that non-adiabatic electrons can add significant additional drive to ITG turbulence and lower the ITG mode critical gradient (for example, see Refs.([21, 41, 46, 47]), and it is the 95 version of the IFS-PPPL model that has been more widely compared with experiments. As described in Sec. 5, the simulations in Fig.(3) used a simplified concentric-circle equilibrium where the linear critical gradient was about 4. The 94 and 95 IFS-PPPL models used a more realistic equilibrium and so the 94 IFS-PPPL model predicts a lower R/L_{T crit} = 3.1. The 95 IFS-PPPL model predicts that for the parameters of the DIII-D base case, trapped electrons will lower R/L_{T crit} further to 2.1, while the stabilizing influence of impurities (Z_{eff} = 2.37), beams (n_{beam}/n_{e} = 0.05), and T_{i}/T_{e} = 1.183 will raise the linear R/L_{T crit} back up to 4.2 (the most important of these 3 factors is the impurities). These various stabilizing influences will be weaker in ITER than in DIII-D, but there may be various possible sources of uncertainty in the IFS-PPPL parameterizations of these complicated effects that could be investigated further. This is beyond the scope of the present paper, but Fig.(8) can be used as a rough guide to the sensitivity of the predicted Q to changes in the critical gradient R/L_{T crit} of 2 and to changes in the coefficient of χ_{i} of a factor of 2. The GLF23 transport model[7] shares some similarities to the IFS-PPPL model and is normalized to nonlinear gyrofluid simulations, but it uses a quasilinear/mixing-length numerical procedure (similar to techniques used in the Bateman-Weiland multi-mode procedures[10]) to predict the various parametric dependences of χ_{i} instead of the analytic parameterization of the IFS-PPPL model. Predictions of ITER by the GLF23 model show similar trends as the IFS-PPPL model, including a strong dependence on the pedestal temperature[6, 7].

The calculations in Fig.(8) were done with a standard type of transport code; similar ITER-related calculations can be found in Refs.[3, 4, 5, 6, 7]. The plasma parameters used here are taken from the baseline scenario in use for the design of ITER circa 1996: major radius R = 8.14 m, midplane minor radius a = 2.8 m, elongation κ_{95} = 1.6, triangularity δ_{95} = 0.24, magnetic field B_{tor} = 5.68 Tesla, plasma current I_{p} = 21 MA, P_{aux} = 100 MW of auxiliary heating (assumed to be centrally deposited with a Gaussian half-width of r/a = 0.1, all deposited in the ions to maximize T_{i}/T_{e}). beryllium impurities with n_{Be}/n_{e} = 0.02 were assumed, and n_{He}/n_{e} was determined by τ_{p∗He}/τ_{E} = 10. The q profile was chosen so that q on axis is 0.8 (lowering the central q below 1 is favorable in this transport model), and the midplane radius of the q = 1 surface was r/a = 0.43. Sawteeth are ignored. Neoclassical ion transport is included but has little effect. Particle transport and any associated convective heat transport is ignored, and any other transport mechanisms that are sometimes thought to play a role at high β or near the edge of the plasma are neglected. Brehmsstrahlung radiation is included, but any other line radiation or charge exchange losses are ignored. The calculations in Fig.(8) also include a small amount of favorable elongation scaling found in Ref.[28], χ → χ/(1 + ((κ − 1)q/3.6)^{2}). The effects of elongation and shaping will be discussed below in more detail. Equilibrium-scale sheared flows, which are thought to get weaker in larger tokamaks, are neglected.

As in the 1996 ITER baseline scenario, a flat density profile is assumed with n_{e} = 1.3 × 10^{20}/m^{3}, corresponding to 1.5 times the Greenwald density limit. Because of uncertainties about whether ITER could operate at such a high density, later ITER designs considered lower density operating points, and the effects of lowering n_{e} to 0.98 × 10^{20}/m^{3} (corresponding to 1.15 times the Greenwald density) are shown for one case in Fig.(8). This lower density operating point includes 0.16% argon for the operation of a radiative divertor, in addition to the 2% beryllium.

The nominal design goal of ITER was P_{fusion} = 1500 MW, so in the cases in Fig.(8) where Q > 15 it would be possible to lower the auxiliary power and Q would rise further. Because of the stiffness of the temperature profiles, in cases with low Q in Fig.(8) it is sometimes possible to raise the Q by lowering the auxiliary heating power (assuming a fixed pedestal temperature), though this would not help on the power loading and nuclear testing goals which require a certain level of P_{fusion}.

To summarize the results of Fig.(8), at fixed pedestal temperature, the gyrokineticbased model achieves significantly higher Q than the original gyrofluid-based model, but the results are still sensitive to the achievable density and to the assumed pedestal temperature, which is fairly uncertain and could be less than 1.5 keV. Some experiments on the largest tokamaks find that the H-mode pedestal width scales linearly with the (poloidal) gyroradius[59, 60], while other experiments find more optimistic scalings that are weakly dependent on or independent of gyroradius[61]. Some simple theoretical models of H-mode pedestal scaling (such as Refs.[2, 62]) give pedestal widths proportional to the gyroradius ρ, consistent with some of the largest tokamak experiments, while other models give a ρ^{2/3}R^{1/3} or weaker scaling. Pedestal models with a strong dependence on gyroradius tend to predict very low pedestal temperatures when extrapolated to regimes of high density relative to the Greenwald density limit. As stated in one review[63], “While, given the present state of knowledge, we cannot provide a reliable estimate of the pedestal parameters in ITER …, a pedestal temperature less than 1500 eV, perhaps much less, is a distinct possibility.” On the other hand, there are uncertainties both ways, and there is a possibility that the pedestal temperature could be sufficiently high. Also, there are various methods that may be able to improve the temperature at or near the pedestal, such as pellet fueling, RF waves or low voltage beams to drive sheared flows, or stronger plasma shaping. Some of the new ITER-RC designs that are presently under consideration have significantly stronger plasma shaping (higher elongation and triangularity) and may have much higher pedestal temperatures than the earlier 1996 ITER design.

8. SUMMARY AND DISCUSSION

This paper has focused on the differences between the gyrofluid and gyrokinetic simulations illustrated in Fig.(3) (for the simplified case of adiabatic electrons, circular geometry, collisionless ions). The resulting sensitivity of the predicted temperature gradients are shown in Fig.(7), and the sensitivity of the predicted fusion gain in Fig.(8).

A primary result of this paper is that detailed noise and other tests have been carried out for the nonlinear gyrokinetic particle simulations, and they appear to be well-resolved and correct for this parameter regime. For the Cyclone DIII-D parameters, at the fixed temperature gradient of R/L_{T crit} = 6.9, the χ_{i} from gyrofluid flux-tube simulations are a factor of 3.3 higher than the gyrokinetic flux-tube simulations (Fig.(3)). However, turning Fig.(3) around to find the predicted temperature gradient for a fixed amount of heating power, we find that these gyrofluid errors lead to only a 20-33% drop in the predicted local temperature gradient, as shown in Fig.(7). [This relative insensitivity is a general feature of critical-gradient types of models near marginal stability.] While these errors are relatively small in one sense, because the fusion cross-section scales as T^{2}, and because the resulting fusion power feeds back to give more heating, the performance of a fusion device at high gain (near ignition) becomes fairly sensitive to the transport model and one would like to have higher accuracy than 20-33%. As shown in Fig.(8), modifying the IFS-PPPL model to better fit the gyrokinetic simulations causes the fusion gain Q to rise significantly at fixed pedestal temperature. But the results are still sensitive to the assumed pedestal temperature, which is fairly uncertain. There is still a risk of low Q, particularly if high density can not be achieved or if the pedestal temperature is low. Of course, the uncertainties go both ways, and it remains possible that the original ITER design may be adequate to achieve ignition. Other sources of uncertainty which need better understanding, in addition to the issues of the pedestal temperature and the achievable density, include the effects of elongation and plasma shaping, collisional damping of zonal flows, plasma rotation, the possibility of density peaking, and fully electromagnetic simulations with non-adiabatic electrons.

It is important to note that the present gyrokinetic/gyrofluid comparisons have been done in a simplified case of collisionless electrostatic ITG turbulence with adiabatic electrons. Additional physics has the potential to reduce the differences. One of the main causes of the gyrofluid/gyrokinetic differences is that in the nonlinear gyrokinetic simulations there is a significant nonlinear upshift in the effective critical gradient, R/L_{T eff} ≈ R/L_{T crit} + 2, due to the generation of an undamped component of zonal flow[43], which then suppresses further turbulence. This nonlinear upshift is probably not a constant value of 2 and future work should investigate how it depends on various parameters. Linear gyrokinetic calculations by Dorland have found that this undamped component of the zonal flows may increase with elongation (and thus make the gyrofluid/gyrokinetic disagreement worse, though the other effects described next will tend to reduce the disagreement). Recent work by Z. Lin et.al.[64] shows that including a small but experimentally-relevant amount of collisions can lead to enough damping of the zonal flows to cause a significant increase in the heat conductivity χ_{i}, particularly near marginal stability. This will tend to reduce the effective critical gradient back towards its linear value and will cause the gyrokinetic χ_{i} to rise to be closer to the gyrofluid-based IFS-PPPL χ_{i}, although the resulting χ_{i} near marginal stability will probably depend strongly on ion-ion collisionality (an effect that is not in the IFS-PPPL model or any other current transport model). Recent work by Diamond and others[65, 64]indicates that the zonal flow generation observed near marginal stability is related to inverse-cascades and a resulting negative turbulent viscosity. The collisionless nonlinear gyrokinetic simulations presented here demonstrate that far above marginal stability it is possible to reach a steady state where the nonlinear generation of zonal flows is balanced by the turbulent viscosity, so a non-zero χ_{i} can be achieved that presumably does not depend strongly on ion-ion collisions. An interesting topic for future work would map out the transition between these two regimes of strong or weak collisionality dependence. Non-adiabatic electrons are known to limit inverse cascades in some types of plasma turbulence, and so including non-adiabatic electrons in gyrokinetic ITG simulations may further reduce the gyrofluid/gyrokinetic differences. Non-adiabatic electrons may also push the plasma into stronger-turbulence regimes where the differences are less important.

There has been some recent work[44] on developing improved gyrofluid closures that incorporate some of the neoclassical effects needed to model an undamped component of the zonal flows. As shown in Fig.(3), this new gyrofluid closure eliminates about half of the original gyrofluid/gyrokinetic difference, and helps isolate the source of the differences. Possible further improvements in the closures will be investigated.

The IFS-PPPL and gyrofluid results in Fig.(3) are stiffer than the Weiland and gyrokinetic results, that is, they show a steeper increase in transport going above threshold. Global gyrokinetics show yet lower transport levels, but are still limited to values of a/ρ lower than the regime of ITER. It would be interesting to undertake an a/ρ scan about the DIII-D base parameter set. Such a scan is probably now possible, even if not to ITER-like values. The Weiland model gives transport levels in the same range as the flux-tube gyrokinetic results for these parameters (though like the IFS-PPPL model it misses the nonlinear upshift in the critical gradient). It is interesting that the Weiland model agrees better with gyrokinetic simulations in this circular limit, since the Weiland model is based on a simpler fluid theory than the gyrofluid model.

The elongation scaling is an important issue for determining the implications for ITER or other reactor designs. The nonlinear comparisons so far have been in simple circular geometry. There are significant differences in the transport models regarding elongation when scaled to ITER plasmas. The Multi-Mode model has an empirical elongation scaling based on experimental data for the scaling of the confinement time with the current, which leads to χ_{i} scaling asymptotically as κ^{−4}, where κ is the elongation parameter. The IFS-PPPL model has a much weaker elongation dependence, based on initial linear gyrokinetic and nonlinear gyrofluid calculations[28] which saw little effect at moderate elongation (κ < 1.6) typical of the original ITER design. However, with very strong plasma shaping (very high elongation and triangularity, or at high Shafranov shift and/or low aspect ratio) there may be a transition to a regime of significantly improved confinement[46, 66]. The shaping criterion needed for this improvement needs to be studied in more detail. Given the constraints on the models described, the predictions for ITG transport in an ITER plasma based on extrapolations from the gyrokinetic simulations and the Weiland model are less pessimistic than the IFS-PPPL model predictions based on gyrofluid simulations. ITER predictions for the Weiland-based Multi-Mode model are in Ref.[3]. It is interesting to note that there are a number of different transport models with fairly different scalings on various parameters (and different amounts of stiffness) that are nevertheless able to achieve similar levels of comparisons with the experimental profile database that has been developed for ITER[4]. For example, the IFS-PPPL model is able to follow many of the main trends in this database from circular to elongated tokamaks despite its weak elongation dependence because of other factors which are also correlated with elongation (higher edge temperatures, broader q profiles, and toroidal rotation). More detailed studies, such as with controlled rotation scans or with perturbative heat pulse propagation experiments, should help to better distinguish between transport models. More complete three-dimensional nonlinear simulations, building on the simulations done here, will also help in developing transport models that can be used to predict and optimize the design of fusion reactors.

Possible future tasks that may shed more light on the differences between various transport models, and help develop more accurate transport models, include more detailed comparisons between codes and experiments of fluctuation spectra, and of the poloidal flow dynamics. More complete scans of a/ρ with the global codes with the same toroidal resolution as for the flux-tube codes will be of great interest as such scans become possible. Further comparisons which systematically remove the simplifications made in the comparisons made here, including shaped equilibria, equilibrium-scale sheared flows, nonadiabatic electrons, and beam and impurity species are also essential to making quantitative assessments of differences in predictions of ITER performance between the models. Work is also under way to design experiments that the turbulence simulations can model with more certainty.

ACKNOWLEDGMENTS

Several of us have had useful discussions on these topics with C. Bolton, Z.-H. Lin, F. Hinton, M.N. Rosenbluth, and R.E. Waltz. The ITG simulations were performed on the C-90 and T3E computers at NERSC. This work was performed under the auspices of U.S. Department of Energy (DOE) by Lawrence Livermore National Laboratory under contract No. W-7405-ENG-48, by Princeton Plasma Physics Laboratory under contract No. DE-AC02-76CH03073, and by various other DOE contracts and grants, and is part of the Numerical Tokamak Turbulence Project supported by the High Performance Computing and Communications Program in the U.S. DOE.

References

[1] “ITER Project: A Physics and Technology Experiment”, R. Aymar, V. Chuyanov, M. Huguet, R. Parker and Y. Shimomura, Proceedings of the 1996 IAEA meeting, volume 1, page 3.

[2] M. Kotschenreuther, W. Dorland, Q. P. Liu, G. W. Hammett, M. A. Beer, S. A. Smith, A. Bondeson, and S. C. Cowley, in Proc. 16th. Int. Conf. on Fusion Energy, 1996 (IAEA, Vienna, 1997), Vol. 2, p. 371; M.

Kotschenreuther and W. Dorland, “Memorandum on Confinement Projections To: FESAC ITER Confinement reviewers”, (Feb. 14, 1997)

http://w3.pppl.gov/∼hammett/ark/1997/FESAC ITER Review memo.pdf
http://w3.pppl.gov/∼hammett/ark/1997/FESAC ITER Review memo.doc.

(GWH ??: will submit these to the AIP or NTIS online archives so they can be formally cited in Phys. Plasmas.)

[3] G. Bateman, A. H. Kritz, J. E. Kinsey, and A. J. Redd, Phys. Plasmas 5, 2355 (1998).

[4] D.R. Mikkelsen, G. Bateman, D. Boucher, et al. Proc. 17th Int. Conf. on Fusion Energy, (Yokohama, 1998) (IAEA, Vienna, 1999); to be published in Nucl. Fusion 39, (1999).

[5] J.W. Connor, M. Alexander, S.E. Attenberger, et al., in Proc. 16th. Int. Conf. on Fusion Energy, 1996, Vol. 2, p. 935 (IAEA, Vienna, 1997).

[6] J.E. Kinsey, R.E. Waltz, D.P. Schissel, Proc. of the 24th Eur. Phys. Soc. Conf. on Controlled Fusion and Plasma Physics, Contributed Papers, Part III, p. 1081 (1997).

[7] R.E. Waltz, G.M. Staebler, W. Dorland, G.W. Hammett, M. Kotschenreuther, J.A. Konings, Phys. Plasmas 4, 2482 (1997).

[8] W. Dorland, M. Kotschenreuther, M. A. Beer, G. Hammett, et al, “Comparisons of Nonlinear Toroidal Turbulence Simulations with Experiment,” in Proc. 15th. Int. Conf. on Plasma Physics and Controlled Nuclear Fusion, 1994 (1995, IAEA, Vienna), Vol. 3, 463.

[9] M. Kotschenreuther, W. Dorland, M.A. Beer, and G.W. Hammett, Phys. Plasmas 2, 2381 (1995).

[10] G. Bateman, J. Weiland, H. Nordman, J. Kinsey, C. Singer, Physica Scripta 51 591 (1995).

[11] M. A. Beer and G. W. Hammett, Phys. Plasmas 3, 4046 (1996).

[12] A. M. Dimits, T. J. Williams, J. A. Byers, and B. I. Cohen, Phys. Rev. Lett. 77, 71 (1996).

[13] R. D. Sydora, Physica Scripta 52, 474 (1995).

[14] M. Kotschenreuther, G. Rewoldt, and W.M. Tang, Comp. Phys. Comm. 88, 128 (1995), and references therein.

[15] E. A. Frieman and Liu Chen, Phys. Fluids 25, 502 (1982).

[16] C.M. Greenfield, J.C. DeBoo, T.H. Osborne, F.W. Perkins , M.N. Rosenbluth, D. Boucher, Nuclear Fusion 37, 1215 (1997).

[17] B. I. Cohen, D. C Barnes, J. M. Dawson,et al., Comp. Phys. Commun. 87, 1 (1995).

[18] S.E. Parker, W. Dorland, R.A. Santoro, M.A. Beer, Q.P.Liu, W.W. Lee, and G.W. Hammett, Phys. Plasmas 1, 1461 (1994).

[19] R.E. Waltz, G.D. Kerbel, and J. Milovich, Phys. Plasmas 1, 2229 (1994).

[20] M. A. Beer, S. C. Cowley and G. W. Hammett, Phys. Plasmas 2, 2687, 1995.

[21] M.A. Beer, Ph.D. thesis, Princeton University, 1995.

[22] J. Weiland, A. Jarmen, and H. Nordman, Nucl. Fusion 29, 1810 (1989).

[23] J. Weiland, A. Hirose, Nucl. Fusion 32, 151 (1992).

[24] P. N. Guzdar, J. F. Drake, D. McCarthy, A. B. Hassam, and C. S. Liu, Phys. Fluids B5, 3712 (1993).

[25] G. Bateman, A. Kritz, J. Kinsey, A. Redd and J. Weiland, Phys. Plasmas, 5, 1793 (1998).

[26] J. Kinsey, G. Bateman, A. Kritz and A. Redd, Phys. Plasmas, 3, 561 (1996).

[27] S.D. Scott, G.W. Hammett, C.K. Phillips, et al., Proc. 16th Int. Conf. on Fusion Energy, (Montreal, 1996) (IAEA, Vienna, 1997), Vol. 1, p. 573.

[28] W. Dorland, M. Kotschenreuther, Q.P. Liu, M.A. Beer, G.W. Hammett, in Proc. of the Joint Varenna-Lausanne Int. Workshop on Theory of Fusion Plasmas (August 1996), Varenna, Italy 1996 (Editrice Compositori, Bologna, Italy), p. 185.

[29] X. Garbet and R.E. Waltz, Phys. Plasmas 3, 1898 (1996); Phys. Plasmas 5 2836 (1998).

[30] W. W. Lee, J. Comput. Phys. 72, 243 (1987).

[31] A. M. Dimits, Ph.D. thesis, Princeton University, 1988.

[32] M. Kotschenreuther et al., in Proc. 14th. Int. Conf. on Plasma Physics and Controlled Nuclear Fusion, 1992 (1993, IAEA, Vienna), Vol. 2, 11.

[33] S. E. Parker, W. W Lee and R. A. Santoro, Physical Review Letters 71, 2042 (1993).

[34] B. I. Cohen T. J. Williams, A. M. Dimits, and J. A. Byers, Phys. Fluids B5, 2967 (1993).

[35] W. D. Dorland, Ph.D. thesis, Princeton University, 1993.

[36] G. W. Hammett, M. A. Beer, W. Dorland, S. C. Cowley, and S. A. Smith, Plasma Phys. Control. Fusion, 35, 973 (1993).

[37] A. M. Dimits, J. A. Byers, T. J. Williams, B. I. Cohen, et al., Fifteenth International Conference on Plasma Physics and Controlled Nuclear Fusion Research, Sept. 26, 1994, Seville, Spain, IAEA-CN-60/D-P-I-5 (IAEA, Vienna, 1995).

[38] A. M. Dimits, Phys. Rev. E48, 4070 (1993).

[39] Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang, and R. B. White, Science 281, 1835(1998).

[40] T. S. Hahm, W. W. Lee, and A. Brizard, Phys. Fluids 31, 1940 (1988).

[41] R.D. Sydora, V.K. Decyk and J.M. Dawson Plasma Phys. Control. Fusion 38, A281 (1996).

[42] S.E. Parker, C. Kim, Y. Chen, Physics of Plasmas 6, 1709 (1999).

[43] M. N. Rosenbluth and F. L. Hinton, Phys. Rev. Lett. 80, 724 (1998).

[44] M. A. Beer and G. W. Hammett, in Proc. of the Joint Varenna-Lausanne Int. Workshop on Theory of Fusion Plasmas (August 1998), Varenna, Italy 1998 (Editrice Compositori, Bologna, Italy, 1999).

[45] A. M. Dimits, B. I. Cohen, N. Mattor, W. M. Nevins, and D. E. Shumaker, S. E Parker and C. Kim, Proceedings of the 17th. IAEA Fusion Energy Conference (paper IAEA-F1-CN-69/TH1/1, accepted for publicaton in Nuclear Fusion).

[46] M. A. Beer, G. W. Hammett, G. Rewoldt, E. J. Synakowski, M. C. Zarnstorff, and W. Dorland, Phys. Plasmas 4, 1792 (1997).

[47] A. J. Redd, A. H. Kritz, G. Bateman, G. Rewoldt, and W. M. Tang, Phys. Plasmas 6, 1162 (1999).

[48] R.E. Waltz, R.R. Dominguez, and G.W. Hammett, Phys. Fluids B4, 3138 (1992).

[49] M.A. Beer, G.W. Hammett, Phys. Plasmas 3, 4046 (1996).

[50] Z. Chang and J.D. Callen, Phys. Fluids B 5, 1167 (1992).

[51] N. Mattor and S.E. Parker, Phys. Rev. Lett. 79, 3422 (1997)

[52] C. H. Kim and S. E. Parker, Bull. Amer. Phys. Soc. 43, 1721 (1998).

[53] T.S. Hahm and K.H. Burrell, Phys. Plasmas 2, 1648 (1995).

[54] K.H. Burrell, Phys. Plasmas 4, 1499 (1997).

[55] E.J. Synakowski, Plasma Phys. Control. Fusion 40, 581 (1998).

[56] D. R. Ernst, S. Batha, M. Beer, et al., Proc. 17th Int. Conf. on Fusion Energy, (Yokohama, 1998) (IAEA, Vienna, 1999); to be published in Nucl. Fusion 39, (1999).

[57] T. S. Hahm, M.A. Beer, Z. Lin, G.W. Hammett, W.W. Lee, and W.M. Tang, Phys. Plasmas 6, 922 (1999).

[58] The scrambling noise test was proposed originally by M. Kotschenreuther, and valuable variations and refinements were suggested by W. Dorland, G.W. Hammett and M.N. Rosenbluth in private communications dating between 1995 and 1997.

[59] T. Hatae, Y. Kamada, S. Ishida, et al., Plasma Phys. Control. Fusion 40, 1073 (1998).

[60] J.G. Cordey and the JET Team, Proc. 17th Int. Conf. on Fusion Energy, (Yokohama, 1998) (IAEA, Vienna, 1999).

[61] T.H. Obsborne, K. H Burrell, R.J. Groebner, et al., to be published in Plasma Phys. and Control. Fusion.

[62] F.W. Perkins, D. Post, M.N. Rosenbluth, et al., in Proc. 16th. Int. Conf. on Fusion Energy, 1996 (IAEA, Vienna, 1997), Vol. 2, p. 963.

[63] Panel Report To Fusion Energy Sciences Advisory Committee (FESAC) - “Review of the International Thermonuclear Experimental Reactor (ITER) Detailed Design Report”, April 18, 1997, Appendix D.I.

http://wwwofe.er.doe.gov/More HTML/REPORTAP (for Microsoft Word),
http://w3.pppl.gov/∼hammett/ark/1997/REPORTAP.pdf
http://wwwofe.er.doe.gov/More HTML/FESACREP.pdf

[64] Z. Lin, personal communication.

[65] P. Diamond, personal communication.

[66] M. Kotschenreuther, W. Dorland, Q.P. Liu, M. Zarnstorff, R.L. Miller, Y.R. Lin-Liu, Proc. 17th Int. Conf. on Fusion Energy, (Yokohama, 1998) (IAEA, Vienna, 1999); to be published in Nucl. Fusion 39, (1999).

Figure Captions

Fig. (1) γ and ω_{r} vs. k_{θ} for the Sydora (global) and Dimits (flux-tube) nonlinear gyrokinetic codes, for the Kotschenreuther and Rewoldt (Full) linear gyrokinetic codes, the Beer nonlinear gyrofluid code, and the Weiland fluid calculation.

Fig. (2) Residual E×B flow fraction vs. Rosenbluth-Hinton parameter h ≡ \sqrt{ϵ}/q^{2}. The points are results from the Dimits et al. flux-tube gyrokinetic code and the line is the prediction of Ref. [43].

Fig. (3) χ_{i} vs. R/L_{T} from the gyrofluid code using the 1994 “thesis closure” [21], an improved 1998 gyrofluid closure, the 1994 IFS-PPPL model [8], the LLNL and U. Colorado flux-tube and UCLA (Sydora) global gyrokinetic codes, and the MMM model for the DIII-D base case.

Fig. (4) The dependence of χ_{i} on ϵ = r/R from gyrofluid and gyrokinetic simulations.

Fig. (5) Normalized χ_{i} vs. tv_{ti}/L_{T} from gyrokinetic simulations with particle numbers ranging from 5 × 10^{5} to 1.34 × 10^{8}, corresponding to 1 to 128 particles per cell..

Fig. (6) Time histories of (a) χ_{i} and (b) the mean squared E×B shearing rate (L_{T i}S/v_{ti})^{2} associated with the flux-surface-averaged electrostatic potential, both in the absence of scrambling and when restarts with scrambling and gradient reduction are done at three times during the run. Cyclone DIII-D base-case parameters were used.

Fig. (7) The ratio of the temperature gradient predicted by the IFS-PPPL model to the temperature gradient predicted by the gyrokinetic flux-tube simulations, versus the normalized heating power, based on the results in Fig.(3). Thus the gyrofluid/gyrokinetic differences in Fig.(3) correspond to only a 20-33% difference in the predicted temperature gradient at a fixed amount of heating power.

Fig. (8) The predicted fusion gain Q vs. assumed pedestal temperature, for the IFS-PPPL 95 model, for a modified model to fit the gyrokinetic flux-tube results of Fig.3 (“GK fit”), and for a further reduction in χ_{i} by a factor of 2 (“GK fit/2”). These 3 cases are at 1.5 times the Greenwald density. Also shown is a lower density case at 1.15 times the Greenwald density using the “GK fit/2” χ_{i}.

Fig. (1) γ and ω_{r} vs. k_{θ} for the Sydora (global) and Dimits (flux-tube) nonlinear gyrokinetic codes, for the Kotschenreuther and Rewoldt (Full) linear gyrokinetic codes, the Beer nonlinear gyrofluid code, and the Weiland fluid calculation.

Fig. (2) Residual E×B flow fraction vs. Rosenbluth-Hinton parameter h ≡ \sqrt{ϵ}/q^{2}. The points are results from the Dimits et al. flux-tube gyrokinetic code and the line is the prediction of Ref. [43].

Fig. (3) χ_{i} vs. R/L_{T} from the gyrofluid code using the 1994 “thesis closure” [21], an improved 1998 gyrofluid closure, the 1994 IFS-PPPL model [8], the LLNL and U. Colorado flux-tube and UCLA (Sydora) global gyrokinetic codes, and the MMM model for the DIII-D base case.

Fig. (4) The dependence of χ_{i} on ϵ = r/R from gyrofluid and gyrokinetic simulations. 30 Fig. (5) Normalized χ_{i} vs. tv_{ti}/L_{T} from gyrokinetic simulations with particle numbers ranging from 5 × 10^{5} to 1.34 × 10^{8}, corresponding to 1 to 128 particles per cell..

Fig. (6) Time histories of (a) χ_{i} and (b) the mean squared E×B shearing rate (L_{T i}S/v_{ti})^{2} associated with the flux-surface-averaged electrostatic potential, both in the absence of scrambling and when restarts with scrambling and gradient reduction are done at three times during the run. Cyclone DIII-D base-case parameters were used.

Fig. (7) The ratio of the temperature gradient predicted by the IFS-PPPL model to the temperature gradient predicted by the gyrokinetic flux-tube simulations, versus the normalized heating power, based on the results in Fig.(3). Thus the gyrofluid/gyrokinetic differences in Fig.(3) correspond to only a 20-33% difference in the predicted temperature gradient at a fixed amount of heating power.

Fig. (8) The predicted fusion gain Q vs. assumed pedestal temperature, for the IFS-PPPL 95 model, for a modified model to fit the gyrokinetic flux-tube results of Fig.3 (“GK fit”), and for a further reduction in χ_{i} by a factor of 2 (“GK fit/2”). These 3 cases are at 1.5 times the Greenwald density. Also shown is a lower density case at 1.15 times the Greenwald density using the “GK fit/2” χ_{i}.