Validation of nonlinear gyrokinetic simulations of L- and I-mode plasmas on Alcator C-Mod

With a new generation of fusion devices on the horizon, the need for steady state scenarios operating with high energy confinement, but low impurity confinement (to avoid fuel dilution and radiative losses), is apparent. The I-mode¹ is a high energy confinement regime of plasma operation that is regularly obtained on the Alcator C-Mod tokamak.² I-mode plasmas have electron and ion temperature pedestals similar to those of H-mode, but have a density pedestal similar to L-mode. I-mode plasmas achieve energy confinement times comparable to or exceeding H-mode, without as strong of a confinement time degradation with increased input power. Despite a natural lack of ELMs,³ I-mode does not suffer from H-mode-like impurity accumulation, even in a metal-walled machine such as Alcator C-Mod.^4–6 

I-mode is generally run with the ion $B × ∇ B$ drift away from the active X-point (unfavorable $∇ B$ drift), which enables a more robust access to the I-mode confinement regime.⁵ I-mode appears in a power window between L-mode and H-mode, and it has been found that this power window opens up significantly at a higher magnetic field.⁷ Alcator C-Mod can operate between 2.1 and 8 T. I-mode is accessible below 3 T, but even small power increases generally lead to a H-mode transition. The operating window is much wider at 5.4 T (C-Mod's standard field), and no I-H transitions have been observed at 8 T, even with full input power.⁷ Despite this smaller operating window at a lower magnetic field, I-mode has also been observed on both ASDEX Upgrade^8,9 and DIII-D.¹⁰ 

The properties of I-mode outlined here make it quite a favorable candidate for an operating regime on future fusion devices, such as ITER or another high field burning plasma experiment. For this reason, it is vital to understand the core transport in I-mode plasmas (both experimentally and with models) and, in doing so, enable accurate predictions for future devices. Significant progress has been made in recent years toward validating gyrokinetic and gyrofluid models in L- and H-modes, but relatively little work has focused on I-mode (see Ref. 11 for experimental work and Ref. 12 for comparisons with gyrokinetic simulations). This paper focuses on validation of the gyrokinetic model as implemented in the code GYRO¹³ in L- and I-mode plasmas, as well as making comparisons between experiment and the reduced gyrofluid model TGLF.¹⁴ 

This validation study is further motivated by the relatively recent discovery of the importance of “multi-scale” effects in gyrokinetic simulations.¹⁵ To date, most work with gyrokinetic validation has utilized the “ion-scale” simulations, which capture the effects of long-wavelength (⁠ $k θ ρ s < 1.0$⁠) turbulence such as the ion temperature gradient (ITG) mode and the trapped electron mode (TEM), ignoring contributions from electron scales (⁠ $k θ ρ s > 1.0$⁠) such as the electron temperature gradient (ETG) mode. We define $k θ$ as the poloidal wave number and $ρ s = c s / Ω c , i$ as the ion gyroradius evaluated with the ion sound speed (where $c s = T e / m i$ is the ion sound speed, and $Ω c , i = e B / m i$ is the ion gyro-frequency).

Recently however, it has been demonstrated that one cannot simply add the effects of the two single scale simulations linearly in many plasma conditions, and that cross-scale coupling between the different types of turbulence plays an important role in the overall plasma dynamics.^15–18 Proper characterization of these effects requires “multi-scale” simulations (with realistic electron mass), which simultaneously capture the electron and ion scales. It has been shown that in an Alcator C-Mod L-mode plasma, the inclusion of multi-scale effects can resolve discrepancies between values of ion heat flux, electron heat flux, and perturbative thermal diffusivity from experimental measurements and those calculated from ion-scale GYRO simulations.¹⁶ 

Unfortunately, these multi-scale GYRO simulations are incredibly computationally expensive, requiring around 100 million CPU hours for a validation study. Due to these massive resource requirements, the present study will use ion-scale GYRO simulations for the gyrokinetic validation study, and then will investigate the multi-scale effects with the TGLF gyrofluid code, which has recently been tuned to take these into account.¹⁹ 

The remainder of this paper is organized as follows. Section II will outline the experimental methods used in this validation study, as well as present some of the important measurements made in both L- and I-mode. Section III then describes the global, nonlinear, ion-scale GYRO simulations performed as part of this study and compares the GYRO results to the experimental measurements. Section IV investigates the possible impact of multi-scale effects on the plasmas described in this study through comparisons with TGLF simulations. Finally, Section V discusses the implications of the results of this validation study, both for the comparison of L- and I-mode plasmas and the validation of GYRO and TGLF, and presents open questions to motivate further research in this area.

Alcator C-Mod is a high field (2.1 – 8.0 T), high performance, compact (a = 0.22 m, R = 0.67 m), diverted tokamak with high-Z plasma facing components.² All of the experimental data discussed in this paper was taken from a single discharge (shot 1120921008), which undergoes a transition from L-mode to I-mode. The plasma parameters are: $n e = 0.7 − 0.8 × 10 20 m − 3 , I p = 1.1 MA , B t = 5.4 T$⁠, and $q 95 = 3.2$⁠. This shot was run with the ion $B × ∇ B$ drift away from the active X-point (unfavorable $∇ B$ drift). In addition to ohmic heating, auxiliary ion cyclotron range of frequency (ICRF) heating is applied, initially at 1.6 MW, and then stepping up to 3.5 MW, which initiates the transition to I-mode (moderate power for Alcator C-Mod operation).

Only times when the plasma was in a stationary L-mode or stationary I-mode were analyzed. The energy confinement time was approximately 30 ms (⁠ $H 98 , y 2 = 0.65$⁠) in L-mode and 40 ms (⁠ $H 98 , y 2 = 0.95$⁠) in I-mode in this discharge, so the analysis periods span multiple energy confinement times in both regimes. Time traces of this discharge are shown in Figure 1.

This particular discharge was designed to have a clean transition from L-mode to I-mode in the same plasma and to include the large range of diagnostics required for proper power balance calculation and comparison to gyrokinetic simulation. Electron temperature profiles were measured with a grating polychrometer (GPC) electron cyclotron emission (ECE) diagnostic.²⁰ Thomson scattering also measured the electron temperature profile, as well as the electron density profile.²¹ Ion temperature profiles and the toroidal rotation were measured with a high resolution X-ray crystal spectrometer (XRCS).²² 

Electron temperature fluctuations were measured with a correlation electron cyclotron emission (CECE) radiometer, described in Refs. 23 and 24. Tunable bandpass filters and a small beam spot size allow the CECE diagnostic on C-Mod to measure the local turbulent fluctuations below $k θ ρ s < 0.3$ in the region between $r / a = 0.7$ and $r / a = 0.9$ in a wide variety of plasma conditions. Three pairs of two channels measure fluctuations at different radial locations simultaneously.

For the purposes of this study, 150 ms of time averaging was used for both the L-mode and I-mode portions of the shot in question. This was done both to allow for sufficient signal to noise ratio for the CECE diagnostic, and to improve the fidelity of temperature, density, and rotation profiles for input to TRANSP and GYRO. For the remainder of this paper, the L-mode portion will refer to the time range from 0.85 to 1.00 s. Similarly, the I-mode portion will refer to the time range from 1.15 to 1.30 s. These time ranges were chosen such that they were relatively close temporally, to minimize unintended changes in plasma parameters, but also such that both time periods were fairly steady in time.

Magnetic field, plasma current, and the q₉₅ were very steady through the transition from L-mode to I-mode. The density did slightly rise with the increase in RF power, as seen in Figure 1, which must be kept in mind when examining the diagnostic data. This relatively small change in density is not uncommon in I-mode (due to the RF power increase), and is entirely distinct from the sharp rise in density associated with the transition to H-mode.

Power balance analysis for the above plasma was performed with the TRANSP code.²⁵ Input profiles of electron and ion temperature, electron density, and plasma rotation were measured with the diagnostics described above, and fits to the temperature and density profiles are shown in Figure 2. Also shown in Figure 2 are the relevant turbulence drive terms, the normalized temperature and density gradients. TRANSP then takes these and other inputs, and solves for the power balance equilibrium of the shot in question.

This study will focus primarily on the radial region well outside of the sawtooth mixing radius (⁠ $r / a > 0.6$⁠) and inside of the edge-pedestal region (⁠ $r / a < 0.9$⁠), as is indicated on the profiles shown in Figure 2. Evolution of the electron temperature profile tends to be dominated by the sawteeth inside of the sawtooth mixing radius. In addition, careful fits of the edge profiles are required for adequate analysis of transport in the edge, and since this study focuses primarily on core transport, such fits were not performed. Such edge analyses have been carried out in the past, such as in Ref. 4.

Sec. III will show the heat fluxes calculated with TRANSP, and will compare them to the results of GYRO simulations.

Uncertainties on the turbulence drive terms were approximately 22% for $a / L T e$⁠, 40% for $a / L T i$⁠, and 30% for $a / L n e$ over the analysis region, calculated from diagnostic and time uncertainty. In addition, uncertainties on all input profiles were propagated through the governing equations of TRANSP in order to obtain the uncertainties displayed on the heat flux profile outputs.

Even before a detailed quantitative analysis, the profiles shown in Figure 2 imply some qualitative differences between turbulence in L- and I-modes. Generically, one sees that all three turbulence drive terms plotted here (⁠ $a / L T e , a / L T i$⁠, and $a / L n e$⁠), decrease from L- to I-mode. Further analysis also shows that the turbulence suppression term (⁠ $γ E × B$⁠) increases (see Table II). These two effects combined suggest that turbulence should decrease in the core of I-mode compared to that in L-mode. In order to measure this turbulence, this study will utilize the CECE diagnostic described above.

The CECE diagnostic data is analyzed using relatively standard time history analysis techniques, which are outlined in Ref. 26. The coherence, γ, between neighboring radiometer channels is calculated and then used to calculate the temperature fluctuation level of the plasma at that location, $T ̃ / T$⁠.

We define the coherence as the cross-spectral density function normalized by the two autospectral density functions (see page 147 of Ref. 26)

where G_xy is the cross-spectral density function between channels x and y, and G_ii is the autospectral density function of channel i.

The total temperature fluctuation is then

where B_sig is the signal bandwidth (⁠ $f 2 − f 1$⁠), B_IF is the intermediate frequency bandwidth (the bandwidth of the bandpass filters used in the CECE system), and f₁ and f₂ define the frequency range over which to integrate the coherence. This is equivalent to Equation (7) in Ref. 27, rewritten in terms of γ_xy instead of the cross-powers and auto-powers separately. See Refs. 28 and 29 for earlier derivations of this and the following formulae.

We note here that both plasma time periods under consideration here are optically thick³⁰ at all three radial locations, and so the fluctuations in the CECE output can be interpreted as electron temperature fluctuations.

As described in Refs. 23, 24, 27, and 28, there exists an instrumental sensitivity limit for CECE diagnostics. This equation places a limit on $T ̃ / T$ below which temperature fluctuations cannot be differentiated from thermal noise. This limit is given by

where

Again, B_IF is the IF bandwidth, and B_sig is the signal bandwidth. In addition, B_vid is the video bandwidth (the bandwidth of the digitized signal), and $Δ t$ is the measurement time window. This instrumental sensitivity limit is based on the thermal noise inherent in the plasma (photon statistics).

In addition to the instrumental sensitivity limit, there also exists a statistical sensitivity limit, which is inherent in the statistical techniques used here (the integration of coherence). This limit is described on page 333 of Ref. 26 as the bias error. The limit is

where

in which n_d is the number of independent ensemble averaging windows.

It is the larger of the two sensitivity limits that determines the smallest measurable temperature fluctuation. Both of these limits depend on the time period over which the measurement is performed, and can be decreased by expanding this period. This study, however, was constrained in the time period due to the desire to look at an L-mode and an I-mode in the same discharge (shot-to-shot comparisons introduce other possible difficulties). For the parameters of this study, the instrumental sensitivity limit is 0.24%, while the statistical sensitivity limit is 0.69%. Therefore, only the measured temperature fluctuations above 0.69% are reported to be true temperature fluctuations.

The fluctuations from the L- and I-mode time periods of the shot analyzed in this paper are shown in Figure 3. The calculated $T ̃ / T$ values for each of the three radial locations at which measurements were made in both L- and I-mode are collected in Table I. This data reveals that the fluctuations decrease at all three radial locations (⁠ $ρ t =$ 0.80, 0.84, and 0.89, where ρ_t is the square root of the normalized toroidal flux) after the transition from L-mode to I-mode. In the two inner radii, the signal drops below the sensitivity limit, described above, in I-mode. This decrease in turbulence at all radii in this region is consistent with the decrease in turbulence drive terms and increase in suppression, as shown in Figure 2, and is also consistent with the previous work.^11,12 The turbulent spectra shown here will be compared to the results of gyrokinetic simulations in Sec. III.

Finally, the electron perturbative thermal diffusivity, $χ e p e r t$⁠, of the plasma was measured via the partial sawtooth heat pulse propagation method.³¹ The perturbative thermal diffusivity differs both in definition and value (in many plasmas) from the power balance thermal diffusivity, $χ e P B$⁠.³² The power balance thermal diffusivity governs the steady state transport of heat through the plasma, and is defined as

where Q_e is the electron heat flux, $∇ T e$ is the electron temperature gradient, and n_e is the electron density.³² Note that later in this paper, the subscript e may be dropped, but all diffusivities refer to electron thermal diffusivities.

On the other hand, the perturbative (also called incremental in some literature) thermal diffusivity governs the propagation of heat pulses and is defined as³² 

This definition reveals that $χ e p e r t$ is related to the profile stiffness of the plasma, where stiffness is the incremental change in flux for an incremental change in gradient, above the critical gradient.³³ 

In order to measure the perturbative thermal diffusivity, heat pulses generated by the partial sawtooth crash are tracked as they propagate radially outward in the plasma on the various channels of an ECE diagnostic. The use of partial sawteeth avoids contamination of the heat pulse data with the “ballistic effect” that affects full sawtooth-generated heat pulses.³⁴ The perturbative thermal diffusivity is then calculated using the “Extended-Time-to-Peak” method.^35,36 This relates the perturbative diffusivity to the velocity of the peak of the heat pulse, as well as the decrease in amplitude of the heat pulse as it moves radially outward.

The values of perturbative diffusivity measured with this method are inherently a radial average over the region in which the heat pulse propagates. For the study here, this is in the region between r/a = 0.64 and 0.84. This radial range was chosen to avoid the region too close to the partial sawtooth mixing radius, while also ensuring sufficient signal-to-noise ratio, as the heat pulse magnitude is too small near the pedestal to make reliable measurements of the pulse propagation.

For the plasma presented here, the measured perturbative thermal diffusivity was $χ e p e r t = 4.0 ± 0.6 m 2 / s$ in L-mode and $χ e p e r t = 4.3 ± 0.9 m 2 / s$ in I-mode. While these values are the same within experimental uncertainty for this plasma, past work has revealed that I-mode generally has a higher perturbative diffusivity than L-mode.³¹ 

This does not contradict the fact that I-mode is a high confinement regime, since the perturbative diffusivity measures the profile stiffness, and not the overall confinement. The power balance diffusivity is lower in I-mode than in L-mode (⁠ $χ e P B = 1.4 ± 0.4 m 2 / s$ in L-mode and $χ e P B = 1.0 ± 0.3 m 2 / s$ in I-mode), though again, not quite outside of uncertainty. Even though I-mode is a high confinement regime (⁠ $H 98 , y 2 = 0.95$⁠), one would not expect the factor of 2 decrease that one observes in $χ e P B$ after the transition from L-mode to H-mode. $χ e P B$ depends inversely on density and H-mode tends to have a strong density increase while I-mode does not.¹² 

This section compares the experimental ion heat flux, Q_i (from TRANSP), electron heat flux Q_e (from TRANSP), electron temperature fluctuations (from CECE), and perturbative thermal diffusivity, $χ e p e r t$ (see Equation (8), from partial sawteeth), to global (r/a = 0.65 – 0.9), nonlinear, ion-scale GYRO simulations. After presenting these comparisons, this section will summarize the conclusions that one can draw from these results.

GYRO is an Eulerian gyrokinetic code exhibiting high physics fidelity.¹³ All inputs to the simulations were obtained directly from the experiment. These simulations included three kinetic species (gyrokinetic ions and impurities, and drift-kinetic electrons), realistic geometry with Miller parameterization, rotation effects (including E × B shear), electron-ion and ion-ion collisions, and electrostatic turbulence (electromagnetic effects are ignored due to the very low β of this plasma). The simulation box size was approximately 105 by 120 ρ_s in the radial and binormal directions, and included 28 toroidal modes and approximately 500 radial grid points (for a grid spacing of $∼ 0.25 ρ s$⁠). These simulations captured the long wavelength turbulence (ITG/TEM) up to $k θ ρ s$ of approximately 1.3.

A summary of the inputs to GYRO at $r / a = 0.8$ for both L- and I-mode is given in Table II.

In order to perform global simulations using the GYRO code, benign buffer regions were employed. For the simulations reported here, these buffers were approximately 13 (8) ρ_s wide in the inner (outer) regions of the I-mode simulation and approximately 11 (8) ρ_s wide for the L-mode simulation. Default values for the buffer region sources were used, with the source annihilation rate, $ν s o u r c e = 0.1 c s / a$⁠. For more details on the buffer regions, see Ref. 37.

To determine the dominant turbulent drive terms, a linear growth rate analysis was performed, the results of which are shown in Figure 4 for L-mode (I-mode results are similar). At long-wavelengths, the ITG mode is dominant, which implies that $a / L T i$ is the dominant turbulence drive term in this regime. In order to obtain the ion heat flux-matched simulations, the experimentally measured $a / L T i$ was varied within the diagnosed experimental uncertainty (40%) to match the power balance value of the ion heat flux. As will be discussed later, it may be that in future validation studies, multiple input parameters will need to be varied simultaneously, but limited computational resources make a variation of only the dominant drive term most economical. This method is consistent with the past validation studies of gyrokinetic codes on Alcator C-Mod, where ITG is typically dominant in most L-, I-, and H-mode plasmas.^{11,12,15,38,39}

Once the ion heat flux-matched simulations were obtained, scans of $a / L T e$ enabled the measurement of the perturbative thermal diffusivity, calculated with Equation (8).

While particle transport was not the focus of this study, analysis was performed for both the L- and I-mode plasmas considered here. Following the approach outlined in Ref. 40, two trace (⁠ $0.0001 n e$⁠) impurity species (Z = 18, A = 40) were introduced into the ion heat flux-matched global simulations, and the impurity diffusion and convection were extracted. Similar levels of diffusion, convection, and peaking were found over the simulated radial region in both L- and I-mode, consistent with the measured impurity confinement times in L- and I-mode.¹ More detailed comparison between experimental particle transport levels and simulation will be addressed in the future work, and must take into consideration recent advances in the techniques for constraining experimental measurements of the diffusive and convective transport coefficients.⁴¹ 

The first comparison between experiment and GYRO is the electron and ion heat fluxes, characterizing the plasma power balance. The heat flux comparison is shown in Figure 5. The experimental heat fluxes in the figure are from TRANSP output and are represented in blue for L-mode and red for I-mode. The GYRO results are shown in black. Shaded regions represent uncertainties.

Figure 5 reveals that the GYRO simulations are able to match the ion heat flux within experimental uncertainties throughout the radial range of the global simulation in both L- and I-mode. In contrast, GYRO robustly underpredicts the electron heat flux at all radii in L-mode. In I-mode, GYRO underpredicts the electron heat flux between r/a of 0.65 and 0.8, but is able to match the experimental values within uncertainty between r/a of 0.8 and 0.88.

As an additional check of the accuracy of these global simulation heat fluxes, local nonlinear GYRO simulations were performed at various radii within the global radial range, and exhibit an excellent agreement (within 10%) with the global results. Results of local GYRO runs at r/a = 0.85 are shown in Figure 5.

The systematic underprediction of electron heat flux in Alcator C-Mod L-mode plasmas is well documented.^15,16 The only previous set of GYRO simulations of I-mode plasmas in Alcator C-Mod found a similar trend with the electron heat flux. A local simulation at r/a of 0.6 underpredicted the electron heat flux, but a local simulation at r/a = 0.8 agreed within the experimental uncertainty.¹² 

One additional source of electron heat flux is high-k ETG turbulence, which the linear analysis in Figure 4 shows to be unstable. For this reason, local nonlinear electron-scale GYRO simulations were performed for the L- and I-modes at r/a = 0.8. The electron-scale simulation box size was approximately 9 by 6 ρ_s in the radial and binormal directions with a grid spacing of $2 ρ e$⁠. These simulations captured a short wavelength turbulence (ETG) up to $k θ ρ s$ of approximately 56. Other inputs are the same as in Table II.

These simulations reveal that ETG drives a non-negligible amount of electron heat flux, but that it only accounts for approximately 10%–20% of the experimental heat flux, not nearly enough to resolve the discrepancy between the simulation and the experiment. Specifically, in L-mode, the experimental electron heat flux at $r / a = 0.8$ was approximately $0.31 MW / m 2$⁠. The ion-scale contribution, shown in Figure 5, was approximately $0.14 MW / m 2$⁠, while the electron-scale contribution to the electron heat flux was only $0.03 MW / m 2$⁠. Added linearly, the ion- and electron-scale simulations underpredict the experimental electron flux by more than 45% at this radius, which already has a closer agreement to experiment than most other radii. I-mode shows similar results.

As stated above, such a linear addition of single-scale simulations is generally not a valid approach to combining these two scales, as cross-scale coupling can lead to heat fluxes that differ significantly (both higher and lower) from a linear addition of two scales.¹⁵ The exact results are highly dependent on the specific plasma conditions and can only be assessed accurately with multi-scale gyrokinetic simulations.

In order to compare the experimental electron temperature fluctuation spectra to the output of GYRO, one must first apply a synthetic diagnostic to the simulation results. This synthetic diagnostic accounts for the finite spot size (assumed Gaussian with $L r = 1.2 cm$ and $L z = 0.64 cm$ (Ref. 24)) and k-sensitivity of the physical diagnostic, as well as the effects of plasma rotation. In addition, the synthetic diagnostic includes the effects of the physical diagnostic k-spectrum sensitivity. The specific implementation of the synthetic CECE diagnostic used in this study is described in Ref. 42.

The results of the fluctuation spectra comparison at $ρ t = 0.80$ are shown in Figure 6. GYRO underpredicts the electron temperature fluctuations fairly significantly in L-mode, predicting a level below the sensitivity limit, $T ̃ / T < 0.69 %$ (calculated by integrating over all frequencies), as opposed to the experimental level of $T ̃ / T = 0.83$⁠. This is in contrast to the past work with L-mode plasmas that showed that even when the electron heat flux was underpredicted, GYRO still agreed quite closely with the electron temperature fluctuation spectrum.⁴² 

The plasma rotation was varied within the experimentally measured uncertainty for the L-mode case in order to obtain the best possible match of the spectrum shape. As described above, this is done as part of the synthetic diagnostic. This variation was not, however, able to obtain any closer agreement with the fluctuation spectrum than what is shown in Figure 6.

In I-mode, GYRO also predicts a fluctuation level below the sensitivity limit. Since the experiment also measured a temperature fluctuation level below the sensitivity limit, the most that can be said about the I-mode is that the GYRO results are not inconsistent with the experiment.

This study will now compare the experimental and simulated perturbative electron thermal diffusivities. As described above, the perturbative diffusivity was measured experimentally with partial sawtooth heat pulses. In GYRO, $a / L T e$ was scanned about the experimental value, using the ion heat flux-matched simulations as the base case. To clarify, heat pulses were not simulated in GYRO. The experimental heat pulse analysis and the $a / L T e$ scans performed in GYRO measure the same physical quantity, $χ e p e r t$⁠, albeit in very different ways.^31,34,43,44

Figure 7 shows the results of the GYRO scans. The figure is plotted as Q_e vs. $n e ∇ T e$ in order to have the physical units of $χ e p e r t$⁠, which are $m 2 / s$⁠. The circled points are the ion heat flux-matched simulations, and the points to either side are the scans in $a / L T e$⁠. The slope of the solid line (a fit to the GYRO data) is the GYRO perturbative diffusivity $χ G Y R O p e r t$⁠. The slope of the dotted line is the experimental perturbative diffusivity $χ E x p p e r t$⁠. In other words, if the GYRO results agreed with the experiment, one would expect the GYRO scan points to lie along the dotted line, rather than the solid line. The figure makes clear that GYRO significantly underpredicts the perturbative diffusivity in both L- and I-modes. The comparison is shown more quantitatively in Table III.

The values of perturbative diffusivity reported in Figure 7 and Table III are radially averaged over the GYRO simulation domain. This was done in order to match most closely with the experimental measurements, which are radially averaged over approximately the same domain. This is inherent in the measurement of perturbative diffusivity via the propagation of heat pulses.³¹ 

As a brief summary of these simulation results, GYRO is able to match the ion heat flux in both L- and I-mode across the global simulation domain. GYRO underpredicts the electron heat flux in both L- and I-modes, though it can match the I-mode electron heat flux at the outer radii. In addition, GYRO also underpredicts the electron temperature fluctuation level, as measured with CECE, in L-mode. In I-mode, both the experimental measurement and the GYRO prediction are below the diagnostic sensitivity limit, so the most that can be said is that they do not disagree. Finally, GYRO underpredicts the perturbative electron thermal diffusivity in both L-mode and I-mode, though the underprediction is more severe in L-mode (approximately an order of magnitude compared to a factor of 4).

These discrepancies between the simulation and experiment may seem rather severe, but the past work with L-mode plasmas revealed that multi-scale physics (the coupling of electron- and ion-scale turbulence) was able to resolve the disagreement.^15,31 Multi-scale simulations of GYRO were able to simultaneously match Q_i, Q_e, and $χ e p e r t$ within the experimental uncertainties for an Alcator C-Mod L-mode plasma.¹⁵ Unfortunately, these multi-scale GYRO runs are incredibly computationally expensive, so this study will investigate the effects of multi-scale physics using the reduced model TGLF.^14,19

TGLF is a reduced transport model, originally described in Ref. 14, but has recently been upgraded to include the cross-scale coupling found in multi-scale gyrokinetic simulations.¹⁹ As stated above, multi-scale GYRO is computationally expensive, requiring approximately 100 million CPU hours for a validation study,¹⁵ while the multi-scale TGLF runs on the order of minutes.

The TGLF simulations in this study were run with the same experimental inputs as the GYRO simulations, but were run locally, at r/a = 0.64, 0.74, and 0.84. In order to be consistent with the GYRO simulations, only $a / L T i$ was varied within the experimental uncertainty in order to match the experimental ion heat flux. These simulations included the heat flux contributions from $k θ ρ s$ up to 24, and all used the SAT-1 version of TGLF, which includes the multi-scale effects described above.¹⁹ 

As with GYRO, this section will compare the experimental Q_i, Q_e, and $χ e p e r t$ to TGLF simulations.

The TGLF heat fluxes are compared with the experimental heat fluxes in Figure 8. Like GYRO, TGLF is able to match the experimental ion heat flux at all radii in L-mode. In I-mode, the two inner radii are well-matched, and there is only a slight over-prediction at the outermost radius.

Also like GYRO, TGLF underpredicts the electron heat flux everywhere in L-mode and everywhere except the outermost radius in I-mode. The TGLF electron heat fluxes are closer to experiment than those predicted by GYRO, but they still fail to match within the experimental uncertainty.

This underprediction of electron heat flux with multi-scale TGLF SAT-1 contrasts with previous results, where multi-scale TGLF was able to reproduce experimental heat fluxes in an Alcator C-Mod L-mode plasma.¹⁹ We note that the conditions and radial locations studied here differ from the previous investigations. Possible explanations for this discrepancy will be discussed in Sec. V.

In exactly the same manner as with GYRO, the perturbative diffusivity was extracted from TGLF by varying $a / L T e$ around the ion heat flux matched simulations. This was done at all three radial locations, and then the three perturbative diffusivity values were averaged together for comparison with the experiment. The results of these scans are shown in Table III.

As with GYRO, TGLF severely underpredicts the perturbative diffusivity in both L- and I-modes. While the TGLF prediction in L-mode is slightly closer to the experiment than GYRO (75% higher than GYRO), the I-mode prediction is actually further away. This result indicates that even multi-scale TGLF SAT-1 is not properly capturing the plasma stiffness in either L- or I-mode.

As with the ion-scale GYRO, multi-scale TGLF simulations varying only $a / L T i$ were able to match the experimental values of Q_i (with the slight exception of the outermost radius in I-mode), but underpredicted Q_e and perturbative diffusivity in both L- and I-modes. This result was somewhat unexpected, since the past work with GYRO simulations revealed that the inclusion of multi-scale effects was able to resolve the discrepancy in both Q_e and the perturbative thermal diffusivity in a L-mode plasma.¹⁵ Multi-scale TGLF simulating that same plasma was able to reproduce the GYRO results.¹⁹ 

The results of the validation study presented here raise important questions about the role of cross-scale coupling and multi-scale turbulence in high performance plasmas, as well as about current modeling capabilities.

The first of these questions asks why multi-scale TGLF simulations do not resolve the discrepancy in electron heat flux and perturbative thermal diffusivity between ion-scale simulations and experimental measurements.

There are a number of possible explanations for this discrepancy, some of which are already under investigation as part of the future work. First, it is possible that the multi-scale TGLF used in this study does not fully capture the multi-scale effects in these plasmas. In other words, multi-scale TGLF was tuned to the results of multi-scale GYRO for a specific L-mode plasma and the immediately surrounding parameter space (see Ref. 19 for more details on the tuning procedure). It may be that TGLF is well tuned to multi-scale GYRO in the parameter range in Ref. 15, and so was able to reproduce the multi-scale GYRO's results for that plasma, but does not fully capture the multi-scale effects in the parameter space of the plasmas used in this study. Unfortunately, investigation of this possible explanation for the disagreement between multi-scale TGLF and experiment would require multi-scale GYRO runs, which, as was stated above, are enormously demanding in terms of computational resources. If this is the case, then this study may help identify which regimes of plasma operation require multi-scale GYRO simulations, and which can be sufficiently modeled using multi-scale TGLF.

Another possible explanation for the discrepancy is that one must simultaneously vary multiple TGLF (or GYRO) input parameters, rather than just $a / L T i$⁠, in order to correctly reproduce the experimental results. Other possible inputs that may vary are $a / L T e , a / L n e$⁠, Z_eff, $T e / T i$⁠, E × B shear, v_tor, $ν *$⁠, etc., which each have their own experimental uncertainty and are known to play important roles in determining the local stability of both ITG and TEM turbulence. These values were not varied in this study in order to be consistent with GYRO, where only $a / L T i$ varied (multi-parameter scans are possible in GYRO, but require significantly more computation resources). In general, however, it may be that only the multi-parameter scans can correctly reproduce experimental results. One must be careful that these scans are done self-consistently (i.e., one cannot change Z_eff without also changing collisionality). This work is already in progress, and the initial results are promising.

Finally, it may be that the additional effects must be included in the GYRO or TGLF simulations in order to properly capture all of the relevant physics. For example, several recent experiments have shown that fast ions may have a significant impact on tokamak microturbulence, and must be included in gyrokinetic simulations to properly reproduce the experimental transport.^33,45 It may be that the simulations described in this study require the effects of fast ions to match the experimental measurements. The inclusion of fast ions in gyrokinetic simulations for the L- and I-mode described here is now in progress.

In addition to exploring the possible reasons for the discrepancy between multi-scale TGLF and experimental measurements described above, another avenue of exploring the transport in I-mode involves expansion of this study to additional machines. To this end, a new CECE diagnostic has recently been installed on the ASDEX Upgrade tokamak.⁴⁶ This diagnostic will also be able to operate as an “n-T phase” diagnostic: CECE coupled with a Doppler reflectometer that can measure the relative phase between density and temperature fluctuations at the same location in the plasma.^47,48 Future work will involve the study of L- and I-mode plasmas on ASDEX Upgrade and comparing experimental results to gyrokinetic simulations with GENE. Preliminary results with this diagnostic have shown a good agreement between experimental heat flux and fluctuation measurements and nonlinear, ion-scale GENE simulations, which are intriguing given the disagreements on Alcator C-Mod presented in this paper. Similarly, it would also be possible to run GENE on the Alcator C-Mod plasmas described in this paper.

Collecting these observations, future validation work will involve additional multi-scale GYRO runs for tuning of TGLF, self-consistent variation of multiple input parameters with TGLF (and possibly also ion-scale GYRO), inclusion of additional physics in GYRO, and expansion of this validation effort to ASDEX Upgrade (including comparison of power balance, fluctuations, and perturbative diffusivity).

In addition to motivating future work concerning the specific validation results presented here, this study also motivates further discussion and investigation of I-mode in general. Experimentally (now that C-Mod has completed operation), it would be worthwhile to continue to study the properties and operation of I-mode on ASDEX Upgrade and DIII-D, on which I-mode has already been accessed. In addition, expansion of I-mode to additional machines, such as TCV, JET, EAST, KSTAR, and WEST, would significantly improve our understanding of its operation, and would enable a more accurate extrapolation to a future reactor. Work on these machines may shed light on the fundamental differences between L-, I-, and H-modes, and thereby deepen our fundamental understanding of plasma turbulence.

To summarize how the results presented in this paper fit into the bigger picture, past ion-scale GYRO simulations of Alcator C-Mod L-mode plasmas match the experimental ion heat flux, but robustly underpredict the electron heat flux and perturbative diffusivity. Past work was able to resolve these discrepancies in L-mode, however, with multi-scale GYRO simulations. The new validation work presented in this paper, using global, nonlinear, ion-scale GYRO simulations, also underpredicted the electron heat flux and perturbative diffusivity in both L- and I-modes, but was surprisingly unable to resolve these discrepancies by including the multi-scale effects with TGLF. These results motivate further work with multi-scale GYRO simulations, a more thorough investigation with TGLF, and multi-machine comparisons of I-mode validation.

This work was supported by the U.S. DOE under Grant Nos. DE-SC0006419 and DE-FC02-99ER54512-CMOD, and by the U.S. DOD under the NDSEG Fellowship. GYRO simulations were performed with MP224 Repo at the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. TGLF simulations used resources of the MIT PSFC parallel AMD Opteron/Infiniband cluster Loki.