Energy transport analysis of NSTX plasmas with the TGLF turbulent and NEO neoclassical transport models

This work presents a study of plasma transport at low aspect ratio on the National Spherical Torus Experiment tokamak, where the turbulent and neoclassical energy fluxes calculated by the quasilinear Trapped Gyro Landau Fluid (TGLF) model and the multi species drift-kinetic Neoclassical solver (NEO) are validated against experimental data. The turbulent energy transport of two plasma discharges, one in the L-mode confinement regime and another in the H-mode regime, is dominated by electrostatic drift-wave instabilities, while the ion heat transport has a significant neoclassical contribution. The data analysis workflow is described in detail to understand how the variations of mapping and fitting of experimental data affect the power balance solution and subsequent flux-matching plasma profile predictions with the TGYRO solver. On average, the predicted plasma profiles are consistent with experimental data. However, the solutions are sensitive to various input parameters, including boundary conditions, and the electron-ion coupling. Linear gyrokinetic stability analysis demonstrates close agreement of the real frequencies of unstable modes between TGLF and CGYRO gyrokinetic simulations, but higher growth rates are predicted by TGLF, especially for the H-mode case. Estimates of the low-k modes’ contributions to the total flux are consistent with linear stability analysis and the E × B suppression of turbulence in TGLF simulations with the SAT1 saturation model, while the SAT2 saturation model over-predicts the low-k modes’ contribution in the H-mode case. Moreover, the results with SAT1 model are consistent with power balance analysis, which indicates only neoclassical ion energy fluxes inside ρ < 0.4 in the L-mode case and ρ⩽0.7 in the H-mode case. The presence of multi-scale turbulence and ion-scale driven zonal flow mixing effects are also observed in TGLF scans of the electron turbulent heat flux over a range of temperature gradients and the electron-ion temperature ratio, which could explain the strong model sensitivity to variations of input parameters.


Introduction
Understanding plasma transport is necessary to optimize plasma behavior in magnetic confinement devices and for confidence in transport predictions for current and future fusion devices.Comprehensive numerical transport models suitable for simulation of plasma processes that occur over a wide range of time scales require significant computational resources, which are typically not feasible for real-scale fulldevice modeling.However, with computationally fast reduced numerical models, it is possible to perform scans of various parameters to investigate the sensitivity of a given system to different plasma conditions.Identifying the regimes of validity where such reduced models can successfully reproduce experimental data increases the confidence in the predictive modeling.The design of the next generation fusion machines requires validation of transport models over a range of aspect ratio to determine the optimal parameters.
The Trapped-Gyro-Landau-Fluid (TGLF) model [1] is a reduced turbulence model for the prediction of transport driven by drift-wave instabilities, such as Trapped Electron Modes (TEMs) as well as Ion and Electron Temperature Gradient modes (ITGs, ETGs).TGLF computes turbulent fluxes using the linear eigenmode solution and a model of the saturated intensity multiplied by quasi-linear weights to complete the flux calculations.The saturation intensity model, which is often referred to as a saturation rule, is constructed to fit a large set of nonlinear gyrokinetic turbulence simulations with kinetic electrons.The original model for the Saturated Fluctuation Intensity (SAT0) [1] was 1-D with the intensity of each poloidal wavenumber only depending on the linear growthrate for that wavenumber.The saturation model was later extended to be a 2D model with multi-scale coupling between poloidal wavenumbers and radial wavenumbers through zonal flow mixing (SAT1) [2].Recently a 3D saturation model (SAT2) [3], that includes the poloidal angle dependence of the intensity due to geometry effects, has been fit to CGYRO turbulence spectra.
The TGLF model has been extensively utilized for transport modeling of conventional tokamaks with aspect ratio (R/a ≳ 3), where R is the major radius and a the minor radius [4][5][6][7][8][9][10][11][12].Using the Miller geometry approximation, the TGLF model takes into account effects of elongation and triangularity; therefore, it is not limited to a circular geometry and can be applied for the transport prediction of lower aspect ratio tokamaks with shaped geometries.However, transport mechanisms in low aspect ratio (R/a < 2) Spherical Tokamaks (STs) are different [13] compared to conventional tokamaks.
One of the main distinguishing features of STs is that the ion thermal transport is close to the neoclassical level as ion-scale low-k (k θ ρ s < 1, where k θ is the binormal wavenumber and ρ s is deuterium sound speed over gyrofrequency) turbulence is stabilized by the strong E × B shear induced by Neutral Beam Injection (NBI) ( [14] and references therein).Similar types of plasmas are observed on conventional tokamaks only in particular scenarios, such as the DIII-D high plasma β p scenario [15].In addition, the naturally high-β operation at low aspect ratio can introduce electromagnetic effects, driving Micro Tearing Modes (MTMs) and Kinetic Ballooning Modes instabilities (KBMs).
The validation of reduced transport models at low aspect ratio plasmas is still in its early stages and only a few studies have been reported [16][17][18].The under-prediction of high-k (k θ ρ s > 1) electron turbulent transport, in the first demonstration of TGLF on STs [17], has been determined to be due to limitations of the implementation of the collisional model and the saturation rule.Modeling performed later on the low aspect ratio National Spherical Torus Experiment (NSTX) [18], with the updated TGLF model, demonstrated better agreement with experimental data when simulations were performed with the more comprehensive SAT1 model to capture low and high-k modes' interactions, compared to simulations with the simplified SAT0 model.
In the work presented herein, an extension of previous studies is conducted to verify if turbulent heat transport predicted by TGLF and neoclassical heat transport, predicted by the multi species drift-kinetic neoclassical solver NEO [19,20], are consistent with experimental data from the NSTX low aspect ratio tokamak.Two plasma discharges, one in the L-mode confinement regime and another in the H-mode regime, were selected specifically for this study because their electron turbulent transport is predominantly driven by electrostatic drift-wave instabilities, according to previous gyrokinetic investigations [21].It was found that the growth rates of electromagnetic instabilities are below the E × B shearing rate.In both discharges, the ion heat transport has a significant neoclassical contribution typical for NSTX plasmas [22].A sensitivity study to understand how the experimental data analysis might affect the accuracy of predicted results is performed by comparing simulations with small variations of input plasma profiles.Linear stability analysis and scans of the turbulent fluxes over a range of plasma temperature gradients and electron ion temperature ratio performed with TGLF are analyzed to determine if a qualitative picture of mechanisms driving the transport are captured correctly by the TGLF and NEO models in these NSTX experimental conditions being considered.
The L-mode discharge studied here is 141716, an NBI heated plasma, which has been extensively analyzed via turbulence diagnostic measurements of a high-k scattering system and gyrokinetic simulations [14], and an initial validation study with the TGLF model [18].Furthermore, at low-β and density, as in this plasma, the drift-wave turbulence is electrostatic in nature [21], which avoids the complexity of electromagnetic effects.In addition, the E × B shear is relatively low in this plasma, providing conditions close to conventional tokamaks where both low-k ion scale and high-k electron scale turbulence are present.Gyrokinetic simulations used for the calibration of the TGLF saturation rules typically correspond to such plasma conditions.
The second discharge 129017 selected for this transport study is a typical H-mode scenario on NSTX or STs in general [23] with high-β and relatively high density.The set of shots 129015-129020, with very similar experimental parameters, have been investigated previously [21,24].The electron transport is expected to be dominated by high-k instabilities, according to gyrokinetic simulations [21].The ion transport is neoclassical from the plasma core to at least ρ ∼ 0.7, where ρ is the square root of normalized toroidal flux.The ionscale turbulence is suppressed by E × B shear, which is higher than the growth rates of these unstable modes.
The rest of this paper is organized as follows.Section 2 provides some experimental details of these selected NSTX discharges and the data analysis workflow.The description of the numerical set-up and the comparison of the predicted plasma profiles with the experimental data are presented in section 3. Results of the linear stability analysis and scans of the turbulent fluxes obtained with the TGLF model are discussed in section 4. Finally, a summary of the main results and discussion of possible paths toward predictive modeling of low aspect ratio plasma with the TGLF model are provided in section 5.

Experimental details and data analysis workflow
NSTX is a medium-size ST with the following main parameters: major radius R of 0.89 m, minor radius a of 0.61 m, with plasma current (I p ) up to 1.5 MA, toroidal field (B t ) up to 0.55 T and auxiliary Neutral Beam Injected (NBI) heating power up to 6 MW [23].Only NBI-heated shots have been analyzed in this work.
One of the main criteria for selection of the experimental data for TGLF model validation is that the thermal transport is not dominated or strongly affected by Magnetohydrodynamic (MHD) instabilities.The time of interest for the transport modeling in this paper is selected to avoid the influence of low frequency (<100 kHz) large amplitude instabilities (such as sawtooth, kinks, fishbones, tearing modes) on the heat transport.In all selected shots, the minimum safety factor (q) is above unity and no sawtooth activity has been detected on the measured neutron rate signals or the Mirnov coils signals.
Based on the time traces of Mirnov coils, which correspond to instabilities with toroidal mode numbers n = 1-4, the time intervals with amplitudes of magnetic fluctuations below 5 Gauss are selected for the transport analysis.
High frequency (>100 kHz) MHD instabilities (such as compressional and global Alfvén eigenmodes (GAEs)or toroidicity-induced shear Alfven eigenmodes) possibly can impact the electron core transport on NSTX [23,[25][26][27], in particular in high-power discharges.Although not shown here, the preliminary analysis of a Fourier spectrum of the high frequency magnetic fluctuations shows activity with the signature of GAEs for the H-mode discharge and some high frequency modes in the L-mode discharge.However, the experimental neutron and stored energy measurements are reasonably matched by transport code predictions with classical fast ions.An estimate of the amplitude of such modes requires a complicated numerical analysis, which is beyond the scope of this paper.These modes are localized deep inside the core, so they might affect only the small inner portion of the domain analyzed with the TGLF and NEO transport models.
Time traces of some relevant quantities for the analyzed discharges are shown in figure 1; the main plasma parameters (averaged over the times of interest) are presented in table 1 and in the following text.Experimental data is analyzed during the plasma current flat-top phase over a time window of 40 ms, which are shown by the vertical shaded regions (green for the L-mode 380-420 ms and blue for the H-mode discharges 430-470 ms) in figure 1(a).We do not analyze the Lmode plasma after 450 ms, corresponding to a more steadystate stored energy, due to strong MHD instabilities in this time interval, which could dominate the heat transport and are not predicted by the TGLF model.The flat-top plasma current and the magnetic field on-axis are higher in the L-mode discharge, compared to the H-mode (table 1).In both discharges, the lineaveraged plasma density is evolving in time due to the features of the density control system on NSTX, with the L-mode plasma changing significantly faster ( 1 Comparing the two discharges, the averaged density is approximately two times lower in the L-mode than the H-mode shot (figure 1(b)).The normalized plasma beta (β N ) and the injected power (P NBI ) in the Hmode discharge are more than double that of the L-mode discharge (table 1, figure 1(c)).In the L-mode plasma, the stored energy rapidly increases with time as 1 , however, it is lower than the inverse energy confinement time (16.6 1 s ), indicating that the speed of plasma energy change is slower than the plasma transport timescale.The plasma stored energy in the H-mode discharge is nearly stationary at the time of interest (figure 1(c)).The safety factor at the 95% normalized poloidal flux surface q 95 is 10 and 9 for the L-mode and H-mode plasmas, respectively (table 1).A standard normalized H-mode factor H 98y,2 [28] is below 1 for the L-mode case and around 1 in the H-mode discharge (table 1).
Figure 1(e) shows the plasma boundary (normalized poloidal flux ρ = 1) for these two selected shots, which have different magnetic configurations.The L-mode discharge is an inner-wall-limited discharge, while the H-mode shot has  a lower single null diverted configuration.The radial plasma profiles are measured by different diagnostic systems.The electron plasma temperature and electron density are from the multi-point Thomson scattering system, which measures from the high field side to the low field side, along the plasma midplane [29].The carbon temperature, density, and toroidal rotation profiles are determined by the Charge Exchange Recombination Spectroscopy diagnostic (CHERS) [30].The impurity concentration is estimated by the effective plasma charge (Z eff ) based on the CHERS measurements of carbon density, assuming that carbon is the main impurity, as is the case for the shots analyzed here.Radiated power profiles are measured by the bolometer diagnostic.The equilibrium reconstructions produced in this work are obtained through an EFIT [31,32]-TRANSP [33,34] integrated workflow using the OMFIT framework [35].Such integration provides a full kinetic equilibrium reconstruction, including information of total plasma pressure and current density profiles.The self-consistent full kinetic equilibrium reconstruction is obtained through several iterations of a cyclic workflow that first includes producing an initial equilibrium reconstruction with the EFIT code based on the magnetic diagnostics data, motional Stark effect measurements [36] and kinetic profiles measured by the Thomson scattering system.Then experimental profiles are mapped onto the flux space and together with the equilibrium are provided as input parameters to the TRANSP code, which computes the fast ion pressure contribution to the total plasma pressure with the NUBEAM module [37] as well as the plasma current density profile, including the bootstrap and beam driven currents.Outputs of transport codes improve the accuracy of the equilibrium reconstruction by providing constraints on the internal pressure and current density profiles to apply on the next iteration of the equilibrium reconstruction.In this work, a selfconsistent solution is usually obtained after 3-4 iterations of this workflow.It should be mentioned that in order to estimate how sensitive the predictions of the temperature profiles are to the equilibrium reconstruction, cases based only on magnetic constraints have been included in the analysis as discussed in the section 3.
To increase the reliability of experimental data by statistical averaging while keeping the temporal treatment of both the equilibrium and experimental measurements, the data is processed with the OMFITprofiles tool [38].First, the experimental profiles are mapped to the equilibrium and then convolved with a window function to the time frame of interest common for all diagnostic measurements.In the current study, the time frame is determined by the available input files for the EFIT code.The window function depends on the time resolution of the diagnostic measurements, and it has been set to 20 ms for the L-mode case and 40 ms for the H-mode discharge in order to include at least one experimental measurement inside this window.To fit the experimental data, the scale-length fitting method is used.This fitting method was determined to best fit the experimental measurement data in these discharges to provide the monotonic profiles required for a transport analysis and avoid spikes in transport coefficients, as well as negative gradients of plasma profiles, which are not physical in the conditions of present experiments.This fitting method is based on a least squares optimization of representations of the data using integrals of linearly connected scale lengths (the free parameters) at the user defined knots.To reach the optimal results between smooth fitting and the representation of the experimental profiles shape we limited the number of knots to three for the temperature profiles and to five for the density profile, which are equally spaced in both cases.In addition, we applied monotonic constraints to allow only positive scale lengths.The resulting fit corresponds to the case with the smallest least squared error among all fits with variations of the free parameters.
Data analyzed with the OMFITprofiles tool is shown in figure 2, where all profiles corresponding to available EFIT reconstructions during the times of interest (shown in figure 1) are plotted to demonstrate the evolution of the plasmas during these intervals.Circles represent the experimental measurements with their associated uncertainty and solid lines are the scale-length fitting results.A noticeable increase in time, in addition to relatively big uncertainty of measurements, is observed for the ion temperature in the L-mode case (figure 2(a)).The latter is explained by the relatively low amount of injected neutral beam power and consequent low intensity of plasma emission registered by the CHERS diagnostic.The electron density and temperature profiles are measured with low measurement uncertainty; however, they slowly evolve in time as the discharges are not strictly stationary (figure 2).The bump in the outer radius of the electron density in the H-mode at 469 ms (figure 2(g)) is attributed to the experimental measurements occurring during a different phase of an Edge Localized Mode (ELM) cycle.Toroidal rotation frequency profiles have low uncertainty and are stationary in time (figures 2(b) and (f ))).To improve the accuracy of the input data for the transport analysis, the profiles are averaged over the interval of interest of 40 ms and are referred to as experimental data for the remainder of the paper.
The comparison of the representative experimental profiles between the two selected discharges is shown in figure 3, where the shaded areas represent the variation of data inside the time averaging window.These profiles are an output of the time-dependent power balance analysis with the TRANSP code with the NUBEAM module, where the density and energy conservation equations are solved to determine the power deposition and transport coefficients.The experimental profiles: electron temperature (T e ), ion temperature (T i ), electron density (n e ), and toroidal plasma rotation frequency (Ω) are provided as an input, but since the TRANSP code does not support the propagation of the uncertainty, only the nominal fitted profiles shown on figure 2 are used as inputs for TRANSP.Therefore, the shaded regions depicted in figures 3(a)-(c) and (f ), which represent the evolution of plasma parameters in the time-averaged window are smaller compared to the range of values shown in figure 2. In order to estimate how the variations of the fitting method can affect the prediction of the temperature profiles, cases with a fitting method that represent the non-monotonic behavior of plasma profiles have been included in the analysis and some of the representative cases are presented and discussed in section 3.
The L-mode case has higher core temperatures, and steeper core density and temperature profiles compared to the H-mode case (figures 3(a)-(c)).The electron temperature is higher than the ion temperature at ρ ⩾ 0.6 in the H-mode case, and in the region 0.4 ⩽ ρ ⩽ 0.6 in the L-mode case.The plasma density profile in the L-mode case is almost half that of the Hmode case across the entire radial domain.The fast-ion density (not shown) is negligibly small compared to the main ion density across the radial domain.An effective plasma charge Z eff (profiles are not shown here), has an average value of 1.2 and increases toward the core in the L-mode case, signifying that the impurity density concentration is much smaller compared to the main ion density.In the H-mode case, the highest value of Z eff is around 2.3 at ρ ∼ 0.8 and decreases toward the plasma core.At the maximum Z eff , the ratio of the main ion density to the impurity density is n D /n C ⩾ 15. Between ρ = 0.2-0.4,values of angular momentum velocity are approximately 20% higher for the H-mode case compared to the Lmode (figure 3(f )).
Based on the power balance analysis with the TRANSP code with the NUBEAM module, the total electron energy transport (which is the integral of the divergence of the energy flux which represents the convection + conduction heat transport) is more than double the total ion power in both discharges (figures 3(d) and (e)).The convection losses are negligibly small compared to the conduction losses in all cases.The heat exchange term (Q ei ∼ neni(Ti−Te) T 3/2 e [39]) due to ion-electron collisions is comparable to the ion heat flux.Large uncertainty of the time-averaged fluxes is caused by the non-stationary plasma conditions and variations of the electron/ion temperature ratio in time, which is extremely sensitive to the temperature profiles' fitting, especially in the region where the absolute values are low.

Description of predictive modeling set-up
Due to a strong sensitivity of turbulent fluxes to profiles gradients when the plasma is above the instabilities threshold or the critical gradient, the results predicted by turbulent transport models are highly affected by the variations of the experimental measurements and mapping of the plasma profiles [40].The propagation of the uncertainty of profiles gradients into the transport model might produce large scattering of resulting fluxes and artificially large sensitivity to radial location.To address this issue, a flux-matching gradient is identified via the TGYRO solver [41].The TGYRO solver reconstructs plasma profiles based on gradients that satisfy the transport equation, which can be written in the following form for each plasma species where V ′ is the derivative of the flux surface volume with respect to radial coordinate r, Q is an energy flux and S is an integrated source, which has a dimension of power/volume.The sources are obtained in this paper from power balance analysis performed with TRANSP and a more detailed explanation of the solved equations and explicit form of all source terms can be found in [33,34].In the current study, where only the energy transport is being investigated, the sources for the electron transport include contributions from Ohmic P OH and auxiliary NBI heating P NBI , radiation P rad and ionization losses P ioniz and the collisional electron-ion heat exchange [39] minus energy gain (which is a time derivative of the energy) 3  2 d(n e T i )/dt, where n is the density and T is the temperature of the corresponding plasma species.Taken together they are equal to the divergence of the energy flux ∇Γ, which is purely diffusive due to a small contribution of the convection term.
S e = ∇Γ e = P OH + P NBI −P rad − P ioniz − P ei − 3 2 d (n e T e ) /dt , (2) Since the plasma profiles are provided as inputs to the power balance code, the energy gain contribution is calculated from the experimental data.A similar equation describes the sources for the ion heat transport with the main differences that the sources do not include contributions from the Ohmic heating, instead of radiation there are losses due to charge-exchange processes, and the collisional electron-ion heat exchange term has a different sign.
The energy heat flux Q is a sum of the neoclassical and turbulent contributions, which are calculated by the NEO and TGLF transport models, respectively, based on the local gradients of the temperature profiles and necessary input information: equilibrium parameters and other plasma profiles that are not being predicted.The TGYRO solver applies an iterative scheme varying the local gradients at each radial location to minimize the difference between the predicted fluxes (Q model ) and the target flux (Q target ).This target flux is not equal to the initial input source from the power balance (Q PB ) due to the self-consistent recalculation of the electron-ion collisional exchange term based on the temperature ratio.Plasma profiles are obtained by the integral of the local gradients at the radial locations ) where T * is the temperature at the boundary point r * .
The benefit of this flux-matching solution analysis is that in a situation with stiff profiles, the uncertainty in the fluxes gives a smaller error to the solution than the uncertainty in the gradients [40].However, it should be mentioned that sources can not be directly measured but are instead obtained through a power balance analysis (here using the TRANSP code) and therefore are subject to model bias.
In this work, TGYRO simulations are performed using eight radii (where each radius is a flux tube) in the range 0.2 ⩽ ρ ⩽ 0.7, unless otherwise specified.As will be shown below, the sensitivity of the plasma gradients to the changes in mapping rapidly increase toward the edge from ρ ∼ 0.75; therefore, the boundary location is chosen at ρ = 0.7.Setting the boundary at this location also minimizes the effects of neutrals on the results, where the model inside the TRANSP calculation of the neutrals density is not as high fidelity.An additional radial point at the magnetic axis is used as an inner boundary point where the fluxes are equal to zero.For simplicity, only the temperature profiles are predicted in these simulations; the density and rotation profiles are taken from experiment.Two ion species, deuterium and carbon, are used in these simulations with the assumption that carbon is the main impurity based on experimental data.Simulations are performed with the TGLF saturation model SAT1, unless otherwise specified, which is physically more comprehensive than SAT0 and more stable in terms of convergence than SAT2.Initial simulations of the selected NSTX discharges indicated that in some cases, especially relevant to STs, where the low-k turbulence is suppressed, simulations with SAT2 over-predict the contribution of the low-k instabilities in the turbulent flux and these results are shown in the section 4.1.Future work will be devoted to a more systematic validation of the SAT2 model for the low aspect ratio plasma, as this model has a potential to give more accurate results at the plasma edge, especially for strongly shaped plasma [3].Electromagnetic effects in the TGLF model are not included here because such settings do not significantly affect the predicted profiles in this analysis but instead increase the convergence time for the TGYRO solutions.

Comparison of predicted profiles with experimental data
For the prediction of plasma profiles, TGYRO requires the following inputs: magnetic topology (equilibrium), temperatures at the boundary point, sources that should be matched, and the experimental profiles that are held constant, such as densities and rotation.The interpretive integrated modeling of NSTX plasmas has demonstrated that increasing the complexity of the full kinetic equilibrium reconstruction and adding more constraints into the equilibrium solution increases the number of possible solutions.Even though data validation and verification are performed on many steps, there is still not enough experimental data with sufficient time and space resolution, and accuracy of calibration measurements to achieve a unique solution.Since the turbulence transport modeling is extremely sensitive to the local plasma values (equilibrium and plasma profiles), more attention to the validation metrics of the equilibrium solution and the power balance analysis is required, which will be the focus of future work.
To understand how the variations of experimental measurements (and equilibrium reconstruction) can affect the numerical solution, the following sensitivity studies have been performed.Instead of carrying out separate scans of some of the most important variables that determine the transport solution-safety factor, magnetic shear, rotation profiles, gradients of plasma density and temperatures-the predictions of plasma profiles with the TGYRO solver have been performed for different TRANSP runs (each with their own unique ID).These different TRANSP runs are the result of a self-consistent integrated solution, where small variations of the equilibrium reconstruction (which are obtained by various combinations of constraints or weight of different diagnostic measurements provided to the equilibrium solver) or fitting of experimental profiles has modified all dependent variables, through the integrated workflow described in section 2. Such self-consistent modifications of the input variables provide a more realistic picture of the model sensitivity, as they include coupling of all variables affecting the results of the equilibrium reconstruction, data analysis, and power balance calculations, which are the input parameters for the predictive models.One of the TRANSP runs (A63) published in previous investigations [14,18] is included in the analysis as well, verifying if our results are consistent with this old run, even though we have pursued analysis via different toolchains.
A selected set of input profiles, which are processed with the GACODE suite of tools [42] from TRANSP runs, are shown in figure 4 for the L-mode (left) and H-mode (right) cases.Each profile is labeled with its unique TRANSP ID and the vertical lines indicate the radial range where the predictive simulations are performed.The largest deviation among the experimental profiles (T e , T i , n e , Ω) is for the toroidal rotation in the region closest to the magnetic axis.The large oscillations of the rotation profile are due to a neoclassical correction to the impurity rotation profile applied in the TRANSP code.The workflow has been organized in the such way that the GACODE suite of tools derives the rotation profile from the radial potential of plasma obtained in the TRANSP code as a solution of radial force balance equation.Large variations of the ion heat flux between different IDs in the L-mode discharge (figure 4(e)) are a result of its very low absolute value: smaller than the energy gain term and comparable to the electron-ion heat exchange, which are very sensitive to the profiles' shapes (especially T e -T i ) as has been mentioned above.
The equilibrium quantities, namely safety factor and magnetic shear, which are expected to be some of the most sensitive parameters for the heat transport prediction, are shown in figure 5.
The predicted profiles and their flux-matching gradients corresponding to the input profiles from a particular TRANSP runID are in figure 6, where the gray area encompasses all values of input profiles.To compare the flux-matching gradients with the experimental values, the inverse scale lengths a/L Te , a/L Ti defined as a/L x = −(a/X)dX/dr minor are used.In most of the cases, the flux matching gradients (figures 6(c) and (g)) of the electron temperature are above the experimental level at the outer boundary point ρ = 0.7, especially for the L-mode case.This potentially indicates that TGLF underestimates the turbulent flux, as it needs much higher gradients to match the target flux level.Results in figure 6 also demonstrate that the TGYRO solution is sensitive to small variations of the input parameters inside their uncertainty intervals, such that predicted profiles can be lower or higher compared to experimental values.The ion temperature prediction is in good agreement with experimental values in the region ρ = 0.4-0.7 for both the L-(figure 6(b)) and H-mode cases (figure 6(f )), except for the one TRANSP runID A63, which had been obtained through a different workflow previously, where the prediction is below the experimental level.Inside ρ ⩽ 0.4,  across the entire radial domain (Z79 in figure 6(a); Z20, Z50 in figure 6(e)).
To verify that the toroidal rotation profile would not determine the results of temperature predictions, simulations were performed with all profiles corresponding to Z36 but the rotation profile from Z50 (figure 4(e)); they demonstrated no significant difference in the computed temperature profiles.This observation is consistent with results shown in section 4.1 that at the core radial locations the plasma gradients are very small and the growth rates of the unstable turbulent modes are much lower than the E × B shearing rate.Therefore, even variations of approximately 20% of the toroidal rotation velocity at ρ ∼ 0.2 do not significantly affect the predictions.A similar validation study has been performed for the variation of the safety factor profile, which demonstrated no significant effect on the predicted temperature profiles.This might be explained by the dominance of the turbulent transport in the predicted channels, which is not so sensitive to the variations of the safety factor at least in the range of parameters used in this study.
A correlation between the ion and electron temperature predictions inside ρ ⩽ 0.4 can be seen in figures 6(a)-(b) and (e)-(f ): the lower the electron temperature, the lower the ion temperature prediction.The electron and ion channels in modeling with TGYRO are coupled through the dynamic heat exchange term, which adjusts the target fluxes based on the electron to ion temperature ratio.In addition, the non-linear dependence of the electron and ion temperature gradients is included in the saturation rule of the TGLF model.While it is not possible to completely exclude this non-linear coupling, the prediction of each temperature profile separately will eliminate the impact of the deviation of one of the predicted profiles from the experiment on the prediction of another profile through the heat exchange term.It should be noted that in these simulations TGYRO still recalculates the dynamic exchange term on each iteration based on the shape of the predicted profile on the previous step.Therefore, there is always a moving target that TGYRO is trying to match, and the interpretation of the results becomes more complicated if several profiles are involved.
The results of the electron temperature profile prediction while the ion temperature is held fixed are shown in figures 7(a), (c), (e), (g) and vice versa for the ion temperature figures 7(b), (d), (f ) and (h).The variations of the predicted ion temperature profiles in both the L-mode and H-mode cases become more uniform: solutions obtained with different inputs are similar to each other and closer to experimental values.Such results indicate that the over-prediction of the ion temperature inside ρ ⩽ 0.4 in the H-mode plasma (figure 6(f )) and under-prediction in the L-mode plasma (figure 6(b)) are likely caused by the heat channels coupling in such a way that the transport models cannot accurately predict the temperature profile in this region.Although the diversity of predicted electron temperature profiles (figures 7(a) and (e)) visually looks similar to the previous simulations on figures 6(a) and (e), there is also an effect of the electrons-ions channels coupling: some of L-mode case profiles (Z58, Z79) are reversed from over-prediction as shown in figure 6(a) to under-prediction as shown in figure 7(a).
To quantify the total deviation of predicted profiles from experimental data, quantities introduced in [43]-RMS error σ T and offset f T -are calculated as follows where ϵ j = T s,j − T x,j is the deviation between the j th radial simulation point T s,j and the corresponding experimental point T x,j , where T is the local ion or electron temperature and N is the total number of simulation points; the boundary points are excluded.
Results are presented in figure 8 where bars labeled as T e (T e &T i ) and T i (T e &T i ) correspond to the profiles shown in figure 6, and bars labeled as T e and T i to the profiles shown in figure 7. Horizontal lines indicate the 25% level.The Hmode case has the tendency to over-predict the ion temperature  stronger than the electron temperature (figure 8).RMS error varies from 6% to 50% depending on the simulation setup.More under-predicted cases are observed in the L -mode discharge than in the H-mode (figure 8).The prediction of each temperature separately demonstrates lower error compared to the prediction of two profiles together.
To illustrate how the dynamic heat exchange term or deviation of the TGYRO target from the power balance source affects the predicted solution, the total heat fluxes are shown in figure 11.The electron heat flux obtained with TGYRO is entirely turbulent and in most of the cases, it is close to the power balance calculations (figures 11(a), (c), (e) and (g)).
If TGYRO found the target flux to be significantly higher or lower than the power balance (figures 11(c) Z84, (g) Z20), the temperature is under-predicted (figure 7(a) Z84)) or overpredicted, correspondingly (figure 7(e) Z20)).These results indicate that the TGYRO solution is predominantly affected by the accuracy of the sources obtained from the power balance analysis and the ion-electron flux ratio found to be in the TGYRO solution.For the electron heat flux an energy gain term due to the variation of the stored energy in the L-mode case (figure 1(d)) is around 25% compared to the conduction losses; therefore it might be a reason for additional concern in these TGYRO simulations aiming to predict the time averaged plasma profiles.However as seen from figures 11(c) and (b), there are already variations of the target fluxes in the range of about 25% in the present set of simulations, which can demonstrate the breadth of spreading of possible predicted solutions (figure 7).
As an additional validation exercise, TGYRO simulations where both temperature profiles are predicted for each time slice inside the time-averaged window have been performed.Offset and RMS errors, shown in figure 9, for these simulations are calculated based on the experimental and predicted profiles for each time slice.As an example, the predicted temperatures for TRANSP ID Z79 for each time slice inside the averaging window is shown in figure 10.The absolute values of error for the same time slices varies between predictions based on the different TRANSP IDs, but a common trend is an increase of the errors toward the end of the time-averaged interval can be seen.Additional simulations have been performed with TRANSP in the predictive setup (PTRANSP [44]) with TGLF and NCLASS modeling of the turbulent and neoclassical transport, respectively.This is a time-dependent  simulation where the flux solution depends on the previous time step, which explains the nearly same error for each time slice (figure 9(b)).The main advantage of using the TGYRO solver instead of the time dependent PTRANSP modeling are a much shorter computation time if only the selected time slices are analyzed and flexibility to adjust settings for the iterative scheme to control the convergence for each time slice.Another difference is the opportunity to use the NEO model, which is a more comprehensive neoclassical model compared to the NCLASS model [45].
The ion heat flux in the H-mode case is neoclassical, with a small turbulent contribution at ρ > 0.4 (figures 11(f ) and (h)), except in the TRANSP runID Z20 (figure 11(c)), where the turbulent contribution at ρ = 0.7 is significant.In the Lmode case, the turbulent contribution is dominant at radial locations ρ > 0.5 and decreases toward the core, where the ion heat flux becomes neoclassical (figures 11(b) and (d)).These results are consistent with the power balance analysis shown in figure 12, where the ion heat diffusivity inferred from experimental profiles (χ i ) is compared with the neoclassical diffusivity calculated by NCLASS within the TRANSP code.In the H-mode case, the ion heat diffusivity is at approximately the neoclassical level (figure 12(b)).For the L-mode case, it is difficult to make a confident estimate of the level of experimental diffusivity compared to the neoclassical level due to the high variation in the evolution of plasma parameters in time, but it is reasonable to expect that the ion heat diffusivity might be significantly above the neoclassical level at ρ > 0.4 (figure 12(a)).Dimensionless collisionality for both discharges is around 0.2-0.8 and based on [46] we can conclude that it corresponds to a plateau regime.
To estimate the influence of plasma effective charge on the prediction of the temperature profiles a scan of Z eff has been performed by multiplying the carbon density profile by a constant value and recalculating the ion and impurity concentration in the input files to keep the quasi-neutrality of the plasma.Results presented on figure 13 demonstrate increase of the predicted temperatures with increase of Z eff and this dependency is stronger for the H-mode case.A weak dependency in the L-mode case can be explained by the low absolute value of Z eff in this discharge.Need to be noticed, that in these simulations two temperature profiles are predicted together and the

Linear stability analysis
To better understand the turbulent transport shown in the previous section, a linear stability analysis has been performed with the TGLF model at various radial locations with an assumption of a flux tube geometry in the Miller approximation.The linear stability analysis provides linear growth rates and real frequencies of unstable modes as a function of the wave-number.The growth rates (γ) are in normalized units c s /a unit , where c s = √ (T e /m i ) is the ion sound speed and a unit is the outboard midplane minor radius of the separatrix which is used as the unit of length.The sign of the real frequency (ω) identifies the direction of the mode propagation.In the TGLF convention, a negative frequency corresponds to propagation in the ion diamagnetic direction and a positive real frequency means propagation in the electron diamagnetic direction.Modes that are unstable in the wave-number range k θ ρ s ⩽ 1 are labeled in this article as low-k modes.Typically, they are associated with ion-scale turbulence.Unstable modes with wavenumbers k θ ρ s ⩾ 1 are labeled as high-k modes and typically associated with electron scale turbulence.Different types of modes, namely ITG (Ion Temperature Gradient), TEM (Trapped Electron Modes), MTM (Microtearing Mode), KBM (Kinetic Ballooning Mode), and ETG (Electron Temperature Gradient) can occur within these wavenumbers from low-k to high-k and have been observed in gyrokinetic or gyrofluid simulations.A justified identification of different types of modes is beyond the scope of this paper as it would require a lot of scans of gyrokinetic simulations due to complexity of the mode spectra on NSTX and existing of hybrid modes such ITG/KBM and TEM/MTM/KBM [13,47].Therefore, only the relative contribution of low-k versus high-k of the total turbulent flux is discussed The linear stability analysis also identifies the range of the wavenumbers where the modes are unstable and their propagation direction.
The results of the TGLF linear stability analysis performed on the TGYRO flux-matched profiles for one of the selected TRANSP runIDs in the L-and H-mode cases are shown on figure 14 for radii from ρ = 0.3−0.7.The amplitude of growth rates increases toward the edge (figures 14(e)-(h)) in both discharges indicating that instabilities become stronger as the radii increase.Larger amplitude of growth rates in the L-mode case can explain lower temperature in the L-mode discharge, compared to the H-mode in the radial range ρ > 0.4 shown in figure 3. Inside ρ ≲ 0.39, the amplitude of growth rates is relatively small, particularly for the range of low wavenumbers k θ ρ s < 1 (figure 14(e)).The plasma gradients are too small to drive sufficient turbulent transport at this radial location and plasma conditions are close to the stability limit of electrostatic drift-wave turbulence (shown in section 4.1.2);the TGLF model (as probably any turbulent model) is less accurate in such plasma conditions.
The growth rate amplitudes are similar between the L-mode and H-mode cases at radial locations up to ρ ⩽ 0.3 and modes propagate predominantly in the electron diamagnetic drift direction.Moving toward the edge of the plasma, the growth rates of modes at low-k propagating in the ion diamagnetic direction appear, first in the L-mode plasma at ρ > 0.3 then at higher radius at ρ > 0.4 in the H-mode plasma (figures 14(e)-(h)).In the L-mode discharge, the amplitude of these modes exceeds the E × B shearing rate (γ E×B = −(r/q)dω 0 /dr, where ω 0 = c s E r /RB pol is an angular frequency) at ρ ⩾ 0.4 (figure 14(f )), while for the H-mode discharge it is observed only at ρ ⩾ 0.6  (figure 14(g)).Therefore, up to these locations, unstable modes in the low-k portion of the spectra are expected to be weak or completely suppressed.The growth rates of the high-k modes propagating in the electron diamagnetic direction are above the E × B shearing rate over the entire radial domain, and possibly dominate the transport at least at the locations where the low-k modes are suppressed.Starting from ρ ⩾ 0.45, modes propagate deeper (in terms of wavenumber) in the ion diamagnetic direction in the H-mode case compared to the L-mode (figures 14(c) and (d)).

Comparison with CGYRO analysis.
For more detailed verification of the accuracy of the linear eigenmodes solver of the TGLF model for this work, the results are compared with linear gyrokinetic simulations performed with the CGYRO code [48].Many such validation studies have been done before for conventional tokamaks [3,49,50], while here the first comparison with gyrokinetic simulations for a low aspect ratio NSTX plasma is presented.The comparison is done for the H-mode discharge corresponding to TRANSP RUNID Z37 and L-mode corresponding to TRANSP RUNID Z79.
No CGYRO solutions have been found inside ρ ⩽ 0.4 for the H-mode case, which confirms that plasma gradients are close to the threshold of a stability limit at this region.The comparison of growth rates and real frequencies obtained on the experimental profiles and TGYRO solution profiles at two different radial locations for each discharge are shown in figures 15 and 16.
In both discharges, a good agreement is observed between real frequencies in electrostatic simulations.Including the electromagnetic effects (parallel vector potential (δA || ) and parallel magnetic field (δB || ), in addition to the electrostatic potential (δϕ) that is present in electrostatic simulations) does not change the computed frequencies or growth rates in the L-mode case, except for a small range of low wave-numbers k θ ρ s ≲ 0.4 in the TGLF simulations (figure 15).For the H-mode case, electromagnetic modes propagating in the electron drift direction are present only in CGYRO simulations at the range of wave-numbers 0.5 ≲ k θ ρ s ≲ 2 (figure 16).The TGLF electromagnetic simulations can properly resolve only (δA || ) perturbations, but discrepancy with CGYRO simulations with similar settings are observed in all cases of the Hmode simulations.
Regarding the suppression of low-k instabilities, the CGYRO simulations of different radial locations show results similar to those TGLF: for the H-mode case, growth rates of low-k modes are below or close to the value of E × B shear at radii up to ρ ∼ 0.6 (figure 16), while for the L-mode case they are above of E × B shear (figure 15).In both discharges TGLF model predicts higher growth rates than CGYRO and this effect is especially pronounced in the H-mode case.This result of over-prediction of growth rates in TGLF compared to CGYRO is similar to what was obtained in [50] for the DIII-D H-mode plasmas, where TGLF systematically predicted higher growth rates compared to simulations with the gyrokinetic code GYRO [51].A recent work comparing TGLF linear growth rates with CGYRO linear simulations on a large database of DIII-D cases also shows this similar trend: growth rates predicted by TGLF are higher than CGYRO simulations [52].
Increase of the growth rates in the range of low-k wavenumbers and up to k θ ρ s ∼ 2 in CGYRO electromagnetic simulation compared to electrostatic simulations (figures 16(g) and (h)) indicates the destabilization of electromagnetic turbulent modes.This is a similar behavior reported in previous works [13,47] demonstrating the range of hybrid modes like TEM/KBM and ITG/KBM observed in NSTX plasmas.In addition, as seen on the figures 16(g) and (h) growth rate is underestimated in the low-k range when only (δA || ) is included, which confirms the importance of full electromagnetic simulations for the STs plasmas revealed in [13,47].
Electromagnetic effects start playing important roles as the plasma β e increases.In CGYRO, the relevant quantity is the normalized electron beta β e,unit = 8π n e T e /B 2 unit where B unit = (q/r)dψ /dr [48].Values of β e,unit are shown in the table 2 for two different radial locations for both experimental and TGYRO solution.As expected, the H-mode has much higher β e,unit values and, therefore, electromagnetic effects are more important as described above.However, the threshold of the importance of electromagnetic effects is likely dependent on various other parameters that depend on ρ such as collisionality, safety factor, and gradients as demonstrated through various scans in previous linear stability analysis of different NSTX shots [13,21,53].

Analysis of turbulent fluxes.
To calculate the total turbulent flux in TGLF a saturation rule is applied, which multiplies the quasi-linear fluxes at each wave-number with a saturation intensity whose free parameters are adjusted to fit the database of non-linear gyrokinetic simulations.Investigation of the contribution of low-k and high-k modes into the total turbulent flux is performed to verify if these results are consistent  17.Fluxes are measured in GyroBohm normalized units (Q GB = n e T e c s (ρ s /a) ); therefore from ρ ∼ 0.6, the amplitude of fluxes rapidly increases with radius, especially in the L-mode case.
Inside ρ ⩽ 0.4, the turbulent electron heat flux is driven by high-k modes in both the L-mode and H-mode cases (figure 17 , and the ion turbulent heat flux predicted by TGLF is close to zero at all radial locations (figures 17(f )-(h)) in simulations with SAT1 model.In both cases, for the L-mode and H-mode discharge, the contribution of low-k modes is bigger in simulations with SAT2, which even predicts significant turbulent ion heat flux at ρ > 0.55 in the H-mode case.
Figure 18 shows the contribution of low-k modes into the total turbulent electron and ion heat fluxes calculated as a percentage of the flux integrated over the range of low wavenumbers k θ ρ s ⩽ 1 to the total flux integrated over the entire wave-spectrum.Values in the range 0.2 ⩽ ρ < 0.3 as well as the ion heat flux for the H-mode case are not shown, because the amplitude of the total predicted flux is too low to confidently utilize the results of the turbulent model in this range.As seen from this plot, there is a large contribution (up to 80%) of the low-k modes into the total electron heat flux in the L-mode case inside of 0.45 < ρ < 0.7.Large fluctuations of the percentage of the total contribution might indicate that instabilities driven by low-k modes are close to the stability threshold or the TGLF model is more sensitive to plasma parameters if multi-scale turbulence is present.The turbulent ion heat flux in the L-mode case is predominantly driven by low-k ion scale instabilities at all radial locations.
Overall, the picture of the low-k modes contribution to the total flux is consistent with a linear stability analysis showing the amplitude of growth rates in comparison to the E × B  shearing rate.Moreover, these results are consistent with a power balance analysis demonstrating the neoclassical level of the ion heat transport inside ρ < 0.4 in the L-mode case and ρ ⩽ 0.7 in the H-mode case.Results from [18], where local gyrokinetic simulations show a large contribution of the ion scale turbulence into the total electron heat flux at 0.55 ≲ ρ ≲ 0.7 in the L-mode case are also in line with observations presented here.No contribution of the low-k modes into the electron heat flux in the H-mode case at 0.6 < ρ < 0.69 might be the result of the strong dominance of high-k instabilities in these plasmas.Another possible reason might be an underestimation of low-k modes' contributions by TGLF at these plasma conditions, but future dedicated systematic comparisons with non-linear gyrokinetic simulations are required in this case.
A scan of the plasma temperature gradients and the electron-ion temperature ratio can demonstrate how close the plasma conditions are to the stability threshold and how steep the dependency of the flux is on these parameters.Such results are useful to understand how sensitive the model is to the uncertainty of experimental measurements.The scans of the electron turbulent heat flux are shown in figure 19 for the L-mode plasma (a)-(d) and for the H-mode plasma (e)-(h).The vertical line shows the experimental conditions and the x-axis represents the percentages of the variation of the scanned parameter from the experimental value.The lack of a dependency of the electron heat flux on the ion temperature gradients (a/L Ti ) or the ion-electron temperature ratio (T e /T i ) indicates that plasma conditions are below the instability threshold, as seen for the H-mode case at ρ = 0.35 (figure 19(e)).Instabilities driven by the electron temperature gradient (a/L Te ) are destabilized only at gradients above 25% of experimental level, which leads to a sharp increase of the turbulent flux with increasing a/L Te .At all other radial locations for both H-mode and L-mode, the dependency of the flux on plasma gradients or temperatures is more complicated: the electron heat flux can sharply increase/decrease with both a/L Te and T e /T i and not so strongly depend on a/L Ti (figures 19(a), (f ) and (g)).At these radial locations, the linear stability analysis shows the suppression of the low-k modes by E × B shear in the corresponding L-or H-mode scenario (figures 14(e) and (f )).A drop of Q e with the increase of a/L Ti (figures 19(b)-(d) and (h)) is a result of ion-scale driven zonal flow mixing effect which is in line with an evidence of unstable low-k modes from the linear stability analysis (figures 14(g) and (h))).
Overall, the scans demonstrate the complicated dependency of the electron turbulent flux on the plasma parameters, in particular if multi-scale turbulence is present.For a confident prediction of such plasmas, more detailed validations against the non-linear gyrokinetic simulations would be required to validate the correct representation of multi-scale turbulence iterations, the transition from low-k to high-k turbulence and the zonal flow mixing effects.
An additional scan has been performed to investigate the influence of Z eff on turbulent fluxes, which is presented in figure 20.As seen in the figure, no strong dependency of turbulent fluxes on Z eff in the range close to the power balance estimates has been observed.

Conclusion
The focus of this work is to understand how accurately the TGLF and NEO models reproduce experimental results on the NSTX low aspect ratio tokamak.Two specifically selected shots, where the electron transport is dominated by electrostatic drift wave instabilities and the ion heat transport has a significant neoclassical contribution have been selected for this study.The description of the data analysis workflow demonstrates that the non-stationary plasma conditions and a large sensitivity of the electron-ion heat exchange term to the fitting of the temperature and density profiles are the main reasons for the uncertainty in the heat flux calculations.
A comparison of the predicted temperature profiles with experimental data shows that statistically the results are consistent with experimental measurements; however, the numerical results are sensitive to input parameters.Therefore, the RMS error of the deviation of the predicted profiles from the experiment varies from 6% to 50%, depending on the simulation setup.The shortcomings of the TGYRO simulations are amplified when both temperature profiles are modeled instead of a single temperature profile.One of the parameters that affects the results of the prediction is the electron-ion channels' coupling and the modification of the target flux from the power balance estimates.In particular, it strongly affects the ion temperature prediction described by the NEO model.As a result, the over-prediction of the ion temperature inside ρ ⩽ 0.4 in the H-mode and under-prediction in the L-mode are likely caused by the heat channels coupling, in such a way that the transport models cannot accurately predict the temperature profile in this region.For the electron temperature prediction, a significant deviation from the experimental data is observed when there is a large discrepancy between the target and the power balance fluxes.
The linear stability analysis demonstrates good agreement in real frequencies of unstable modes between TGLF results and CGYRO gyrokinetic electrostatic simulations for both discharges; however, higher growth rates are predicted by TGLF.Electromagnetic effects do not affect computed frequencies or growth rates in TGLF simulations, except for a small range of low wave-numbers k θ ρ s ≲ 0.4.The electromagnetic effects in CGYRO simulations do not change the computed frequencies or growth rates in the L-mode case.In contrast, in the simulations of the H-mode case, CGYRO predicts higher growth rates, which become comparable with the TGLF results and electromagnetic modes propagating in the electron drift direction are present at the range of wave-numbers 0.5 ≲ k θ ρ s ≲ 2.
The estimate of the low-k modes' contribution to the total flux is consistent with a linear stability analysis and the suppression of turbulence by E × B shearing in TGLF simulations with SAT1 model.Moreover, the results obtained with SAT1 model are consistent with a power balance analysis demonstrating the neoclassical level of the ion heat transport inside ρ < 0.4 in the L-mode case and ρ ⩽ 0.7 in the H-mode case.Results obtained with the SAT2 saturation model over-predict the low-k modes' contribution in the H-mode case.A presence of multi-scale turbulence and ion-scale driven zonal flow mixing effects are observed on the scans of the electron turbulent heat flux over a range of the temperature gradients and the electron-ion temperature ratio.
Paths toward reliable predictive modeling of low aspect ratio plasmas with the TGLF model should include more accurate power balance analyses, which determines the target fluxes of the flux-matching solution.In addition, more accurate power balance analysis is required for the analysis and prediction of non-inductive scenarios on NSTX, where the bootstrap current plays a large role, and is dependent on the details of the predicted profiles.More systematic calibrations of the saturation model on non-linear gyrokinetic simulations for a large set of experimental data from low aspect ratio tokamaks would verify the applicability of the saturation rule models for calculations of the total turbulent flux for low aspect ratio plasmas, especially if multi-scale effects are present.Higher plasma gradients and plasma-β, which are expected to be present in plasma of future ST power plants, will increase the influence of electromagnetic instabilities on thermal transport; therefore, models with comprehensive physics of electromagnetic instabilities would be required.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Figure 1 .
Figure 1.Time traces of (a) plasma current Ip, (b) line-averaged density ⟨ne⟩, (c) NBI injected power P NBI (d) stored energy W for the selected NSTX discharges: L-mode (solid green), H-mode (blue dashed dotted).Shaded areas indicate the times of interest where transport analyses have been performed: green for the L-mode and blue for the H-mode.(e) A representative plasma boundary obtained through the full kinetic equilibrium reconstruction workflow.

Figure 2 .
Figure 2. Representative experimental profiles for both cases, L-mode on the left (a)-(d) and H-mode on the right (e)-(h), analyzed with the OMFITprofiles workflow demonstrating the variation of plasma quantities during the time of interest: (a), (e) ion temperature profiles T i , (b), (f ) toroidal plasma rotation frequency Ω, (c), (g) electron density ne, (d), (h) electron temperature Te.Circles represent the experimental measurements with their associated uncertainty and the solid lines are the scale-length fitting method fits.The presented time frame corresponds to the available EFIT reconstructions.

Figure 3 .
Figure 3. Kinetic profiles of electron temperature (a), ion temperature (b), electron density (c), total energy transport of electrons (d) and ions (e), and rotation profile (f ).The profiles in (d) and (e) are the results of power balance analysis with the TRANSP code and NUBEAM module and averaging over the interval of interest.Variations of data due to the evolution in time inside the time-averaged window is shown by the shaded area.

Figure 4 .
Figure 4. TGYRO input profiles processed with the GACODE suite of tools from TRANSP results: electron temperature (a), (g), ion temperature (b), (h), electron density (c), (i), total electron heat power (d), (j), total ion heat power (e), (k), and toroidal rotation profile (f ), (l), are shown for the L-mode (left) and H-mode (right) cases.Each profile corresponds to a unique TRANSP ID and the vertical lines indicate the radial range withing which the predictive modeling is performed.

Figure 5 .
Figure 5. TGYRO equilibrium input profiles processed with the GACODE suite of tools from TRANSP results: safety factor (a), (c), and magnetic shear (b), (d) are shown for the L-mode (left) and H-mode (right) cases.Each profile corresponds to a unique TRANSP ID and the vertical lines indicate the radial range within which the predictive modeling is performed.

Figure 6 .
Figure 6.Predicted plasma temperature profiles (a), (b), (e), (f ) and flux matching gradients (c), (d), (g), (h) obtained with TGYRO (TGLF + NEO) solver for the L-mode (left) and H-mode (right) discharges, when the ion and electron temperature profiles are predicted together.The gray area encompasses values of all the input profiles corresponding to different TRANSP runIDs from figure 4.

Figure 7 .
Figure 7. Predicted plasma temperature profiles (a), (b) and flux matching gradients (c), (d) obtained with TGYRO (TGLF + NEO) solver for the L-mode (left) and H-mode (right) discharges, when ion and electron temperature profiles are predicted separately (holding the other temperature profile fixed).The gray area encompasses values of all the input profiles corresponding to different TRANSP runIDs from figure 4.

Figure 8 .
Figure 8. Deviation of the predicted temperature profiles from the experimental values measured as RMS error and offset.Bars labeled as Te(Te&T i ) and T i (Te&T i ) correspond to the profiles shown in figure 6, and bars labeled as Te and T i to the profiles shown in figure 7.

Figure 9 .
Figure 9. Deviation of the predicted temperature profiles from experimental values measured as RMS error and offset obtained in simulations of separate time slices for the L-mode discharge for selected TRANSP IDs.Temperature profiles have been predicted with TGYRO solver (a) and in TRANSP predictive simulations (PTRANSP) (b).

Figure 10 .
Figure 10.Predicted plasma temperature profiles obtained with TGYRO solver for TRANSP runID Z79 (a)-(d) and TRANSP predictive simulations equivalent to Z79 (e)-(h) for each time slice inside the averaged window.

Figure 11 .
Figure 11.TGYRO total energy flux (solid) and the turbulent contribution (dashed) into total energy flux are shown for the simulations when two temperatures have been predicted-first row (a), (b), (e), (f ) and for simulations where each temperature profile has been predicted separately-second row (c), (g), (h).The different colors correspond to different TRANSP run inputs as given in legend.The gray bands correspond to the range of power balance fluxes across all TRANSP runs.

Figure 12 .
Figure 12.Ion heat diffusivity inferred from experimental profiles (χ i ) in comparison with neoclassical values calculated by NCLASS during the power balance analysis with TRANSP code.Results correspond to the TRANSP ID Z79 for the L-mode (a) and Z37 for the H-mode (b).The bands correspond to the variation in given quantity.

Figure 13 .
Figure 13.Temperature profiles predicted with TGYRO solver by scaling the carbon density and keeping the quasi-neutrality of plasma.
(a)), and the predicted TGLF turbulent ion heat flux is negligibly small (figure 17(e)).The low-k modes start contributing to the total heat flux at ρ ∼ 0.4 and their contribution increases with radius moving towards the edge for both the ion and electron fluxes in the L mode case (figures 17(b)-(d)), while the ion heat flux is driven entirely by the low-k modes (figures 17(f )-(h)).In contrast, in the H-mode case the electron flux is driven predominantly by the high-k modes up to ρ ∼ 0.69 (figures 17(b)-(d))

Figure 17 .
Figure 17.Turbulent electron Qe (left) and ion Q i (right) heat fluxes obtained with the TGLF model (SAT1 and SAT2) on the selected TGYRO predicted profiles for L-mode (green) and H-mode cases (blue) at different radial location.

Figure 18 .
Figure 18.Contribution of low-k modes into the total electron and ion turbulent heat flux for the L-mode and H-mode cases.(a) Simulations obtained with SAT1 (b) Simulations obtained with SAT2.

Figure 19 .
Figure 19.The turbulent electron heat flux as a function of electron temperature gradient (a/L Te ), ion temperature gradient (a/L Ti ) and the electron-ion temperature ratio (Te/T i ).The scans for the L-mode case on different radial locations are shown on the first line (a)-(d) and scans for the H-mode case are shown on the second line (e)-(h).The vertical line shows the experimental conditions and the x-axis represents the percentages of the variation of the scanned parameter from the experimental value.Stars are the power balance estimation of fluxes.

Figure 20 .
Figure 20.The turbulent electron and ion (only (d)) heat fluxes as a function of the scan of Z eff value.The scans for the L-mode case on different radial locations are shown on the first line (a)-(d) and scans for the H-mode case on the same radial locations are shown on the second line (e)-(h).The vertical line shows the experimental conditions and the x-axis represents the percentages of the variation of the scanned parameter from the experimental value.Stars are the power balance estimation of fluxes.

Table 1 .
Average plasma parameters at the time of interest in the selected NSTX shots.

Table 2 .
β e,unit values for L and H cases on two different radial locations.B shear and zonal flow mixing caused by interaction of the ion and electron scale modes is captured correctly by the TGLF model in the considered NSTX plasma conditions.Spectra of the total turbulent electron and ion heat fluxes obtained with the TGLF model on the TGYRO flux-matched profiles are shown for different radial locations in figure