Validation studies of gyrokinetic ITG and TEM turbulence simulations in a JT-60U tokamak using multiple flux matching

Quantitative validation studies of flux-tube gyrokinetic Vlasov simulations on ion and electron heat transport are carried out for the JT-60U tokamak experiment. The ion temperature gradient (ITG) and/or trapped electron modes (TEM) driven turbulent transport and zonal flow generations are investigated for an L-mode plasma in the local turbulence limit with a sufficiently small normalized ion thermal gyroradius and weak mean radial electric fields. Nonlinear turbulence simulations by the GKV code successfully reproduce radial profiles of the ion and electron energy fluxes in the core region. The numerical results show that the TEM-driven zonal flow generation in the outer region is more significant than that in the core region with ITG- and ITG–TEM-dominated turbulence, leading to moderate transport shortfall of the ion energy flux. Error levels in the prediction of the ion and electron temperature gradient profiles in the core region are estimated as less than ±30%, based on a multiple flux matching technique, where the simulated ion and electron energy fluxes are simultaneously matched to the experimental values.


Introduction
Performance of particle and energy confinement in magnetic fusion plasmas is mainly influenced by turbulent transport driven by electrostatic and electromagnetic microinstabilities. Physical understanding and quantitative prediction of the turbulent transport are central issues in fusion plasma researches. The five-dimensional nonlinear turbulence simulation based on the gyrokinetic theory (see, e.g. [1,2] and references therein) is widely recognized as a promising way to address the turbulent transport issues, and one can find comprehensive reviews on the gyrokinetic simulations in [3,4]. Great efforts have been devoted so far to the development of the radially local simulation approaches and their extension to a global model, and then the ion temperature gradient (ITG) and the trapped electron mode (TEM) driven turbulent transport has been extensively investigated. Systematic crosscode benchmark tests have confirmed that the ITG driven turbulent heat transport level in the global simulation well converges to the local flux-tube one in a limit of sufficiently small ion thermal gyroradius (ρ ti ) compared to the plasma size (a) and/or the characteristic profile width of the logarithmic ion temperature gradient (∆ ∇ T ln i ), i.e. / / ρ ρ = < * a 1 300 ti or / / ρ ρ = ∆ < * ∇ 1 300 T eff ti ln i [7,8]. The local gyrokinetic simulations are, therefore, well applicable to the next-generation of large tokamak devices, such as ITER and DEMO with

Nuclear Fusion
Validation studies of gyrokinetic ITG and TEM turbulence simulations in a JT-60U tokamak using multiple flux matching / ρ < * 1 500, and the evaluation of their prediction capability through quantitative comparisons with the existing plasma experiments is a crucial requirement.
In the last decade, validation studies of gyrokinetic simulation codes have been actively carried out for various tokamaks (e.g. DIII-D [9][10][11][12], ASDEX Upgrade [13], Alcator C-MOD [14], etc) and helical (e.g. LHD [15,16]) devices. In the tokamak validation studies, a significant underprediction of the heat flux, which is the so-called 'transport shortfall' [10], has been identified in the outer region around / ⩾ ρ = r a 0.7, and the numerical and physical effects on the transport shortfall are investigated for the DIII-D L-mode plasma as a representative case [9,11,17]. Particularly, Görler et al has recently pointed out that the significant ion transport shortfall is resolved within the measurement errors of the ion temperature gradient, based on the flux matching technique [18], taking into account the stiff dependence of the turbulent heat flux on the ion temperature gradient [11]. Also, the strong impact of the mean radial electric field shearing on the simulated ion heat flux has been confirmed.
On the other hand, characteristics of the ITG and/or TEM driven zonal flows, which depend on the radial position in the L-mode plasma, and the relation with the transport shortfall, have not been investigated in earlier works. Indeed, since the transport suppression by zonal flows is expected to be more significant in future large tokamaks such as ITER than the mean radial electric field shear, which is roughly proportional to ρ * , gyrokinetic-simulation-based studies of the zonal flow properties are indispensable for improving the prediction capability. As a fundamental study on the nonlinear interactions of zonal flows and turbulence, the gyrokinetic entropy transfer analyses have been carried out for the ITG, TEM, and electron temperature gradient (ETG) turbulence [19][20][21][22].
In this paper, the first quantitative validation study on the ion and electron heat transport, including the zonal flow analysis for a JT-60U tokamak plasma, is presented. Here, the JT-60U L-mode discharge is carefully chosen such that the local limit condition of / ρ < * 1 300 i is well satisfied, i.e. the mean radial electric field shearing rate is negligibly small in comparison to the linear-mode growth rate. The nonlinear ITG and/or TEM turbulence simulations under the experimental conditions are performed by using an electromagnetic gyrokinetic Vlasov code GKV [23,24]. Then, the experimental equilibrium profiles, power and particle balances are provided by an integrated tokamak transport solver TOPICS [25,26]. Furthermore, the conventional ion heat flux matched by adjusting a single parameter of the ion temperature gradient [11] is extended to a multiple flux matching associated with the nonlinear dependence on the ion and electron temperature gradients, where the simulated ion and electron energy fluxes are simultaneously matched to the experimental values. Then, we evaluate the accuracy of the temperature-gradient profile prediction in the core region.
The rest of the paper is organized as follows. The calculation model in GKV under the experimental condition produced by TOPICS is presented in section 2. Then, plasma parameters, equilibrium profiles, and the linear instabilities of the present L-mode plasma are shown in section 3. In section 4, nonlinear simulation results on the ITG and/or TEM driven turbulent transport and zonal flows at several radial positions are presented. The quantitative comparisons with the experimental results and the multiple flux matching are also discussed. Finally, a summary of this paper is given in section 5.

Flux-tube gyrokinetic simulation model for realistic tokamak equilibria
In this section, we briefly summarize the gyrokinetic simulation model used in the electromagnetic gyrokinetic Vlasov code GKV and an interface with the integrated tokamak transport solver TOPICS providing the experimental equilibrium profiles. The detailed descriptions are also given in [24]. One of the governing equations is the electromagnetic gyrokinetic equation describing the time evolution of the perturbed gyrocenter distribution function ( ) δ f s g on the five-dimensional , where the subscript 's' is the index of the particle species. The Fourier representation with the perpend icular wavenumber ⊥ k is given by with the perpendicular velocity ⊥ v . The gyro-averaged potential fluctuation is denoted by 0s 0 s is the zeroth-order Bessel function, and the former and latter terms mean the electrostatic and electromagnetic parts, respectively. Since we focus here on the finite but low-β L-mode plasmas, the parallel magnetic field fluctuation ∥ δ ⊥ B k (so called the compressional component) is ignored. The equilibrium part of the distribution function is given by the local Maxwellian distribution, i.e.
( / ) Ms s s s 3 2 s 2 s , where n s represents the equilibrium density. The symbol ∑ ∆ appearing in the nonlinear term of equation (1) means the double summations with respect to k ′ ⊥ and k″ ⊥ , which satisfy the triadinteraction condition of k k k″ = + ′ ⊥ ⊥ ⊥ . Collisional effects are introduced in terms of a linearized model collision operator C s , where a gyro-averaged Lenard-Bernstein model [27] is applied here.
The gyrokinetic equation shown in equation (1) is solved in the local flux-tube coordinates (x, y, z) [28] defined as , , , where a and ( ) ρ q 0 denote the plasma minor radius and the safety factor on the flux surface label of interest ρ 0 , respectively. In these coordinates, the magnetic field vector is given by = ∇ × ∇ B B x y ax with the field strength on the magnetic axis B ax . The perpendicular wavenumber vector and the parallel derivative operator are given as ij 1 2 is the Jacobian given by the contravariant components of the metric tensor g ij for (i, j) = (x, y, z). The magnetic and diamagnetic drift frequencies, ω Ds and ω * T s , are then given by . The geometric coefficients K x and K y are defined as follows: ln , x xz xy xx yz 2 ax 2 (4) ln .
Note that the low-β approximation is applied to ω Ds so that the finite-β effect is ignored in the curvature drift. The electromagnetic potential fluctuations are determined by the Poisson and Ampère equations: where λ π = ∑ noted that, in the adiabatic electron limit with ρ ⊥ k 1 te , the gyrocenter density fluctuation for electrons is approximated by , and the field-line average is defined by . By using the gyrokinetic and the Poisson-Ampère equations, one can derive another important equation describing the balance and transfer of the entropy variable where the superscripts '(trb)' and '(zf )' mean the non-zonal and zonal components in the wavenumber space, respec- . Each term represents (i) the variation of the entropy variable, (ii) the variation of the field energy, (iii) the nonlinear entropy transfer from non-zonal to zonal modes, (iv) the entropy production by turbulent particle and heat fluxes, and (v) the collisional dissipation, respectively (see, e.g. [29] for their definitions). Note that, by definition, the turbulent-flux driven entropy production term does not appear for the zonal modes in equation (9). The entropy balance/transfer equation provides us with a good measure for the turbulence simulation accuracy as well as useful physical insights associated with the turbulence saturation mechanisms. The detailed numerical analyses of the entropy balance and the transfer processes in ITG, TEM, and ETG turbulence are shown in, e.g. [20][21][22][23].

Interface with integrated-transport solver
In this section, we present the framework of combined analyses by means of GKV and the integrated-transport solver TOPICS with experimental density and temperature profiles and magnetic fields. In GKV, effects of magnetic field geometries are incorporated into the field intensity B, the parallel derivative ⋅ ∇ b , the magnetic drift ω Ds , the perpendicular wavenumber ⊥ k , and the field-line-average operator ⟨ ⟩ z shown in section 2.1, through the contravariant components of metric tensor g ij and the Jacobian g xyz . Figure 1 shows a schematic procedure to make the realistic MHD equilibrium data for turbulence simulations with GKV. First, the MHD equilibrium based on the experimentally measured density, temperature, and current profiles is reconstructed by TOPICS. The radial profiles of the ion and electron energy fluxes are simultaneously provided by the power balance analysis for, e.g. Ohmic, NBI, and EC heatings, where the power depositions and losses are calculated by a fast-ion orbit following Monte-Carlo Fokker-Planck solver OFMC [25,30]. After that, a flux coordinate generator IGS [24] makes high-accuracy interpolations of the flux surfaces and constructs the straight-field-line coordinates, such as the axisymmetric (or natural), Boozer, and Hamada coordinates. Using these coordinates data, the metric tensor and the Jacobian for geometry-dependent quantities and operators are calculated in GKV. It is noted that the up-down asymmetry resulting from the divertor configuration is consistently included in this framework.
The combined analyses among GKV and TOPICS via IGS enable us to make not only experimental analyses, but also predictive studies for future devices, e.g. shape optim ization on the microinstability and turbulent transport. Actually, in [24], GKV with realistic tokamak equilibria is numerically verified through the cross-code benchmark test. It is then applied to two types of shaped plasmas expected in the JT-60SA tokamak device, i.e. ITER-like and highly-shaped plasmas, where a decrease in the ITG-mode growth rate and enhancement of the residual zonal flow level in the highlyshaped case are identified.

JT-60U L-mode equilibrium and linear instability analyses
Using the framework of GKV-TOPICS, gyrokinetic simulation studies for the JT-60U tokamak experiments are carried out. Then, the prediction capability of GKV is examined through comparisons between the simulation and experimental results on the turbulent ion and electron energy fluxes. For a validation study of the local flux-tube gyrokinetic simulation, we have chosen a JT-60U deuterium plasma on the L-mode discharge with a positive magnetic shear profile [31], where  T eff ti i for ρ < 0.8 well satisfy the local limit condition of / ρ < * 1 300 [5][6][7][8]. The configuration parameters of the experiment are summarized in table 1.
The plasma shape and equilibrium profiles are shown in figure 2. It is noted that the electron temperature is slightly higher than the ion one for ⩽ ρ 0.7. The mean toroidal rotation and the mean poloidal × E B rotation are negligibly small in the whole radial region, i.e. | |∼ One also finds almost monotonically increasing profiles of the normalized density/temperature gradient, with the characteristic collision time τ ′ ss . Besides, the mean radial electric field shearing rate γ Er is about ten times smaller than the linear-mode growth rate γ lin . The local simulation approach is, thus, well justified for the plasma investigated here.
The wavenumber dependence of the linear instability growth rate γ lin (k y ) calculated by the linear GKV simulations is shown in figure 3, where γ lin (k y ) for three radial positions, ρ = 0.26, 0.50, 0.76 are plotted. Since, as shown in figure 2(g), the normalized gradients rapidly increase towards the plasma edge, the different linear microinstabilities appear at each radial position, where the ITG, ITG-TEM, and TEM instabilities are dominant in the ion-scale modes of ⩽ ρ k 2 y ti at ρ = 0.26, 0.50, 0.76, respectively. The ETG modes are also unstable at ρ = 0.50 and 0.76. The ITG-and ETG-mode growth rates calculated with adiabatic electrons or ions are also plotted. One can see that the adiabatic-electron approximation is only valid in the deep core region (ρ ∼ 0.26) where a fraction of the trapped electrons becomes less significant in comparison to that in the outer region. It is also noteworthy that the maximum growth rate in the ion-scale modes with ρ < k 1 y ti increases towards the outer region.

Numerical settings and the entropy balance relations
Turbulence simulation results and comparisons with the experimental measurements are discussed in this section. The ITG and/or TEM turbulence simulations at several points in     (8) and (9), respectively.

Comparisons of turbulence simulation results with exper imental measurements
The simulation results on the turbulent energy fluxes for ions and electrons are compared with experimental measurements. Figure 5 show the spectro-temporal evolutions of the radial turbulent energy flux density  . Indeed, the wavenumber at the spectral peak is similar to or lower than that giving the maximum linear growth rate (shown by the horizontal dashed line in the contours). It is also found that contribution of the higher-wavenumber modes is more significant for the electron energy flux, where the amplitude | | ⊥ Q k e is about ten times larger than that | | The radial dependencies of the turbulence part W trb with ≠ k 0 y and the zonal flow part W zf with k y = 0 in the generalized flow energy W total are shown in figures 6(a) and (b), where W total is defined by representing the nonlinear source for the zonal flow generation are also plotted in the figures. We see that the monotonic increase in the turbulence energy W trb is well correlated with the mixing-length estimate. As for the zonal flow energy normalized by the turbulence part / W W zf trb , a weak radial dependence is found for the core region of ⩽ ρ 0.50 indicating ITG-and ITG-TEM-dominated turbulence, while more significant zonal flow generation is observed for the outer region with TEM-dominated turbulence. This is consistent with the tendency of the entropy transfer rate in figure 6(b). The understanding of the correlation among the turbulence energy, zonal flow energy, and the entropy transfer rate is useful for constructing a simplified turbulent transport model [15] which can be applied to the integrated transport simulations [32]. Figure 7 show the comparisons between the GKV simulation results and the experimental measurements on the radial profiles of the ion energy flux P i , the electron energy flux P e , and the convective part resulting from the turbulent particle flux ( / ) Γ T 5 2 e e in the unit of MW. The fluxsurface-integrated radial energy and particle fluxes are defined as   S eq i e ), radiation loss (S rad ), and charge exchange loss (S CX ). Note that the radiation loss is not taken into account in this analysis. For comparison, the results with the adiabatic-electron ('ad.-elec.') approximation and with a gyrofluid-based transport model TGLF [33] are also plotted. Note that, since the L-mode plasma chosen here shows negligibly small radial electric fields (see figure 2(e)), weak impact of the radial × E B shearing effect is confirmed in the TGLF cases. One finds that the GKV simulations with kinetic electrons reproduce transport levels relevant to the experimental measurements in the core region ( ⩽ ρ 0.50) for both the ion and electron energy fluxes, whereas the adiabatic-electron case predicts lower transport levels. Also, TGLF and GKV results show similar energy transport levels for the ITG-dominated region of ⩽ ρ 0.4, but some differences are found in the outer region with the TEM-dominated turbulence. As for the particle flux Γ e , a slightly lower transport level is found in the GKV results, but still the same order of magnitude in comparison to that in the TGLF cases. Note that the experimental evaluation of Γ e is not available due to large uncertainty in the particle source and sink depositions, which are strongly influenced by the particle recycling and the ionization/recombination processes in the plasma edge.
It is stressed that the ion energy flux in GKV decreases towards the outer region of ρ > 0.50 despite the monotonically increasing tendency of the linear growth rate shown in figure 6(a). Indeed, the stronger zonal flow generation from the TEM turbulence is associated with the transport reduction in the outer region. Consequently we see a moderate transport shortfall with relatively larger deviations from the experimental value for P i in comparison to those in the core region of ⩽ ρ 0.50. The parameter sensitivity of the transport shortfall associated with the strong zonal flows are discussed in section 4.3, as well as the flux matching evaluation of the profile prediction accuracy.

Ion and electron temperature gradient dependencies and multiple flux matching
As pointed out in [11], since the experimental plasmas often show a stiff response of temperature profiles against the auxiliary heating, the sensitivity of the calculated turbulent fluxes on the equilibrium profiles should be examined to make a more robust validation. In this section, we discuss dependencies of the turbulence simulation results on the numerical and physical parameters. Also, in the later part, the flux matching technique is applied for both P i and P e simultaneously to evaluate the prediction accuracy of the ion and electron temperature gradient profiles.
The numerical convergence with respect to the maximum wavenumber is shown in figure 8, where influences of the faster and finer TEM-unstable fluctuations ( figure 8(a)) are examined for the time evolution of ion and electron energy fluxes at ρ = 0.50 ( figure 8(b)) and ρ = 0.76 (figure 8(c)). We see that the higher maximum wavenumber leads to faster linear growth of the TEM modes for both cases. Their growths, then, saturate at much lower levels than that after the satur ation of the ITG-mode growth. Thus, the mean saturation levels of P i and P e in the steady state are not affected by the TEM modes with higher k y . This is also consistent with the fact that the lower wavenumber modes are dominant in the turbulent energy fluxes shown in figure 5.
The temperature-gradient dependencies of the ion energy flux P i , the electron energy flux P e , the turbulence energy W trb , and the zonal flow energy W zf (normalized by the total  In the present local flux-tube delta-f approach with the fixed background gradients, the steady temperature and density profiles in the power balanced state can be predicted by using the so-called flux matching technique [11,18].
A way of evaluating the prediction accuracy is to measure the deviation between actual gradients observed in the experiment and the input gradient values reproducing the experimental turbulent fluxes in the simulation. For instance, let , , , ,   Actually, the coefficient matrix is evaluated in a similar way to figures 9(a)-(c). In the previous work in [11], the ion heat flux is adjusted by only the ion temper ature gradient parameter, so that the electron heat flux still deviates from the experimental one. For more precise treatments, one needs to consider the × 3 3 matrix approach shown in equations (10)- (13), resulting from the coupling among the heat and particle fluxes at each radial position, but the accurate experimental evaluation of the particle flux is necessary.
In the present study, since there are no experimental data for Γ e as mentioned in section 4.2, the multiple flux matching technique is applied to the GKV simulation results (shown in figure 7) such that the radial profiles of P i and P e simultaneously match the experimental ones, i.e.
Note that although the linearized form with a matrix coefficient in equation (13) (14) and (15) One finds that the experimental fluxes are well reproduced except for the ion flux in the outer region of ρ > 0.58, where the prediction error for the ionand electron-temperature gradients for ⩽ ρ 0.58 is evaluated as less than ±30%. In the outer region, the weak dependence on the temperature-gradient shown in figure 9(c) prevents the flux matching within ±30% range of ( has also been found in [34]. Since it has been confirmed that the significant zonal flow generation in the TEM turbulence is related to the transport suppression in the outer region (see figure 6(b)), the turbulence simulation without zonal flow generations can provide us with useful insights, where the result is shown in the last row in table 2. Then, we see that the simulation results on both the ion and electron energy fluxes are strongly affected by the existence of zonal flows, and the transport shortfall in the ion flux becomes less significant. This is an artificial treatment, but clearly suggests the importance of the accurate zonal flow treatment in the simulation model. Actually, the generation of TEM-driven zonal flows and its impact on the turbulent transport show quite strong dependence on the equilibrium para meters ( / / R L T T , T e i e , etc) as shown in [35,36], and have not been fully clarified yet. Also, collisions with impurities from the wall and divertor and/or the toroidal magnetic field ripple, which are ignored in the present simulation model, may lead to relatively stronger zonal flow damping in the outer region. Some global effects of the radial propagation of the heat flux and the turbulence intensity [8,37,38] can also influence the zonal flow dynamics.

Summary
In this paper, quantitative comparisons of the ion and electron heat transport between gyrokinetic simulation and the JT-60U tokamak experiment are carried out by using a local gyrokinetic code GKV incorporating realistic magnetic geometry and fully-gyrokinetic electrons, where the ITG-and/or TEM-driven turbulent transport and zonal flow generations are investigated. In order to examine the prediction capability of the flux-tube gyrokinetic simulations, an L-mode plasma with sufficiently small ( / ) ρ ∼ * 1 500 is chosen for the quantitative comparisons, where the mean radial electric field and its shearing effect are also negligibly small. Nonlinear ITG and/or TEM simulations by GKV with kinetic electrons successfully reproduce the radial profiles of the ion and electron energy fluxes, which are relevant to the experimental values in the core region, whereas the adiabatic-electron case indicates relatively larger deviations. It is revealed that the zonal flow generation in the outer region (ρ > 0.58) with TEMdominated turbulence is much more significant than that in the core region with ITG-and ITG-TEM-dominated turbulence ( ⩽ ρ 0.58). Then, the transport shortfall for the ion energy flux appears in the outer region.
Extending the conventional ion heat flux matching by adjusting a single parameter of the ion temperature gradient, we performed a multiple flux matching for both the ion and electron energy fluxes. The temperature-gradient variations giving the matched energy fluxes in the core region are simultaneously determined for ions and electrons with a prediction error less than ±30%, while the weak temperature-gradient dependence of turbulent transport and zonal flows in the outer region prevents the flux matching. The Table 2. Temperature-gradient variations in multiple flux matching for P i and P e .  turbulence simulation without the TEM-driven zonal flows indicates the strong influence of both the ion and electron energy fluxes, where the ion transport shortfall becomes less significant.
To further improve the prediction accuracy of gyrokinetic simulations, and also to construct a credible reduced transport model, one needs the quantitative comparison of the particle flux between detailed experimental measurements and turbulence simulations. Also, more detailed analyses of the zonal flow dynamics are required, including its generation and damping mechanisms associated with the collisions with impurities, the toroidal magnetic field ripple, and some global effects, such as the radial propagation of the heat flux and the turbulence intensity. As the fluctuation measurements in JT-60U L-mode plasmas are available only for the very edge region, direct comparisons of the turbulence spectrum with the experimental measurement is out of scope at present, but they will be addressed for other appropriate experimental set-ups in future works.