Simultaneous reproduction of experimental profiles, fluxes, transport coefficients, and turbulence characteristics via nonlinear gyrokinetic profile predictions in a DIII-D ITER similar shape plasma

The worldwide push for the development of a fusion pilot plant (FPP) and recent large investments in fusion via private industry has accelerated the demand for accurate modeling of fusion devices. Modeling and prediction of future burning plasmas has the potential to optimize device design and performance, potentially decreasing the overall capital cost of next-gen fusion devices. Transport in the core of tokamaks is now generally accepted to be the result of plasma turbulence, which is driven unstable by the free energy present in gradients of the plasma density and temperature and ultimately limits the overall performance of tokamak fusion reactors. Turbulence in tokamak plasmas is often characterized by its poloidal wavenumber (⁠ $k θ$⁠) multiplied by the sound speed gyroradius [ρ_s, the Larmor radius evaluated with the ion sound speed (c_s)]. Ion-scale turbulence, also referred to as low-k turbulence exists with $k θ ρ s < 1.0$⁠, while electron-scale turbulence or high-k exists with $k θ ρ s > 1.0$⁠. Ion temperature Gradient (ITG) and Trapped Electron Mode (TEM) turbulence often playing dominant roles at the ion-scale and high-k TEM and Electron Temperature Gradient (ETG) turbulence playing important roles at the electron scale. Years of theoretical and computational work led to the development of the gyrokinetic model and its implementation in a wide range of codes that are capable of solving the gyrokinetic equation, typically using high performance computing. Since the development of gyrokinetic codes, a worldwide effort has been underway to validate the gyrokinetic model against experiment.^1–15 Most of these efforts have focused on comparisons of simulation with single transport channels. Comparisons with turbulent fluctuations have been comparatively limited due to limited number of devices fielding the relevant diagnostics and the complexity associated with synthetic modeling of gyrokinetic fluctuations. In this work, we present a multi-channel, multi-field validation of the nonlinear gyrokinetic model in a DIII-D ITER similar shape (ISS) discharge. Comparisons with experiment span a large range of experimentally inferred values, measurements, and radial locations including: ion and electron heat fluxes, electron and impurity (Li, C, and Ca) fluxes and transport coefficients, and turbulent characteristics (i.e., low-k $δ n e / n e$ and $δ T e / T e$ fluctuation levels, fluctuation level trends with radius, spectral shapes, and cross power), spanning from approximately ρ (square root of normalized toroidal flux) 0.3–0.8. This validation effort is therefore unique in that spans multiple levels of the primacy hierarchy^16,17 as it compares gyrokinetic modeling with turbulence driven fluctuations in background profiles, inferred heat and particle fluxes, which themselves result from the fluctuations, as well as measured density and temperature profiles that result from a combination of sources, sinks, and transport.

The remainder of this paper is organized as follows. Section II covers the experimental setup of the DIII-D ISS discharge, the simulation setup used in this work, and a description of the surrogate technique used for profile prediction. Section III describes the flux-matched profiles obtained via nonlinear gyrokinetic simulation. Section IV describes comparisons of the simulations results with turbulence measurements and other fluctuation characteristics. Section V compares the gyrokinetic results with impurity transport, including an investigation of multi-Z impurity transport in this condition. Finally, Sec. VI summarizes the findings and includes some discussion of the results.

The discharge studied in this work was performed on the DIII-D tokamak,¹⁸ and the experiments analyzed in this work were motivated by the FY20 US Joint Research Target (JRT) focused on studying core to edge impurity transport. The conditions were chosen because they allowed for excellent diagnostic coverage across the profile in discharge conditions that were relevant for future ITER operation. More specifically, the discharge of interest was an Resonant Magnetic Perturbation (RMP) edge-localized mode (ELM)-suppressed H-mode performed with $B T = 2.0 T , I p = 1.6 MA , q 95 = 3.45$⁠, line averaged $n e ∼ 3.0 × 10 19 m − 3$⁠, and $β N ∼ 1.35$⁠, with mixed neutral beam injection (⁠ $P NBI ∼ 4.0$ MW) and electron cyclotron heating (ECH) (P_ECRH = 1.6 MW) deposited around $ρ = 0.25$⁠. Time traces from key experimental quantities are plotted in Fig. 1. In addition to utilizing a similar shape and q₉₅ to the projected ITER baseline scenario, the use of both counter and co-current beams was employed to reduce the total input torque and obtain a more ITER-relevant rotation. It is also notable that the heating mix in this condition resulted in experimental heat fluxes that are similar to those expected in the ITER baseline condition, with an ion to electron heat flux ratio $( Q i / Q e )$ that generally exceeds 1.0. This aspect is particularly important to the relevance to ITER baseline operation as it implies similar turbulence will be present in both devices. More description of ITER flux ratios and can be found here.¹⁹ 

High resolution profiles, turbulence data, and impurity transport data were obtained through repeat discharges of the target shot 183 185. The collection of validation quality datasets in this discharge represented a key component that motivated further investigations presented in this work. Carbon ion temperature, carbon density, and carbon rotation were measured using Charge Exchange Spectroscopy (CER).²⁰ Electron temperature was provided by Thomson scattering,²¹ and electron density was primarily obtained via the use of a density profile reflectometer²² and radiated power was measured with DIII-D's bolometer arrays. The combined high quality profile data allowed for accurate power balance calculation with the TRANSP code²³ which provides the experimentally inferred heat and particle fluxes utilized for the profile prediction described later in this work. A wide range of impurities were introduced into the studied conditions ranging from He to W (Z = 2–74). He was introduced via gas puffing, Li was introduced via the use of the impurity powder dropper,²⁴ C is intrinsic in DIII-D, and F, Al, Ca, and W were all introduced via the DIII-D laser blow-off (LBO) system.²⁵ All impurities other than W were measured via CER resulting in profiles of impurity density in time, with W being primarily measured with soft x-ray (SXR) arrays and Ca measured by both CER and SXR measurements. Although all of these impurities were introduced, the data quality generally limited our investigations in this work to Li, C, and Ca. However, we note that still represents a relatively large span of impurity Z and that additional impurity investigations may be the part of future work. All experimental measurements were used as inputs to TRANSP for calculation of the inferred experimental heat and particle fluxes that were used for comparison with the gyrokinetic modeling. Analysis focused on the time range from 3600 to 4000 ms in discharge 183 185 where MHD activity is relatively low and the discharge is steady. Infrequent sawteeth are present in the discharge, but the input profiles that are utilized for the remainder of this work are obtained from an average of the above time range, and the experimental heat fluxes were determined via TRANSP modeling performed using the OMFIT framework.^26,27

All turbulence simulations performed in this paper utilized the continuum gyrokinetic code, CGYRO.²⁸ CGYRO is a mature gyrokinetic code that has been optimized for modern computing architectures [such as hybrid central processing unit (CPU)/graphics processing unit (GPU) systems] and is optimized for efficient calculation of fusion relevant plasmas, including collisional, near-edge conditions, and multi-scale plasma turbulence. All of the simulations run in this paper were performed on the CORI-Knight's Landing (KNL) partition at the NERSC computing facility. Simulations performed here were local simulations utilizing periodic boundary conditions, allowing us to avoid the potential impact of boundary condition choice on results as shown in recent global simulations.²⁹ Simulations utilized all experimental inputs from DIII-D discharge 183 185 including measured profiles (density, temperature, etc.), rotation, realistic geometry,³⁰ three gyrokinetic species (D, C, and e-), and included high physics fidelity including: E × B and rotation effects, electromagnetic fluctuations (⁠ $δ ϕ$ and $δ A | |$⁠), and the Sugama collision operator.³¹ Each simulation domain was represented by 20 toroidal modes (n_n) and approximately 350 radial modes (n_x) with a radial and binormal box size of $[ L x , L y ] ∼ [ 100 × 100 ρ s ]$ which captured turbulence up to $k θ ρ s ∼ 1.2$⁠, using n_energy = 8 (energies), $n θ$ = 24 (points in theta), and $n ξ$ = 24 (pitch angles). We note that CGYRO uses the effective field strength $B unit = 1 / r × d χ t / d r$⁠, where χ_t is the toroidal flux divided by 2π when evaluating the value of ρ_s throughout this work. These resolutions are sufficient for resolving long wavelength turbulence such as microtearing modes (MTM), ion temperature gradient (ITG) driven modes, and trapped electron modes (TEM) that are generally found to result in significant heat and particle loss in current tokamaks. The objective of the simulations in this work was to iterate to simultaneously obtain ion and electron heat flux and electron particle flux-matched profiles. In this context, flux-matched profiles refer to profiles where the heat and particle fluxes obtained via nonlinear gyrokinetic simulation match the target (experimental) fluxes. This procedure is described in more detail in Sec. II C, but it should be noted that this process required approximately 64 total nonlinear gyrokinetic simulations for completion (4 radial locations × 16 iterations). Therefore, there was significant desire to reduce the overall computational time used. Although simulations were typically run for approximately 500–600  $a / c s$⁠, time averaging windows utilized in this work were typically ∼300  $a / c s$ to obtain time averaged quantities such as heat and particle fluxes. It was found that these averages were sufficiently long for the conditions simulated to provide accurate time averages without wasting computing resources.

It is notable that the simulations performed in this paper did not capture contributions from high-k, electron-scale turbulence such as high-k TEM and ETG turbulence. Accurate modeling of conditions with significant low and high-k turbulence unstable would require the use of multi-scale gyrokinetic simulations,^32–34 capable of simultaneously simulating both turbulent scales and their interactions, and is, in general, computationally intractable. More importantly, the ion to electron heat flux ratios at the radial locations studied in this plasmas point a negligible or subdominant role of electron-scale turbulence. As discussed in Refs. 35, 36, and 19, plasma conditions exhibiting $Q i / Q e ≳ 1.0$ should be dominated by long wavelength modes such as ITG, TEM, and MTM as described by the transport “fingerprints” argument laid out by Kotschenreuther and colleagues.³⁷ This conclusion is supported by the linear stability results reported in Fig. 2. Although electron-scale turbulence is found to be unstable at outer radii, its linear growth rates compared the ion-scale turbulence growth rates are small relative to the low-k [as quantified by max (⁠ $γ high − k$⁠)/max(⁠ $γ low − k ) ≪ 60.0$]. This criteria^38,39 and similar ones⁴⁰ have been identified as indicating when electron-scale turbulence, and the associated cross-scale coupling should be negligible. Therefore, we proceed with analysis in this work focusing on ion-scale turbulence.

The framework used for the profile predictions in this work was described recently in Ref. 36. Details of the profile prediction method, including numerical techniques, and considerations for burning plasma prediction, can be found in Ref. 41 and are out of the scope this paper. However, the basic approach used is briefly described here for completeness. A total of four radial locations in the plasma were chosen for simulation that allowed for an accurate prediction of the inductive discharges studied here: $ρ =$ 0.3, 0.43, 0.62, and 0.8. The innermost radii was limited due to proximity to the q = 1 surface. The outermost radii was limited to $ρ = 0.8$ with the density and temperatures fixed to their experimental values. This location was chosen to (1) avoid effects of the RMPs in the edge which have been shown to effect near-edge transport^42,43 and (2) to avoid investigation of the pedestal region which is out of the scope of this paper. Due to the computational expense of the calculations, the sensitivity of the profile predictions to boundary condition location and value was not investigated in this work. Profiles inside of the innermost radial location simulated are performed by linearly interpolating the normalized scale lengths to zero on axis. Therefore, $a / L x = 0$ on axis, where x = n_e, T_e, or T_i.

To accelerate the convergence to steady-state, flux-matched conditions we use the new PORTALS framework,³⁶ as opposed to standard Newton or gradient-based methods. PORTALS has demonstrated that self-consistent, multi-channel, flux-matched profiles can be obtained with a reduced number of transport model evaluations, without any loss of accuracy (as the final solution is flux-matched on the actual transport model, even if surrogates are used to find such conditions). PORTALS utilizes the auxiliary heating profiles calculated by TRANSP and calculates the radiation and collisional exchange self-consistently for the profile predictions at each iteration. The nonlinear gyrokinetic code CGYRO was used to evaluate the turbulent contributions to Q_i, Q_e, and Γ_e, while the neoclassical code, NEO,⁴⁴ was utilized to calculate the corresponding neoclassical contributions. These tools are then used to solve the set of time-independent, coupled heat and particle transport equations. The initial magnetic geometry utilized was generated via a kinetic EFIT (equilibrium fitting) and is approximated as fixed during the iteration procedure. In principle, a self-consistent update of the equilibrium could be performed at each iteration but the small profile variations would likely result in only small changes to the equilibrium and such a procedure is left for future work. Rotation effects were included in each simulation but the rotation profile itself was set to the experimentally measured value and not evolved during the profile predictions. Finally, although carbon was included in the gyrokinetic simulations, the profiles of this impurity were not evolved during the convergence process, and the experimental concentration (⁠ $n c / n e$⁠) was held constant during electron density profile predictions.

To start the process, the experimental profile + 4 initial profiles based on random variations of input gradients (⁠ $a / L T i , a / L T e , a / L n$⁠) were simulated using nonlinear gyrokinetics (at four radial locations). This initial set of simulations is used to train surrogate models for each of the transport fluxes, thus reproducing the nonlinear gyrokinetic dependencies at each radial location. In this context, a surrogate model is a statistical model, capable of rapid evaluation, that is meant to accurately approximate the results from the nonlinear gyrokinetic simulations; for example, to emulate the response of heat and particle fluxes to changes in normalized gradient scale lengths. The surrogate model is then used to attempt to identify the set of gradients $a / L T e , a / L T i$⁠, and $a / L n e$ that result in matched (i.e., transport and targets) Q_i, Q_e, and Γ_e at each radial location. The predicted solution is then tested with fully nonlinear gyrokinetic simulations. If the nonlinear gyrokinetic simulations are flux matched, the process is complete. If not, the data from the high fidelity simulations are reported back to further train the surrogates, and this process iterates until convergence is achieved. The code iterates trying to find a solution where the profiles simultaneously match the ion and electron heat transport and electron particle transport by varying the normalized gradient scale lengths thus predicting the T_e, T_i, and n_e profiles. The benefit of this approach is a more rapid convergence to a converged solution, which exceeds those of more traditionally used Newton solvers by a factor of 4–6. As a result, these three channel predictions, with three gyrokinetic species and relatively high physics fidelity utilized only $∼ 1.5$ M CPU hours on the Knight's Landing (KNL) partition of the NERSC Cori supercomputer. This amount of computing time is comparable to what might typically be used to investigate a single radial location with nonlinear gyrokinetics by performing independent, single parameter scans.

The objective of this work was to perform profile predictions of the target DIII-D discharge with high physics fidelity and to compare the predicted profiles with those measured in experiment. As described above, the prediction of the flux-matched profiles required 16 total iterations with the first 5 being purely random variations utilized to train the surrogates to begin the optimization.

To initialize these random variations, gradient scale lengths of T_e and T_i were allowed to vary by 50% at each radial location relative to their experimental value, while values of $a / L n e$ were allowed to vary by 100%. Figure 3 plots the different profiles (blue) that were attempted by the surrogate model iterations to achieve the flux-matched conditions. The initial profile used for the convergence process is indicated in red with the final converged value is represented by green.

The original and final values of the fluxes (Q_i, Q_e, and convective heat flux $Q conv = 5 / 2 T e Γ e$⁠) are plotted in Fig. 4. As seen in this figure, despite a large disagreement between the target and simulated fluxes for the initial profiles, generally good convergence was achieved in all channels in the final reported profiles. We note that the initial and final target fluxes were very similar in this condition because the auxiliary heating profiles were dominant (relative to the exchange, Ohmic heating, and radiation) and assumed unchanged during the iterations. However, this situation will no longer exist in a burning plasma where the heating profile is dominated by alpha power. We also note that no attempt to match the Q_e at the $ρ = 0.3$ location was made since (1) this location is very close to ECH deposition location (centered around $ρ ∼ 0.25 )$ and therefore has large uncertainty in the electron heat flux and (2) it is located inside the sawtooth inversion radius. Therefore, the electron fluxes at this location have higher uncertainty and are likely effected by periodic sawtooth activity that would not be captured in a gyrokinetic simulation.

The approach detailed above resulted in predicted profiles consistent with heat and particle fluxes obtained from power balance. In Fig. 5, the predicted profiles of T_e, T_i, and n_e are compared directly with the experimental data obtained in the analyzed time window spanning t = 3600–4000 ms. Electron temperature measurements were obtained via Thomson scattering, ion temperature was measured using charge exchange from the intrinsic C impurities, and the density profile data are the result of density profile reflectometer measurement. As seen in Fig. 5, there is generally excellent agreement between the simulated an experimental profiles within the scatter of the data obtained from the experiment. Agreement is, however, not perfect with some slight disagreement in the T_i profile at inner radii and some slight differences in the n_e profile. Despite slight difference, the predicted profiles clearly lie within the experimental scatter. It is also worth noting that the experimental stored energy and energy confinement are reproduced by the modeling within 5%, well within the typical scatter obtained from empirical scaling such as $τ 98 , y 2$⁠.⁴⁵ This level of agreement in zero-dimensional quantities, coupled with the good reproduction of density and temperature profiles, suggests that these profile prediction techniques may serve as a viable alternative to empirical scalings for the design of new fusion devices. Comparisons of nonlinear gyrokinetic predicted profiles with experiment are quite rare,^46,47 as typically nonlinear gyrokinetics is used only to investigate profiles a discrete radial locations with no consideration of changes in the profiles that might occur from changes to the local gradients and changes to the target powers (e.g., via self-consistent energy exchange, radiation, alpha heating, etc). For that reason, this technique has the significant power of both being flux-matched locally in multiple channels but also predicting profiles that allow for comparison with density and temperature profile data directly. Overall, the good agreement with experiment provides confidence in the validity of the gyrokinetic model in accurately describing the ion and electron heat fluxes and electron particle fluxes in the studied experimental plasma conditions and demonstrates the power of the surrogate accelerated profile predictions.

Although validation of gyrokinetics against experimentally inferred heat fluxes is relatively common, direct comparison with measured turbulence is comparatively less common. Turbulent fluxes such as Q_i, Q_e, and Γ_e are themselves the result of the turbulent fluctuation levels and cross-phases between these fluctuations and therefore are at a higher level in the so-called primacy hierarchy.^16,17 Thus, comparisons with fluxes represent a less stringent test of the gyrokinetic model compared with direct comparison with fluctuation properties. In order to ensure that gyrokinetics predicts the correct fluxes, for the right reasons, comparison with more fundamental turbulence characteristics such as the fluctuation amplitude, fluctuation spectra, and cross-phases are highly desired.

The extensive diagnostic suite on DIII-D enabled the measurement of both low-k (ion-scale) fluctuations in density and electron temperature via the use of the Beam Emission Spectroscopy (BES; $δ n e / n e$⁠)⁴⁸ and Correlation Electron Cyclotron Emission (CECE; $δ T e / T e$⁠)⁴⁹ diagnostics, respectively. These measurements are capable of providing radial profiles of fluctuation amplitudes and the frequency spectra between two correlated measurements. Comparison of gyrokinetic simulation with these measurements requires the use of a synthetic diagnostic such as those described in Ref. 5. For this work, the synthetic diagnostic were ported over from their original formulation with the GYRO code² to the CGYRO code²⁸ utilized throughout this work. Accurate comparison of simulated fluctuations with with experiment via synthetic diagnostics required the output of fluctuation data on a more dense poloidal grid (24 points) such that the R, Z fluctuations (where R is major radius and Z is the vertical coordinate) could be more accurately reconstructed. Additionally, the simulations were run for extended time periods (⁠ $∼ 1000 a / c s$⁠) to enhance statistics. These R, Z fluctuations are then convolved with a filter function representing the measurement volume of the diagnostic to mimic its response. This technique enabled direct comparison of measurements with fluctuation levels, radial trends, frequency spectra, and cross-phases as described in this section. For more details on the synthetic diagnostics for BES and CECE, the reader is referred to Ref. 5 for more details.

BES measurements were obtained between $ρ ∼ 0.6$ and $ρ ∼ 0.8$ that allowed for comparison with the outer two radial locations simulated (⁠ $ρ = 0.62$ and 0.8). Unfortunately, data were not available at the inner radii (0.3 and 0.43) but likely would have resulted in fluctuation levels that were below the sensitivity limit of the diagnostic due to the relatively high confinement of the conditions studied and generally low fluctuation levels found in H-mode plasmas. Figure 6 plots the radial profile of BES measurements obtained from experiment compared with those calculated via synthetic diagnostic modeling at $ρ = 0.62$ and 0.8. As seen in the figure, a general increasing trend of fluctuation level moving to the edge is found. This radial trend is captured by the gyrokinetic modeling with a clear increase from 0.24% to 0.35% in $δ n e / n e$ found at $ρ = 0.62$ and 0.8, respectively. Most importantly, the quantitative comparison with the fluctuation levels at both radial locations is in excellent agreement with experimental measurements.

Utilizing two poloidally spaced channels of the BES diagnostic, one is able to extract density fluctuation cross-power and cross-phase information at different radial locations. In Fig. 7, we plot a comparison of the BES cross power spectra at both $ρ = 0.62$ and 0.8 with synthetic gyrokinetic modeling from 70 to 450 kHz. The minimum of 70 kHz was chosen to eliminate any potential contributions from low frequency MHD that might be present. Excellent agreement is once again found when comparing the cross power spectra. At $ρ = 0.62$⁠, there is some difference between the measured and synthetic fluctuation spectra, with the experimental measurements exhibiting a broader spectrum than found in the simulation. However, this difference may be attributed to (1) the slight differences in the measured and simulated radial location (sim: $ρ = 0.62$⁠, exp: $ρ ∼ 0.6$⁠), (2) an artifact of the relatively short time over which the synthetic diagnostic averages (compared to the experiment which averages over multiple confinement times, or (3) explained by relatively uncertainties in the measured toroidal velocity at this location, which can alter the amount of Doppler broadening of the fluctuations. At $ρ = 0.8$⁠, excellent agreement is found across the cross-power spectrum with the peak well reproduced along with the shape of the spectrum by the synthetic modeling. Despite some misalignment in the spectrum peak at $ρ = 0.62$⁠, the fluctuation levels at both radial locations are in excellent agreement with experiment. Peaks are observed in the synthetic spectrum in both the BES (Fig. 7) and, as will shown later, the CECE (Fig. 9). As described in Ref. 5, this structure results from the finite toroidal mode numbers captured in the local, flux-tube simulations and has been shown previously not to effect the overall fluctuation levels.

Figure 8 plots the cross phase vs frequency between two poloidally spaced BES channels for both experimental data (blue) and synthetic data (black) at $ρ = 0.62$ and 0.8. The slope of the cross-phase with frequency is set by both the local rotation and the phase velocity of the turbulence, with the former representing typically the dominant contribution. In the frequency range of interest, 70–450 kHz, the synthetic BES spectrum is in excellent agreement with the measured values. Although the cross-phase between two poloidally separated channels is dominated by the local rotation, if the turbulence simulated via nonlinear gyrokinetics was largely inconsistent with that measured in experiment, one might expect a clear deviation of the synthetic data from the measurement. Therefore, although not a strong test of the model validity, the good agreement between the simulation and experiment suggests that the phase velocity of the simulated turbulence is likely in good agreement with experiment, which provides an additional validation of the modeling against experiment.

The CECE system on DIII-D has recently been upgraded to allow for the measurement of radial profiles of low-k $δ T e / T e$ fluctuations during a single shot.⁴⁹ The data obtained during DIII-D shot 183 185 utilized 16 total CECE channels for eight radial measurements spanning from $∼ ρ = 0.53 − 0.93$⁠. Radial measurements at $ρ = 0.63$ and 0.78 were compared with the nonlinear gyrokinetic simulated results at $ρ = 0.62$ and 0.8. At these radial locations, the poloidal and radial spot sizes (cm) of the CECE measurements were 4.22 and 4.36 cm in the poloidal direction and 1.62 and 1.54 cm in the radial direction. The poloidal and radial spot sizes are utilized as part of the synthetic modeling to enable comparison between experiment and simulation. Analysis of the experimental CECE data was performed using the techniques presented in Ref. 50.

The comparison of simulated and synthetic CECE fluctuations is plotted in Fig. 9 at the radial location of $ρ = 0.8$⁠. Similar to the BES, comparisons are made only over the 70–500 kHz range with lower frequencies being omitted from the comparison to avoid any possible contributions from low frequency MHD.

As seen in this figure, the simulation appears to underestimate the total fluctuation level measured via experiment. At $ρ = 0.8$⁠, the measured $δ T e / T e$ is found to be 1.62% compared with 0.82% found via synthetic modeling. A similar result is found at $ρ ∼ 0.62$ where the experimental level is, again approximately 2× higher than the synthetic value (exp: 1.18%, syn: 0.60%). We note that for the experimental CECE data in Fig. 9, the noise level has been calculated. This value is determined by randomizing the time histories of measured signals and then calculating the power spectrum of the measured signals. For these data, this value is estimated to be approximately 0.7%–0.75% at the locations studied. We note that multiple approaches are used for the analysis of CECE data in the literature,^50,51 and differences in these approaches may effect quantitative calculation of fluctuation levels. Despite the observed disagreement in the fluctuation level, we note that the radial trend of increasing fluctuations with radius is reproduced by the CGYRO simulation (⁠ $ρ = 0.62$⁠; sim: 0.60%, exp: 1.18%) (⁠ $ρ = 0.8$⁠; sim: 0.82%, exp: 1.62%). Additionally, the shape of the spectra is found to be in good agreement, as demonstrated in Fig. 9 for $ρ = 0.8$⁠. In this figure, the fluctuation data has been scaled such that the simulated and experimental data effectively overlap. The purpose of this comparison is to emphasize the similar spectral shapes. Although not plotted, qualitatively good agreement is also found at $ρ = 0.62$⁠, but the spectrum is broader and closer to a noise level at this location make the comparison less definitive. The agreement of nonlinear simulations with the experimental heat and particles fluxes while showing some disagreement with the $δ T e / T e$ fluctuation amplitude may on the surface seem inconsistent. However, it is important to note that the fluxes themselves are set through a combination of the fluctuation levels between multiple quantities (such as $δ ϕ , δ T e , and δ n e$⁠) and the cross phases between these fluctuations. Therefore, it is possible that agreement could occur in fluxes but not in the turbulent fluctuation levels. However, the addition of multiple channels in this validation exercise makes that type of occurrence presumably less likely and could potential point to other sources of error such as systematic uncertainties.

The results from Secs. II–IV demonstrated comparison of gyrokinetic simulation results with ion and electron heat fluxes, electron particle fluxes, and low-k density and temperature fluctuation measurements. Generally, good agreement was observed when simulation was compared with experiment. In this section, we extend our investigations to that of impurity transport spanning low to medium-Z impurities. These investigations are important as it is expected that impurities will play a more important role in next generation devices that operate in burning plasma regimes. In this section, we will compare simulated impurity transport with experimental impurity transport of Li, C, and Ca to span a wide range of Z. More specifically, measured Li peaking profiles will be compared with simulation, measured density profiles of fully stripped carbon (n_c) will be compared with simulation, and D and V profiles obtained from laser blow-off injected Ca impurities will be compared with simulation. Also, a more theoretical study into multi-Z impurity transport is presented that suggests the existence of a Z dependence of impurity diffusion in the discharge studied and provides some insight into its physical origin.

Repeat discharges of the target condition were used to study a range of impurity Z in the ITER similar shape condition. At low-Z (Z = 3), the DIII-D impurity powder dropper²⁴ was used to drop small granules of Li into the plasma and the emission from these granules was measured via charge exchange spectroscopy. Analysis of these impurities was focused on time periods well after the initial injection of the granules so that the impurity profile is essentially stationary and decaying self-consistently in time. During this time window, the impurity transport code Aurora⁵² was utilized to iterate over a 1.5D model of impurity transport to determine the impurity diffusion and convection (D and V) that produced the best match of the experimental data. This was performed under the assumption that the impurity flux can be written in the form: $Γ = − D ∇ n + V n$⁠, with D and V transport coefficients independent of charge state. Experimental analysis utilized Bayesian inference techniques to sample different combinations of D and V profiles and synthetic signals (in this case, Li density profiles) are generated at each iteration. This process is continued until the difference between measured Li signal and those predicted by Aurora is minimized. For more details in the analysis technique, the reader is referred to Refs. 53 and 54. The flux-matched profile simulations that were discussed in Sec. II were then used to evaluate the impurity transport coefficients utilizing the method that is described in Ref. 55. Essentially, two trace impurity species (⁠ $10 − 6 × n e$⁠) are introduced into the gyrokinetic simulation with identical charge and mass but different values of impurity density gradient scale length (⁠ $a / L n z$⁠). For the purpose of this work, values of 0.5 and 2.0 were used. Trace impurity temperatures are set to be equal to the deuterium temperature. The resulting impurity species fluxes normalized to their density (⁠ $Γ / n z$⁠) are fit to − $∇ n z / n z$ to obtain the trace impurity diffusion and convection coefficients from the linear fit. The results from the impurity peaking (−V/D) obtained from this analysis, compared with the experimental values obtained from the Li profile (⁠ $− V / D = 1 / L n L i$ in steady state), are plotted in Fig. 10.

As seen in Fig. 10 for the low-Z Li impurity studied, there is reasonably good agreement between simulation and experiment within uncertainties at 0.3, 0.62, and 0.8. The experimental profile appears to demonstrate some hollowing of the profile near $ρ = 0.35 − 0.45$ that is not captured by the simulation. However, the infrequent sawteeth present in this discharge may contribute to this disagreement as their effects are not captured by gyrokinetic simulation.

Due to the presence of a carbon wall, DIII-D has a significant fraction of carbon present during all discharges and measurements of carbon density profiles are routinely made using charge exchange spectroscopy. To compare gyrokinetic predictions with these CER measured density profiles, the same approach used in the paragraph above was employed. Two trace carbon species were introduced into the flux-matched conditions to determine impurity diffusion and convection and the impurity peaking factor. The equivalence of the impurity peaking factor (−V/D) to the $1 / L n c$ when measurements are in steady state implies that the simulated V/D profile can be integrated to obtain a predicted steady state carbon profile and then compared directly with the measured carbon profile. This comparison is shown in Fig. 11.

In Fig. 11, the raw data obtained during the time averaging window are compared with the gyrokinetic predicted profile. Again, good agreement is found between simulation and experiment within the scatter of the data. However, it is worth noting a slight hollowing of the carbon profile may be present in the experimental data that occurs around 0.35–0.4; similar to observations found in the Li profiles. This feature is again missed by the simulated profiles and likely originates from the sawtooth activity that is cited in the explanation above. Inside of the last gyrokinetic simulation (⁠ $ρ = 0.3$⁠), the profile is extrapolated to a gradient scale length on axis equal to 0.0, so agreement/disagreement in that region should not be attributed to the model itself. However, the general conclusion of this comparison is there is favorable agreement between gyrokinetic simulation and experiment for carbon measurements (Z = 6).

To better understand the transport of both low and medium/high-Z impurities, the laser blow-off technique was applied to introduce trace levels (non-perturbing to background plasma parameters) of calcium (Z = 20) impurities into the conditions studied. The calcium impurities introduced via this technique were measured using multiple soft x-ray arrays⁵⁶ and charge exchange spectroscopy. Calcium is not generally fully stripped at the temperatures studied [ $T e ( ρ = 0.5 ) ∼$ 3 keV] and ionization and recombination rates from ADAS⁵⁷ were used in the impurity transport code Aurora to infer the experimental impurity diffusion and convection. Since the LBO injections are approximately a delta function source in time, with no other presence of calcium impurities typically present in DIII-D, the time evolution of the calcium impurities was modeled and compared with measurements to determine both the impurity convection and diffusion independently. The method and tools used for this analysis is best described in Ref. 53. From the simulation side, the trace species technique described above was used. However, we chose to specify the Ca impurity with a Z = 18 and A = 40 since it is generally the most abundant charge state of Ca over the radial locations studied. The resulting impurity diffusion and convection from experimental inference are compared with the simulation results in Fig. 12.

In Fig. 12, we find some clear disagreements between the simulated impurity transport coefficients and those found from experimental inference. Although the general shape of the diffusion coefficient profile is arguably well reproduced, there is up to a factor of approximately six difference between the magnitude of the impurity diffusion between modeling and experiment at some radial locations. In the convection, good agreement is found at inner radii (⁠ $ρ = 0.3$ and 0.43), but the results completely diverge at larger radii, with a large outward convection inferred from experiment for $ρ > 0.6$ and a large inward pinch predicted by simulation at these radii. These results demonstrate a very clear disagreement between the simulation and modeling in a region of that plasma where excellent agreement was found between fluxes, profiles, turbulence fluctuation measurements, and low-Z impurity transport.

Given the generally good agreement found in other quantities, the poor agreement with medium-Z impurity transport is surprising. However, the transport determined for the Ca impurities is inferred from experimental measurements as direct measurement of a single dominant charge state across the minor radius is not possible (as in the case of Li and C). Inference of impurity transport requires diagnostics with sufficient temporal and spatial resolution, accurate atomic physics data for the ionization, recombination, and photon emission rates, and generally relies on an assumption that D and V are charge state independent. The observed and simulated Ca transport (⁠ $D ∼ 10 m 2 / s$⁠) is quite high and is fairly poorly resolved by the $≥$ 5 ms time resolution of the CER measurement in this discharge. The transport of the impurities from the edge to the core following a laser blow-off injection occurs in only a few CER measurement times, likely leading to larger uncertainties in the inferred transport. Inferences are improved by the inclusion of SXR data, which has higher time resolution. However, SXR data are not charge state resolved and therefore is more reliant on accurate modeling of atomic physics. Unlike low-Z impurities that tend to exhibit a dominant, fully stripped, charge state over a wide radial range in the plasma, medium and high-Z impurities are partially stripped across most/all of the minor radius without a clearly dominant charge state in many regions. As a result, systematic (unquantified) uncertainties in atomic physics data play an important role in the accuracy of transport inferences. These known issues could potentially cause significant uncertainties with the inference of medium and high-Z impurities and are difficult to quantify. In addition to these uncertainties, the dependence of transport on impurity charge may also point to a potential cause of discrepancy and this is examined in more detail in Sec. V B.

To better understand multi-Z impurity transport in these conditions, and to potentially shed some light onto the origin of the discrepancy between calcium simulation and experiment, we present the results from a stand-alone simulation investigation performed on this target discharge. All results presented in this section are representative of turbulent contributions to impurity transport. Neoclassical contributions were evaluated with the NEO code for all impurities but were found to play a negligible role (typically <20× smaller) than the turbulent contributions in the region of interest (⁠ $ρ = 0.3 − 0.8$⁠). Unlike the analysis in Secs. II–IV, this analysis utilized the experimental density and temperature profiles obtained from experiment directly (as opposed to the multi-channel, flux-matched profiles described above), and we performed modifications to the experimental $a / L T i , a / L T e$⁠, and $a / L n$ manually to approximately match the experimentally derived ion heat fluxes and electron particle fluxes. This work was performed in this manner because the analysis predates the flux-matched profile simulations presented above. However, it is presented here because the results from this investigation shed light on the discrepancy found above and points to some interesting observations related to the Z dependence of impurity transport. In each simulation, we introduced trace impurities into each simulation that span a wide range of Z and A including: He (Z = 2, A = 4), Li (Z = 3, A = 6), F (Z = 9, A = 18), Al (Z = 13 A = 26), Ca [He-like (2 electrons remaining on the atom) with Z = 18, A = 40], and W (partially stripped, Z = 45, A = 184). These impurities were chosen because they were all introduced into the experimental condition via gas puffing, LBO, or the powder dropper and they span a wide range of impurity charge. The impurity transport coefficients (D and V) and impurity peaking (V/D) were extracted from simulations at $ρ =$ 0.3, 0.4, 0.5, 0.6, 0.7, and 0.8 and are plotted in Fig. 13. It should be noted that for all impurities except W, the charge to mass ratio was essentially $Z / A ∼ 0.5$⁠. Therefore, these results provide insight into predominately the dependence of impurity transport on charge and not charge to mass ratio.

The results from this analysis reveal some important information about the dependence of D and V on impurity charge, Z. Namely, it is shown in Fig. 13 that a dependence of impurity diffusion on charge exists in these conditions while there is no clear dependence of the impurity convection on Z. The dependence of impurity diffusion on Z and the general lack of dependence for impurity convection leads to a trend in the impurity peaking, V/D, where we find that at outer radii (⁠ $ρ ≥ 0.5$⁠) impurity profiles flatten as the impurity Z is increased. At inner radii, high Z impurities break this trend and appear to become more peaked. The trends in the impurity diffusion are found to vary with the radial location studied. This can best be seen by plotting the impurity diffusion as a function of Z at three different radii studied. In Fig. 14, the simulated impurity diffusion at $ρ =$ 0.3, 0.5, and 0.7 is plotted vs the simulated impurity Z.

As shown in Fig. 14, the impurity diffusion displays an inverse dependence on the impurity charge at inner radii (0.3 and 0.4) and then reverses at outer radii (outside of 0.4) to a positive dependence on the impurity charge. In all situations, the magnitude of the impurity diffusion asymptotes at large values of Z. The implication of this result is that low-Z impurity transport appears to vary significantly with impurity charge. In contrast, at higher Z, impurity diffusion between medium and high-Z impurities should be nearly indistinguishable. Although in the experiments presented here, we were unable to make a quantitative comparison with the simulated impurity diffusion trends; this clear trend in impurity diffusion may represent an attractive target for future experiments.

Due to operating with temperatures in the few keV range, medium and high-Z impurities in present tokamaks are only partially stripped. As a result, an impurity transport code is generally utilized to determine the transport coefficients and total impurity density for high-Z impurities. These codes (such as STRAHL⁵⁸ and Aurora⁵²) take in density and temperature profile information, utilize calculated ionization and recombination rate data, and calculate the charge state distribution in the plasma with the assumption that the impurity flux is $Γ = − D ∇ n + V n$ and that D and V are charge state independent. However, this assumption of charge state independence of the impurity transport is clearly invalid according to the modeling results plotted in Figs. 13 and 14 as the impurity diffusion is a strong function of Z, particularly at lower charge states that would exist in cooler, outer regions of the plasma. If the impurity diffusion is, in fact, a strong function of impurity charge, this would presumably lead to incorrect inferences of the impurity D and V from experimental data using typical analysis techniques and could contribute to the disagreement between simulated and experimentally inferred transport in the case of calcium. In the case of Li and C, one would expect errors to be smaller or negligible, as these impurities are fully stripped inside of the plasma edge and only a single charge state therefore spans most of the plasma minor radius. In the case of carbon, the density profile was measured directly and only steady state comparisons were made; therefore, there was no need for impurity transport inference to determine profile. It is therefore plausible that the disagreement in calcium impurity transport is, in part, due to the flawed assumption on impurity diffusion. We note that Aurora has been generalized to enable simulations with arbitrary Z dependence of transport coefficients. Reanalysis of these discharges would be required to demonstrate the impact of this effect and these investigations are left for future work.

The origin of the Z dependence of impurity diffusion has also been examined by taking advantage of modeled (non-experimental) plasma conditions. To investigate the dependence of impurity diffusion on Z for different turbulent types (ITG, $∇ n$ driven TEM, and $∇ T e$ driven TEM), we setup three different conditions based on the CYCLONE base case parameters.⁵⁹ This condition is well-known to exhibit strong ITG turbulence. The gradients use for this ITG dominated condition were $a / L T i = 2.5$⁠, $a / L T e$ = 0.1, and $a / L n = 0.8$ with an inverse aspect ratio of 0.36. However, to simulate different types of turbulence, modifications were made to this condition to destabilized primarily $∇ n$ driven TEM by setting $a / L T i = a / L T e$ = 0.1, $a / L n = 2.5$⁠, and $∇ T e$ driven TEM by setting $a / L T i = 0.1 , a / L T e$ = 2.5, and $a / L n = 0.8$⁠. A total of four trace impurities species were introduced into these simulations (Li, C, Al, W) to evaluate the impurity diffusion, convection, and peaking as a function of Z. The results of this exercise can be found in Fig. 15.

The results from this exercise shed light on the physical origins of the scaling of impurity diffusion with Z found in Fig. 14. As shown in Fig. 15, the impurity diffusion coefficient displays fundamentally different behaviors for ITG vs $∇ n$ driven TEM turbulence. ITG turbulence displays an inverse relationship of D with Z, while $∇ n$ driven TEM displays a positive dependence of D on Z. Both of these behaviors were observed in Fig. 14. The interpretation of these results in the context of the DIII-D experimental analysis is as follows: At inner radii, the plasma conditions are dominated by ITG turbulence leading to the inverse dependence of D on Z, while the balance of ITG to TEM shifts in the outer half of this condition and becomes dominated by $∇ n$ driven TEM, leading to the observed positive dependence of D on Z. Generally, a weak or small dependence of D on $∇ T e$ driven TEM is found. We note that this overall picture is consistent with the linear stability analysis shown in Fig. 2 and from results of analytical modeling of impurity transport in Ref. 60. The results from this reference link the Z dependence of the impurity diffusion to both the propagation direction of the turbulence and the charge dependence of the impurity drift frequency. An in-depth investigation into the physical mechanisms of the charge dependence is part of future work and are out of the scope of this paper.

Sections II–V presented an extensive validation of the nonlinear gyrokinetic model in a DIII-D ITER Similar Shape (ISS) discharge. A new surrogate modeling technique was coupled with nonlinear gyrokinetic simulations to perform some of the highest fidelity (including a wide range of physics effects and species) predictions of density and temperature profiles to date with a much reduced computation cost compared with traditional Newton solver based methods. Using this technique, Q_i, Q_e, and Γ_e were simultaneously matched by nonlinear gyrokinetic simulation over much of the plasma minor radius which allowed for direct comparisons of predicted T_i, T_e, and n_e profiles with experimental measurements inside of $ρ = 0.8$⁠. Excellent agreement between the simulated profiles and experiment was found within the scatter of the experimental data, suggesting that the ion-scale gyrokinetic model used to predict the profiles is an accurate model for the turbulence and transport found in the plasma conditions.

Comparison with experiment was extended beyond just heat and particle fluxes and made with more fundamental properties of the plasma turbulence. Measurements of low-k density fluctuations obtained from Beam Emission Spectroscopy and low-k temperature fluctuations obtained from Correlation Electron Cyclotron Emission (CECE) were compared directly with simulation at $ρ = 0.62$ and 0.8 utilizing a synthetic diagnostic developed for the CGYRO code. Agreement between simulation and experiment in density fluctuations was found to be excellent. Fluctuation levels between simulation and experiment were well within uncertainties and fluctuation cross-power spectra also demonstrated good agreement, with potential disagreements in the spectral peak at 0.62 resulting potential from uncertainties in the measured rotation at that location. The cross phase of the turbulence was also found to be in good agreement with measured BES. CECE measurements were also compared both in the measured fluctuation levels and the cross-power spectra. Fluctuation levels evaluated using synthetic modeling were found be under-predicted (by $∼ 2 x$⁠) relative to the experiment when using analysis techniques described here.⁵⁰ However, fluctuation spectra were in good agreement at 0.8 with less definitive comparisons possible at 0.62.

The final sections of this manuscript focused on multi-Z impurity transport. Impurity transport was obtained from the experiment and took the form of profiles of impurity peaking factors (−V/D) for Li impurities, measured density profiles for intrinsic C impurities, and experimentally inferred transport coefficient profiles (D and V) for Ca impurities. Values obtained from experiment were compared with transport determined by trace impurity transport obtained from the flux-matched CGYRO simulations. Low-Z impurities (Li and C) demonstrated excellent agreement with experiment with only some slight disagreement in the profile shape near the q = 1 surface. In contrast, significant disagreement was found comparing simulated diffusion and convective coefficients with experimentally inferred values for Calcium (He-like, Z = 18, A = 40). Although there are known sources of uncertainties in the inference of medium and high-Z impurity transport, part of this disagreement is potentially explained by examining the Z dependence of impurity diffusion and convection in these discharge conditions. By scanning the impurity Z from He to W using trace impurities in nonlinear CGYRO simulations, it was shown that this discharge exhibits a strongly Z dependent diffusion coefficient. The dependence of the diffusion coefficient on Z appears to switch from inversely proportional to Z at small radii to proportional to Z at larger radii. The physical origins of this were probed using variations around the CYCLONE base case parameters. It was shown that when ITG turbulence is dominant a inverse dependent of D on Z is found which reverses to a positive dependence when $∇ n$ driven TEM is dominant. The experimental results were therefore explained by a transition between dominant modes as a function of radius. More generally, the strong dependence of D on Z suggests that the assumption made in the impurity inferences, that D is independent of charge state, is flawed for these conditions and therefore could contribute to observed differences between simulated and experimentally inferred Ca transport.

This collection of work represents one of the most complete comparisons of experiment with gyrokinetic simulation performed to date. Comparison spanned fluxes (Q_i, Q_e, Γ_e), profiles (T_i, T_e, n_e), low-k density and temperature fluctuations (Beam Emission Spectroscopy and Correlation Electron Cyclotron Emission), and multi-Z impurity transport (Li, C, and Ca). Although some disagreements were found, in general, the agreement between simulation and experiment was found to be quite good suggesting that nonlinear gyrokinetics is indeed an excellent model for turbulence and transport in these conditions. This work presented the first validation exercise that utilized surrogate accelerated prediction of density and temperature profiles via nonlinear gyrokinetics. The results from this work are directly relevant for ITER itself and have important implications for prediction of ITER and future fusion reactors.

The authors would like to thank the DIII-D team for their support running and executing the experiments simulated in the paper. Specifically, we would like to thank A. Bortolon for operation of the Li powder dropper for these experiments and F. Scotti and T. Abrams for help planning the JRT experiments that were used to collect these data. This material was based upon work supported by the Department of Energy under Award Nos. DE-FC02-04ER54698, DE-SC0014264, DE-AC05-00OR22725, DE-SC0019256, DE-AC52-07NA27344, DE-AC02-05CH11231, DE-SC0017992, DE-SC0018287, and DE-SC0019352. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award <FES>-ERCAP<m224>. Some of the simulations performed in this paper were performed on the MIT PSFC Engaging cluster.

The authors have no conflicts to disclose.

Nathan Thomas Howard: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead). Pablo Rodriguez-Fernandez: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Christopher Holland: Formal analysis (supporting); Writing – review & editing (supporting). Tomas Odstrcil: Formal analysis (supporting); Writing – review & editing (supporting). Brian A. Grierson: Formal analysis (supporting); Writing – review & editing (supporting). Francesco Sciortino: Formal analysis (supporting); Writing – review & editing (supporting). George R. McKee: Formal analysis (supporting); Writing – review & editing (supporting). Zheng Yan: Formal analysis (supporting). Guiding Wang: Formal analysis (supporting). Terry Lee Rhodes: Formal analysis (supporting); Writing – review & editing (supporting). Anne E. White: Methodology (supporting); Writing – review & editing (supporting). Jeff Candy: Software (supporting). Colin Chrystal: Formal analysis (supporting).

The data that support the findings of this study will be made available through the PSFC library upon publication of this manuscript.