Study of the L-mode tokamak plasma “shortfall” with local and global nonlinear gyrokinetic $δf$ particle-in-cell simulation

Gyrokinetic simulations are widely used to study tokamak plasmas and the models are now comprehensive enough that the simulation predictions can be directly compared to experimental results. Despite the success of the validation and verification studies using gyrokinetics, there are still open questions remain, such as the “transport shortfall”.^1,2 In GYRO³ simulations of a DIII-D L-mode discharge dominated by only micro-instabilities, i.e., no magnetohydrodynamic (MHD) instabilities in the model, the ion heat flux transport, Q_i, from simulation agrees well with the experimental result predicted by the transport code ONETWO at the inner core for $ρ=0.5$⁠, where ρ is the square root of the normalized toroidal flux. However, the simulation Q_i is much smaller than the experimental result at outer core, for $ρ=0.75$⁠. It seems this underprediction cannot be explained by the uncertainties of the experimental parameters used in the simulations. The problem has even raised questions on the validity of the gyrokinetic model,^4,5 and also raised questions regarding whether ONETWO is a good model for experimental transport predictions.⁶ 

On the other hand, a GYRO simulation of the outer core region of an Alcator C-Mod L-mode observes underprediction for Q_e but not for Q_i.⁷ In another systematic transport analysis of Alcator C-Mod L-mode using GYRO at different radial locations for two discharges differing in terms of the amount of auxiliary ion heating (termed as low and high power discharges),^8,9 no DIII-D L-mode type Q_i shortfall is found whatsoever and Q_e underprediction is observed in the case of the low power discharge. This Q_e underprediction is likely due to the lack of the high k electron scale physics in the simulations. A more recent multi-scale simulation including electron temperature gradient (ETG) scale fluctuations for these cases^10,11 finds the Q_e underprediction is diminished. Furthermore, simulations of the L-mode plasmas in the ASDEX Upgrade tokamak using the nonlinear gyrokinetic code GENE^12,13 have shown good agreement with experimental results and there is no shortfall within the uncertainties of experimental profiles in both the inner and outer core regions.¹⁴ Then, a new study of the same DIII-D shortfall case in Refs. 1 and 2 using the GENE code shows no shortfall by varying the pressure profiles.¹⁵ 

In this paper, we address the shortfall problem with another gyrokinetic code, GEM.^16,17 Unlike the GYRO and GENE code which are both continuum (Eulerian) models, GEM is a $δ f$ particle-in-cell code, and it has been well benchmarked with continuum codes such as GYRO and GS2 in different studies.^18,19 A particle code provides a critical verification of the previously reported nonlinear calculations using a completely different numerical approach. In Ref. 2, an earlier version of GEM has been used to simulate the same DIII-D shortfall case and has similar results as GYRO. An alternative flux-tube model has been implemented recently using a complete local equilibrium solution of the Grad-Shafranov equation,¹⁹ and we will apply this model in the current study.

Additionally, the GYRO and GENE models used in previous shortfall studies are local, i.e., flux tube simulations and the heat flux transport might be affected by this restriction. Adding edge physics, such as turbulence spreading^20–24 and avalanches propagating from the edge,²⁵ could increase the transport in the outer core region. GEM is an intrinsically global code and it has been successfully used to study both core^26–28 and edge physics.^29–31 Here, in addition to flux tube simulations, we also apply the global GEM model to study the edge turbulence spreading effects on the outer core transport.

The purpose of this paper is to document a fairly comprehensive GEM study of the relevant shortfall cases, and explore the mechanism of global edge turbulence spreading that has the potential to explain the shortfall. We first analyze the L-mode discharges from DIII-D and Alcator C-Mod tokamaks of previous GYRO studies in Refs. 1, 2, 8, 9, and 32, and compare ion and electron heat fluxes obtained from our simulations with those estimated in experimental measurements. In addition to these nonlinear results, we will identify the dominant linear instabilities at different radial locations across the minor radius of the two tokamaks. A global calculation is then carried out to examine the global profile effects and the effect of edge turbulence on the shortfall for the DIII-D L-mode discharge. We find that GEM has reconfirmed the GYRO results, i.e., the shortfall exists in the DIII-D case and there is no DIII-D L-mode like shortfall in Alcator C-Mod case, and although the edge physics has significant effects on the outer core transport, the effect is not big enough to explain the shortfall.

In this section, we document the $δ f$ particle simulation model in the GEM code. The main features include electromagnetics with full kinetic electron physics, $δ f$ particle in cell method, $p∥$ formulation, a split weight scheme for electrons, high β Ampere algorithm, a field line following coordinate system covering $0≤θ<2π$ and global profile effects. The model includes collisions, equilibrium flow, an arbitrary shaped tokamak equilibrium, and impurity species. In terms of the canonical momentum $p∥s=v∥s+qsms〈A∥〉$⁠, where s stands for particle species, the gyrokinetic equation is given as

where,

The guiding center velocity is given as $vGs=v∥sb̃+vds+vE$⁠, where $b̃=b+〈δB⊥〉/B, vE=〈E〉×b/B$⁠, and for low β plasmas

Note that in Eq. (1), $C(fs)$ is the collision operator. It is considered that ion collision operator $C(fi)=0$⁠. The Lorentzian operator is considered for electrons so that $C(fe)=CL(fe)$⁠, where

with $λ=v∥/v$ is the pitch angle. Also

where

The equilibrium distribution function is considered to be Maxwellian for all species, given by, $f0s(x)=n0s(x)/((2π)1.5vTs3(x))e−εs/msvTs2$⁠, where $εs=ms(v⊥s2+p∥s2)/2$ is the total energy for species, s. Note that $εs$ is a function of $p∥s$ instead of $v∥s$⁠.

A $δ f$ scheme is considered for ions, for which the evolution can be expressed as

A split weight scheme^33,34 is implemented for electrons, where the adiabatic component of the perturbed distribution function is split into two parts, as $fe=f0e−ϵgϕ(∂f0e/∂εe)+h$⁠. Then the perturbed electron distribution function evolves as

In the last two equations $ε̇s=μsvGs·∇B+msp∥sṗ∥s$⁠. The weight equations for ions and electrons w_i and w_e can be derived from the evolution equations for $δ fi$ and h and are as follows:

The electrostatic potential $ϕ$ is calculated using the gyrokinetic Poisson equation

where $ϕ̃=∑kΓ0(k⊥2ρti2)ϕkeik·x, ϕ=∑kϕkeik·x$⁠, and $dv=v⊥dv⊥dv∥dξ$⁠, with ξ and $ρ$⁠, respectively, are the gyroangle and the vector from particles' gyrocenter to the actual positions of the particles. The vector potential, similarly, is given by Ampere's law

It is to be noted that $∂ϕ/∂t$ required in the calculation of $dh/dt$ and $dwe/dt$ in Eqs. (3) and (5) can be obtained by taking time derivative of the Poisson equation expressed in Eq. (6).

The field line following coordinate system has been used in the model. They are given as, $x=r−r0, y=r0q0(∫0θdθq̂(r,θ)−ζ), z=q0R0θ$⁠, where $(r, θ,ζ)$ are the toroidal coordinates, R₀ and r₀, respectively, are the major radius at the magnetic axis and the minor radius at the center of the simulations box, q₀ is the safety factor at r₀, and $q̂(r,θ)=B·∇ζ/B·∇θ$⁠. Thus, $(x, y)$ labels the field line while z is the coordinate along the field line. The particle motion is given by $ẋ=v→G·∇x, ẏ=v→G·∇y$⁠, and $ż=v→G·∇z$⁠. Defining m_p = proton mass and $T0i$ = ion temperature, $vu=T0i/mp, xu=vu/ωci, ωci=eB0/mp, tu=1/ωci$ the normalization scheme can be seen in Tables I and II. For further details of the flux tube and global model of GEM readers are referred to Refs. 16 and 17.

We first study the “standard” DIII-D L-mode shortfall case^1,2 discharge 128913. The parameters corresponding to this case are obtained from Refs. 1 and 35 and delineated in Table III for two radial locations $ρ=0.5$ and $ρ=0.75$⁠, in which R is the major radius and a is the minor radius, R/a is aspect ratio, κ is elongation, $sκ=rdlnκ/dr$⁠, δ is triangularity, $sδ=rdδ/dr$⁠, Δ is Shafranov shift, q is the safety factor at the given values of ρ, $ŝ$ is the shear, T_e and T_i are the temperatures of electron and ion, $ρ∗=ρs/a, ρs=cs/ωi, ωi=eB0/mi$⁠, ν_e is the collision frequency, and $βe=2μ0neTe/B2$ is the electron β. The pressure gradient terms are denoted by $R/Lns=−Rd lnn0s/dx$ and $R/LTs=−Rd lnT0s/dx$ for species $s=e,i$⁠. $γE×B$ is the equilibrium radial electric field shearing rate. The values of β_e and ν_e are estimated using the relations $βe=neTe/(B2/2μ0)$ and $νe=nee4lnΛ/4πϵ02m2vthe3$⁠.

Despite the extensive nonlinear analysis of these two radial locations in previous studies that have found the ion heat flux transport shortfall at $ρ=0.75$ but no shortfall at $ρ=0.5$⁠, a complete linear analysis is not yet available. Here, we first run linear flux tube simulations to identify the dominant instabilities at the two locations. For the sake of simplicity, we remove the effect of E_r shearing rate in these linear simulations. The various linear scans are discussed below in tandem.

Figure 1 shows the dispersion relations obtained from linear GEM simulations for $ρ=0.5$ and $ρ=0.75$⁠. The real frequency and linear growth rate of the dominant instability are plotted as functions of $kθρs$⁠. The two locations have similar results, and we identify the dominant instability to be the ion temperature gradient (ITG) mode. The real frequencies in Figs. 1(a) and 1(c) are negative, corresponding to the ion diamagnetic direction in GEM's convention. The corresponding growth rates are shown in Figs. 1(b) and 1(d). Furthermore, a β_e scan shows that the mode is electrostatic, and a collision frequency scan shows the mode is independent of collision frequency. These are all typical ITG characteristics.

For a further confirmation, we can identify the main drive that provides free energy to these modes at the two locations. We carry out electron and ion temperature gradient scans as represented by $R/LTe$ and $R/LTi$⁠. In Fig. 2, the electron temperature gradient is varied by 10% around the experimental value and apparently the growth rates are insensitive to $R/LTe$ at both locations. In the ion temperature gradient scan shown in Fig. 3, the growth rates increase linearly with $R/LTi$⁠. Therefore, the mode is indeed driven by the ion temperature gradient.

Nonlinear flux tube simulation results of the ion heat flux transport of the two locations are shown in Fig. 4. The Q_i values are time-averaged over the steady state after saturation that included time window of width at least $1.5×105$ in $tωci$⁠. The experimental values (circles) are obtained from Ref. 2 for comparison. The simulation results using the “standard” experimental parameters are smaller than the ONETWO prediction of the experimental results at both locations. At $ρ=0.5$⁠, the simulation result is about half of the experimental value, while at $ρ=0.75$⁠, the simulation result is much smaller, and hence, the shortfall. We note here that our simulation result at $ρ=0.5$ is smaller than the old GEM results presented in Ref. 2, this is because we have updated the flux tube simulation algorithm and collision operator since then. The general trend that $ρ=0.75$ has a significant shortfall, however, is consistent with previous studies.^1,2

Could this under-prediction by simulation be due to uncertainties of the measured plasma parameters? Since the dominant instability is ITG, the Q_i flux is most sensitive to ion temperature gradient. Here, we increase the value of $R/LTi$ by 10%, 25%, 30%, and 50% to test the dependence of ion heat transport on $R/LTi$⁠. These new results are also shown in the same figure. At $ρ=0.5$⁠, the simulation value matches the experimental one when $R/LTi$ is increased by 25% (diamond); while at $ρ=0.75$⁠, the simulation result for ion heat transport agrees with the experimental value when $R/LTi$ is increased by 30% (square). These ranges are apparently beyond the error of the experimental measurement.

Next, we consider the two Alcator C-Mod cases studied in Refs. 8 and 9. The two cases are labeled as “low power” (LP) and “high power” (HP) characterized by the amount of auxiliary heating of ions. The physical parameters obtained from Ref. 8 are listed in Tables IV and V for two locations, $r/a=0.5$ and $r/a=0.8$⁠, where r is the Miller parameter that labels the flux surface and a is the minor radius. We note that in Ref. 8 the driving gradients and suppression terms are slightly adjusted to match the heat flux from GYRO simulation with the experimental heat flux obtained from TRANSP³⁶ within the experimental uncertainty, and thus the simulation study is termed as the “flux matched simulation.”⁸ The GEM flux tube simulations here use these same parameters.

For flux tube linear analysis, we first carry out a $kθρs$ scan for real frequency and growth rate at the two locations for both the low power and high power cases without E_r. It appears the low power case is dominated by the ITG at both locations, as shown in Fig. 5. At $r/a=0.5$ in Figs. 5(a) and 5(b), the real frequency is in the ion diamagnetic direction, and the growth rate increases with $kθρs$ because of magnetic drift resonance which depends on k, then reaches a maximum value and eventually rolls off due to the finite Larmor radius effect. At $r/a=0.8$ as in Figs. 5(c) and 5(d), the real frequency is in the electron direction at lower k and goes into the ion direction at higher k. The growth rate continues to increase with increasing $kθρs$ until $kθρs≈0.7$ for a really high n = 120 and tends to saturate. The ion mode is clearly dominant.

The high power case shown in Fig. 6 has a different story. At $r/a=0.5$⁠, the real frequency is first in the ion direction with a small amplitude when $kθρs$ is smaller than 1.1, then jumps to the electron direction with a larger amplitude for bigger k and then increases with k, shown in Fig. 6(a). It appears the dominant instability here is an ion mode for low k and an electron mode for high k, and the ion mode has higher growth rate, as shown in Fig. 6(b). Since for the high power case $R/LTi$ and ν_e are decreased and T_e is increased, we identify this electron instability as the trapped electron mode (TEM). At $r/a=0.8$ shown in Figs. 6(c) and 6(d), the real frequency is in the electron diamagnetic direction at low k then switches to the ion diamagnetic direction and increase with k. The growth rates increase with $kθρs$⁠, asymptoting at higher k. The absence of a strong finite Larmor radius effect could be due to the presence of trapped electron modes in this discharge. Thus, the high power discharge shows a dominant TEM not found in the low power cases. Also, radially outer locations are more susceptible to TEM activity, partly because of the increased fraction of trapped electrons as one moves from the core towards the edge of the tokamak.

β_e scans show the growth rate is independent of β_e in all cases, indicating the modes are electrostatic. This is expected because the β_e is usually very small in Alcator C-Mod plasmas. A collision scan is also carried out for the two cases. The experimental value of ν_e is increased and decreased by 10%. It is observed that the variation in the growth rate with collision frequency is weak around the experimental value.

To identify the drive of the instabilities, we look at the sensitivity of the mode growth rates with respect to the temperature gradients of electrons and ions. First, we vary the electron temperature gradient by $±10%$ keeping the value of ion temperature gradient $R/LTi$ fixed. The corresponding growth rates are calculated and displayed in Fig. 7 for all the cases. The mode growth rate is virtually independent of $R/LTe$ except at $r/a=0.8$ for the high power case, in which γ increases linearly with $R/LTe$⁠, confirming that the instability here is the TEM. In the corresponding $R/LTi$ scan as in Fig. 8, except for $r/a=0.8$ in the high power case which is also the TEM dominant, the dominant instabilities are driven by the ion temperature gradient for the chosen value of $kθρs$⁠.

For nonlinear flux tube results, the time averaged ion and electron heat fluxes Q_i and Q_e from the GEM simulations (stars) are compared to the experimental values (squares) obtained from Ref. 8, shown in Fig. 9. The simulated heat fluxes via ion and electron channels are in good agreement with those from experiments evaluated using the TRANSP³⁶ code except in the low power case, where electron heat flux is lower than the experimental range at both radial locations. This is in good agreement with GYRO results in Ref. 8.

Thus, in contrast to the results from DIII-D L-mode plasmas where the ion heat flux is low (i.e., the shortfall) in the outer radial locations, the present analysis indicates no such Q_i shortfall exists in these Alcator C-Mod L-mode discharges. However, the electron heat flux is under-predicted at all radial locations in the low power discharge. This could be explained by the fact that the low power discharge is characterized by higher $a/LTi$⁠, higher collisions and lower electron temperatures⁸ and hence electron modes are less dominant, resulting in a lower electron heat flux. In fact, including high k, ρ_e scale resolution in simulations of the low power discharge showed better agreement between experimental electron heat flux and the simulation results.^10,11 The including of high k turbulence may, in principle, also enhance the ion heat fluxes in those simulations, but the effect is within the experimental error bars.

In this section, we study the global effects on the local transport level of DIII-D discharge 128913. Flux tube simulations, though with a finite radial domain, are local studies that assume the pressure and magnetic equilibrium parameters and their gradients are homogeneous across the simulation domain. The only exception is the safety factor q, which is implemented using a linear profile matching the local values of q₀ and $ŝ$⁠, to maintain the magnetic shear. In GEM's global model,²⁹ the realistic pressure and magnetic equilibrium profiles are implemented including the experimental q profile. With the global model different physical instabilities may appear at different radial locations and couple with each other within the simulation domain, and turbulence could grow and spread across the domain.

For the present study, the experimental magnetic equilibrium profile (i.e., the eqdsk file) and pressure profile are used directly in GEM, and all radial pressure gradients are then calculated numerically within the simulation. As a result, the Miller parameterization³⁷ and other simulation parameters used in our global simulations are slightly different from the “standard parameters” used in other studies of the shortfall problem.^1,2 We first compare flux tube GEM simulation results using these parameters to the experimental prediction of ONETWO² in Fig. 10. Here, the simulation value of Q_i at $ρ=0.5$ is quite close to the one we get with the standard parameters, though now smaller than the experimental prediction. At $ρ=0.75$ there is apparently a Q_i shortfall. In addition to $ρ=0.5$⁠, the simulations show a second peak of Q_i at $ρ=0.8$⁠, and the dominant instability is ITG at both locations. The TEM appears to dominate around $ρ=0.9$⁠, with a very low Q_i. At $ρ=0.6$ to 0.65, however, the simulations are almost completely stable.

Although the simulation results at all locations are smaller than the ONETWO prediction, the trend of instabilities is consistent with the pressure profile. Figure 11 shows the variation of η_i and η_e. η_i peaks around $ρ=0.5$ and $ρ=0.8$⁠, the same locations of ITG dominance, and η_i is below 2 around $ρ=0.65$⁠, where ITG is stable. For $ρ>0.9$⁠, η_e is big but η_i is very small, and hence we have a small Q_i. In contrast to the smooth experimental fluxes, the local simulation fluxes have ups-and-downs. One possible mechanism that may reconcile the local simulation and experimental results is that turbulence spreading from the strongly unstable regions could make up for the shortfall in the stable regions. Especially, strong turbulence from the steep gradient region of the edge, which is not included in any local studies, may help increase the overall transport of the outer core region. We need to test this possible mechanism with global simulations.

The GEM global simulations here focus on the outer core to edge region, and cover the radial domain of $0.699≤r/a≤0.999$⁠, or $0.641≤ρ≤0.998$ in terms of ρ. The simulations are electromagnetic and collisional, with gyrokinetic deuterium and carbon ions and drift kinetic electrons. The particle-in-cell scheme has 128, 128, and 32 grids in L_x, L_y, and L_z in terms of the field-line-following coordinates¹⁷ and 32 particles per species per cell. The simulation time step is $Δt=1/Ωi$⁠, where Ω_i is the proton gyro-frequency at the core.

Global linear results are consistent with the flux tube simulations discussed in Sec. III. The dominant instability appears at the edge steep density gradient region and is identified as the TEM. In a simulation with a smaller radial domain that includes the edge region only, we have the same linear growth rate as in the simulation with the full domain. In a simulation that does not include the edge, we have observed the ITG instability that appears in the outer core region. Also same as in the flux tube runs, high-n modes (n > 50) dominate both the ITG and the TEM. In nonlinear global simulations including the low-n modes, we have observed a geodesic acoustic mode (GAM)-like instability that has a saturation level beyond the regime where the $δ f$ method is valid, hence here we keep the $n≥20$ modes only to exclude any low-n instabilities, while the $kr>5·2π/Lx$ zonal flow is included.

To test the edge effects on the core transport, we have two nonlinear simulations: one with the original experimental profile and the other with flattened profiles at the edge. Figure 12 shows the n_e, T_i, and T_e profiles for the two simulations. In the flattened case, the pressure gradient terms are multiplied by a fattening function $exp[−((r/a−0.85)/0.05)2]$ for $r/a>0.85$⁠, and therefore the ITG dominant region, which is inside $r/a<0.85$⁠, is not affected by the flattening. Without the pressure gradient drive, the flattened case is almost linearly stable in the edge region.

The nonlinear evolution of the two simulations is shown in Fig. 13. The root mean square (RMS) value of the overall volume-averaged perturbed electrostatic potential is plotted to indicate the turbulence level. In the case with the original edge, the linear growth is very strong and the simulation saturates very quickly at a high level. In the case with the flattened edge, the linear stage lasts for a very long time and the simulation eventually saturates at a low level. Apparently, the edge instability has significantly increased the total turbulence.

Figures 14 and 15 show the process of turbulence spreading of the two cases. Snapshots of the perturbed electrostatic field $δϕ(x,y)$ turbulence were taken from three times during the evolution. In Fig. 14 with the full edge pressure gradient drive, the TEM instability at edge dominates the whole simulation domain for a long time, as shown in the top panel for $tΩi=2000$⁠. The ITG instability at $0.7<ρ<0.8$ grows very slowly and then turbulence appear to peak at two locations, one at the edge from the TEM and other at the outer core from the ITG, as shown in the second panel. During this period, the edge instability has already saturated and its turbulence is spreading inward towards the core. This phenomenon is quite similar to the turbulence spreading studied in Ref. 22 in which ITG turbulence spreads from the outer core region (around r/a = 0.7) to the core. The TEM and the ITG start independently and grow locally for a while. As the ITG instability saturates, its turbulence also begins to spread both inward and outward. Turbulence from the two instabilities goes across the linearly stable region around $ρ=0.85$ and eventually spreads all over the simulation domain, with visible zonal flow structures, shown in the bottom panel.

In the case with the flattened edge pressure profile as shown in Fig. 15, the edge region is stable and the only instability is the ITG at the outer core. After a very long linear growth, the ITG turbulence spreads across the linearly stable region into the edge. Similar to Ref. 22, large amplitude turbulence is still largely concentrated in the linearly unstable region even after the spreading.

Figure 16 shows the time-averaged turbulence intensity, which is represented by RMS $δϕ(x)$ averaged over y and z, for the two cases at the quasi steady state after saturation. In the case with the original edge pressure profile (solid line), turbulence slightly peaks at the TEM dominant edge region, but it has spread throughout the simulation domain. In the case with the outer core ITG only (dashed line), turbulence reaches the edge region but it still largely concentrates in the outer core. For everywhere including the ITG dominant region of $0.7<ρ<0.8$⁠, the turbulence intensity with the edge instability is clearly higher than that without the edge instability. Therefore, we conclude that the edge instability has increased the fluctuation level at the outer core region. This is indeed a significant global effect that it has even greatly increased the turbulence level at inner boundary of the simulation domain near $ρ=0.65$⁠, which is linearly stable in local simulations.

How does the turbulence spreading from the edge affect the outer core transport? In Fig. 17, we compare Q_i at the saturation stage from our global simulations to the flux tube simulations and the experimental results. Results from the case of the global simulation without the edge drive qualitatively agree with that of the flux tube simulations. Compared to the flux tube results, the global results are higher for $ρ<0.8$ and lower for $ρ>0.8$⁠, but the trend is the same. The transport peaks around $ρ=0.8$ and diminishes in both ends. When the edge drive is included as in the results shown in the solid line, the transport at $ρ=0.75$ is strongly increased to more than twice of the results from flux tube simulations, apparently an effect of turbulence spreading from edge. The TEM edge instability also increases Q_i in the edge region, where the ITG is stable. The edge effect, however, is not strong enough to make up for the Q_i shortfall, especially in the ITG stable regions. The carbon impurity heat flux is included in the calculation of Q_i here but its contribution is very small.

In Fig. 18 for Q_e, edge turbulence also increases the transport at the outer core, but the increase is not as big as for Q_i. With the pressure gradient drive there is a second peak of Q_e at the edge, where Q_e is almost zero in the case with the flattened edge.

Figure 19 shows the ion profile relaxation as represented by $δ ni$ and $δ pi$⁠. In $δ f$ simulations for $f=f0+δ f$⁠, we assume $δ f≪f$⁠, and Fig. 19(a) indicates that this conditions is well satisfied even after a very long time on nonlinear run. Since $δ ni/ni0$ is within 1% and $δ pi/pi0<4%$⁠, our $δ f$ simulations are valid and this scale of relaxation is physical. Figs. 19(b) and 19(c) show that the relaxation tend to decrease the pressure gradients at the location with strongest linear instability (around $ρ=0.75$⁠), and as a result increase the gradients at nearby locations. Lower pressure gradients correspond to lower instability, therefore profile relaxation is a mechanism that lowers the transport at the most unstable region. This phenomenon is consistent with the physical intuition that transport tends to smooth the particle and energy distributions. Taking into account for this effect, one would expect simulations without profile relaxation should have higher transport levels than the experimental values to make up the shortfall, yet we did not see this.

We see similar phenomenon in a different DIII-D discharge as well. Figure 20 shows the profiles of Discharge 145781 as compared to the classical shortfall case of 128913. The 145781 has higher temperature gradients at core, and we may expect the edge turbulence spreading should be less effective in increasing the core transport. This is exactly what we see in the results shown in Fig. 21. The core ITG grows quickly in the flattened edge case, and the increase of the saturation level with the edge pressure gradient drive is not big. At the outer core, Q_i is still doubled by edge turbulence spreading, but Q_e is almost unaffected. It is clear that the edge TEM instability has a smaller effect here.

These global results indicate that although the TEM turbulence spreading from the edge can increase the fluctuation intensity, making it nearly constant throughout the simulation domain, its effect on the transport is not as striking. Heat flux transport is still largely determined by local pressure gradients. The global effects of the edge TEM are not enough to explain the transport shortfall. Nevertheless, we note from the Discharge 145781 results that the enhancement of the core transport from edge turbulence is related to the relative strength of the edge instability. A stronger or additional edge instability that is not modeled in this study, e.g., any low-n or MHD modes, may increase the core transport more. In fact, it has been reported recently that GAMs are observed in DIII-D L-mode plasmas.³⁸ 

Only electron collisions are included in the present simulations, since the electron collision frequency is bigger than that of the ions and important for linear physics.³⁹ Zonal flows are driven by turbulence but also regulated by nonlinear damping.⁴⁰ It has been reported that ion collisions could also impact turbulence and transport via zonal flow damping.^41,42 In the present case, we assume that the nonlinear damping is stronger than the collisional damping, and hence neglect the ion-ion collision. We note that collisional damping could occur over larger time scales and therefore should be considered in future shortfall studies.

To compliment recent continuum validation efforts,^1,2,8,14 we have analyzed L-mode shots from both the DIII-D and Alcator C-Mod tokamaks using the nonlinear gyrokinetic particle code GEM.

For the DIII-D L-mode plasmas, the ion heat flow is comparable albeit smaller than the experimental prediction at $ρ=0.5$⁠, but is much too low at location $ρ=0.75$⁠. This is in agreement with earlier GYRO results.² Linearly, both the core and edge are dominated by the ITG mode, and there is little variation of the growth rate with $R/LTe$⁠, β_e, and ν_e. By increasing the value of $R/LTi$ by 25% and 30% at the locations $ρ=0.5$ and $ρ=0.75$⁠, we find satisfactory agreement with the experimental heat fluxes.

For the Alcator C-Mod L-mode case, the ion heat fluxes are in satisfactory agreement with earlier GYRO results as well,⁸ and also within the range of experimental error. The electron heat fluxes are also in fair agreement except for the low power discharge where the electron heat fluxes are much lower than the experimental values. Linear analysis illustrates that the most unstable modes are ITG driven for the low power discharge, but both the ITG and TEM are present in the high power discharge case.

It was hypothesized in the earlier studies¹ that global simulations may help reduce the discrepancy between experiment and nonlinear gyrokinetic simulation. We carry out a global simulation for the DIII-D L-mode shot with slightly modified profiles with and without the effect of edge turbulence. It appears that the core is dominated by the ITG turbulence in the absence of edge. However, when edge is included, the TEM turbulence spreading from the edge to the core enhances the ion heat flux level in the core substantially, but the overall ion heat flux is still lower than the experimentally predicted level. Therefore, we conclude that the effect of TEM edge turbulence alone is not sufficient to explain the observed shortfall in the DIII-D L-mode shots.

This material is based upon work supported in part by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Award Nos. DE-FC02-04ER54698, DOE-SC0008801, DE-FG02-07ER54917, and DE-FG02-08ER54954. DIII-D data shown in this paper can be obtained in digital format by following the links at https://fusion.gat.com/global/D3D_DMP. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), which is a DOE Office of Science user facility.

We thank C. S. Chang, S. Ku, J. Lang, and R. Hager for helpful discussions. We also acknowledge fruitful discussions on “shortfall” with A. E. White and Tobias Görler.