Successful prediction of tokamak transport in the L-mode regime

A long standing shortfall in the predicted L-mode edge energy transport by reduced quasi-linear models of gyrokinetic turbulent transport has been resolved. The improved model TGLF-SAT2 has higher fidelity to gyrokinetic simulations of the electron-scale contribution to the electron energy transport and the ion-scale flux surface shape dependence of energy transport. The success of TGLF-SAT2 in predicting the L-mode and Ohmic edge profiles is critical to whole pulse simulation and opens the door to prediction of the H-mode power threshold.


Introduction
The International Tokamak Physics Activity (ITPA) Transport and Confinement Topical Group established a joint research activity (TC-10) to investigate the cause of the underprediction of the L-mode transport near the last closed flux surface of tokamaks in 2011.This presentation marks the closing of this activity by documenting the success of non-linear gyrokinetic turbulence simulations and a new saturation model for quasi-linear transport models (SAT2) in predicting the L-mode transport.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
In this paper, a survey of the history of the L-mode shortfall problem will be given.A companion paper by Angioni et al [1] presents detailed validation, analysis and predictions using TGLF-SAT2 on a large range of L-modes from ASDEX Upgrade.In section 2 the general observation of a shortfall in the predicted ion and electron energy transport will be illustrated [2] with the quasi-linear trapped gyro-Landau fluid (TGLF) transport model [3] using the SAT0 saturation model [4].In section 3 it will be shown that the shortfall in ion-scale transport observed in global gyrokinetic turbulence simulations [5,6] appeared to confirm the quasi-linear TGLF-SAT0 results and point towards a missing instability, since the measured density fluctuation intensity was also higher than the gyrokinetic simulations [7][8][9].A shortfall in the TGLF-SAT0 energy transport, just inside the separatrix, in Ohmic plasmas was also found [10,11] to make the prediction of the internal inductance too inaccurate for modeling the plasma control system.However, gyrokinetic turbulence simulations in a periodic flux tube, that are fully spectral in the radial wavenumber, largely resolved the ion energy transport shortfall, giving agreement with experiment within the measurement error [12][13][14].A path to resolving the electron energy transport shortfall was found when gyrokinetic simulations of L-modes, with real electron and ion masses [15,16], showed that the electron scale transport is suppressed by zonal fluctuations driven by the ion scale turbulence, but that this suppression can be overcome, leading to enhanced electron energy transport.Importantly, the quasi-linear approximation [17], that the phase shift of the turbulent energy flux is well captured by the phase shift of the most unstable linear modes, was shown to be satisfied by the L-mode gyrokinetic simulations [12,13], so a reduced quasi-linear transport model could be constructed, if a better fidelity model of the saturated fluctuation intensity could be found.In section 4 the improvements to the saturation model for TGLF that captured the multi-scale interactions (SAT1 [18]) and improved the fit to fully spectral gyrokinetic simulations (SAT2 [19,20]) are briefly summarized.In section 5 a set of stationary L-modes from DIII-D are presented showing the improvement in the predictions with TGLF-SAT2.The final section 6 examines how the resolution of the L-mode edge shortfall presented in this paper impacts the future of whole plasma integrated predictive modeling of tokamak plasmas.

Gyro-Bohm scaling and the L-mode shortfall
It was anticipated that the gyro-Bohm scaling of gyrokinetic turbulence would lead to small transport in the low temperature edge region of L-mode tokamak plasmas [21].The gyro-Bohm normalization factor for energy fluxes is where The choice of reference values [21] used for TGLF, GYRO [22] and CGYRO [23] are: electron density (n 0 = n e ), electron temperature (T 0 = T e ), ion mass (m 0 = m i ) and a magnetic field unit defined as (B 0 = B unit = q∂ψ r∂r ) and length (L 0 = a) where 'a' is the maximum of the minor radius 'r'.Here the subscript 'e' is for electrons and 'i' is for the main ion.With these units the gyro-Bohm energy flux scales like Hence, an energy flux of order the gyro-Bohm unit would drop sharply towards the edge of Ohmic or L-mode plasmas.This trend was indeed found with reduced transport models [2,21] leading to edge temperature gradients that were predicted to be much steeper than experiment in L-mode and Ohmic [10,11] plasmas.
An example of this problem is shown in figure 1 for DIII-D discharge 150 139 [2].The original saturation model (SAT0 red) in TGLF (labeled TGLF-09 in [2]) predicts a steep temperature gradient in the L-mode edge of DIII-D compared to the curve fit to the measured values (blue dashed).
Both electron (left) and ion (right) energy transport is underpredicted in the L-mode edge by TGLF-SAT0.An overprediction of the electron temperature, in inductive L-modes like this one, is common within the MHD interchange unstable region of the plasma (flat q-profile < 1.2).This occurs for r/a < 0.3 for this case.We will come back to this figure when we discuss the improved saturation models SAT1 and SAT2 in section 4.
The ITPA Integrated Operation Scenarios (IOS) group in a joint activity (TC-20) with the Transport and Confinement group conducted a test of the available reduced models in 2015 for predicting the plasma temperature profile evolution during the Ohmic current ramp up phase.Discharges with ITER-like current ramp-ups were run on three tokamaks (AUG, DIII-D, JET).The final report determined that none of the models were accurate enough to be used for simulating the control requirements of the ramp-up phase due to either the edge short-fall or the calibration to experiment required [24].The simulation of the DIII-D current ramp-up discharge with the FASTRAN code is shown in figure 2. Shown are the time evolution of the electron (left) and ion (right) temperatures at three different radial locations.Only the SAT0 model (red) was available in 2015 when the TC-20 test was conducted.The strong increase in electron temperature near the separatrix causes the predicted internal inductance to be too low by an unacceptable amount for control requirements.The SAT1 (green) and SAT2 (blue) models, described in section 4, improve the temperature prediction sufficiently for IOS modeling needs.
Proposed explanations for the cause of the transport missing from the quasi-linear models included the strongly nonlinear electromagnetic turbulence present in the open field line region [25].High collisionality instabilities, like the resistive ballooning modes [26], were proposed to explain the missing energy transport on the closed flux surfaces approaching the separatrix.A model of the 'drift resistive ballooning mode' transport was included in the Multi-Mode transport model [27].The model was able to increase the L-mode edge transport but it is not calibrated to gyrokinetic simulation.Electromagnetic gyrokinetic turbulence simulations, on closed flux surfaces, include all of the ion gyro-radius scale instabilities, including micro-tearing and resistive ballooning modes, so attention turned towards gyrokinetic simulations of the L-mode edge plasma.

Gyro-kinetic turbulence simulations
A shortfall in the predicted energy flux was seen in gyrokinetic simulations well away from the separatrix for some cases [5,8].The gyrokinetic turbulence fluxes from GYRO and GEM [28] are compared with TGLF-SAT0 in figure 3(left).The power balance flux from the ONETWO code [29], computed by integration of the sources, is shown in black.Both the electron (top left) and ion (bottom left) simulated energy fluxes fall below the power balance energy flux for r/a > 0.6, but the codes agree well with each other.The error bars are generated for TGLF by varying the temperature gradients  within experimental uncertainty.On the right side of figure 3 is a detailed comparisons of the predicted electron density (top right) and temperature (bottom right) fluctuation amplitudes from the GYRO simulations (synthetic diagnostic filtered) and measured on DIII-D [8,30].Similar results were obtained from analysis of C-MOD [9] L-modes, showing that both the energy fluxes and the fluctuation intensity were underpredicted in the L-mode edge.This strongly suggests that the transport is due to micro-turbulence but that the saturated intensity of the gyrokineitc simulations is too low.These gyrokinetic results appeared to be consistent with the quasi-linear modeling results and the expectations of gyro-Bohm scaling.These turbulence simulations were run with the GYRO code [22].The GYRO code, like other 'global' gyrokinetic codes, uses a grid in the radial coordinate so that it could include the spatial variation of the plasma gradients over most of the plasma radius.The radial grid requires a boundary condition that makes the turbulence damped at the boundaries.A Fourier transform to a radial wavenumber is required in order to evaluate the gyro-averaging operator (Bessel functions).However, a coarse radial grid evaluation of the Bessel functions, with a few radial wavenumbers, is usually used to save the computational expense of a complete Fourier transform of the radial simulation grid.
The first periodic flux tube simulations that were fully spectral in radial wavenumber, showed much better agreement with the power balance ion energy fluxes for L-modes in ASDEX Upgrade [12] and DIII-D [13] with the GENE code [31,32] and JT-60 with the GKV code [14].The GENE simulations also showed that the non-linear turbulence phase shift that drives energy and particle fluxes is determined by the phase shift of the most unstable linear mode.This is the requirement for the validity of the quasi-linear approximation.Hence, if a better model of the saturated fluctuation spectrum in the spectral flux tube simulations could be found, a reduced quasi-linear model of the transport would be possible.

Improved saturation models
With the successful ion energy flux matching by fully spectral flux tube gyrokinetic simulations, the ion channel L-mode shortfall problem appeared to be more of a problem of resolving the saturation of gyrokinetic turbulence rather than missing instabilities.However, the electron energy flux was still too low for some ion scale only simulations pointing to a role for the electron temperature gradient (ETG) instability.The ETG linear growthrate peaks at the electron gyroradius scale and only produces electron energy transport at this scale.Quasilinear models for the ETG modes were included in GLF23 [21], TGLF and Multi-Mode [33] but there was a lack of realistic multi-scale gyrokinetic turbulence simulations to determine the saturated fluctuation intensity of ETG modes.
The first real mass ratio simulations with kinetic electrons and ions became possible in 2015 [34].A strong suppression of ETG by ITG driven zonal flows was observed.A breakthrough set of multi-scale (ion and electron) simulations of a C-MOD L-mode with the physical masses for ions and electrons revealed enhanced electron scale turbulence for weak ITG modes [15,16].Analysis of the non-linear spectrum of potential fluctuations in these simulations lead to the following picture of how this multi-scale coupling comes about.The fluctuating ExB velocity couples modes with different poloidal (k y ) and radial (k x ) wavenumbers distributing the fluctuation intensity across the two dimensional (2D) spectrum leading to a Lorentzian flux surface averaged potential fluctuation distribution centered at k x = 0 in the absence of parity breaking effects.The zonal (k y = 0) fluctuations couple modes with the same k y across the k x spectrum.The high k x modes are damped by the gyro-averaging Bessel functions, and other effects, so the zonal fluctuations are a catalyst for the self-saturation of each k y slice of the spectrum.The zonal fluctuations are independent of k y so they have the same coupling to both ion and electron scale k y modes with a mixing rate of k y V ZF .The mixing velocity of the zonal fluctuations V ZF is computed from the zonal electric potential fluctuation spectrum ( φkx,ky=0 ) with the formula [18] V ZF = 0.5 The k x width of the zonal fluctuation spectrum is similar to the k x width of the ion-scale peak of the turbulence at finite k y suggesting the zonal fluctuations are driven by the ion scale turbulence.A detailed balance argument for the non-linear zonal mixing rate balancing the linear growthrate gives a saturation rule for the zonal mixing velocity given by where γ ky is the linear growthrate of the most unstable mode at each k y , which usually has k x = 0.This zonal fluctuation saturation rule for V ZF has been found to be well satisfied for a wide range of plasma parameters [20].Note that the zonal mixing rate k y V ZF can compete with the ETG linear growthrate but the zonal flow shear rate cannot [18] since it is independent of k y .It was observed in the multi-scale simulations that when the ETG linear growthrate exceeded the zonal mixing rate there was an enhancement of the saturated intensity in that region of the k y spectrum.This non-linear threshold γ ky > V ZF k y for ETG modes to overcome the zonal mixing has been successfully validated against experiments on C-MOD [35].This zonal mixing paradigm is the main feature of the multi-scale 2D (k x , k y ) saturation model SAT1 [18].The different saturation levels for global and local flux tube gyrokinetic turbulence simulations could have been due to the effects of profile shear (variations of the local plasma parameters) that are known to be stabilizing [36].If this were the case, then the global simulations would have been more accurate and the flux tube simulations simply closer to experiment because they neglected a suppression mechanism.Recent work has shown that the impact of profile curvature (diamagnetic velocity shear) in the diamagnetic drifts is negligibly small [37] consistent with it being an order ρ * 2 effect.
The SAT0 and SAT1 models were calibrated with a large database of non-linear GYRO simulations.These were flux tube runs with a low resolution radial grid and sparse sampling of the Bessel functions.These simulations have normalized energy fluxes that go down with increasing flux surface elongation.This appears reasonable since the linear growthrate decreases with elongation.However, the fully k x spectral CGYRO simulations have a normalized flux that increases with elongation.The different scaling with elongation was found to come from the change in the width of the radial wavenumber spectrum at the outboard midplane (θ = 0).The peak amplitude of the potential fluctuations scales like 1/k RMS x where k RMS x is the RMS width of the spectrum at fixed k y and θ = 0.The RMS width scales like k RMS x ≈ k y B(0)/(B unit |∇r| 0 ) = k y (rB/(qRB p ))| θ=0 for large enough k y [20].This is due to the fact that the perpendicular wavenumber in the Bessel functions is normalized to the ion gyro-radius with the full magnetic field B(θ) rather than the effective magnetic field B unit that is natural to the drift terms [38].Running GYRO in a periodic flux tube, with a full Fourier transform evaluation of the Bessel functions, recovers the spectral CGYRO result.Hence, it was the low resolution of the Bessel function evaluation that impacted the elongation dependence.The SAT2 model [19,20] is designed to fit the three dimensional (3D) spectrum (k x , k y , θ) of potential fluctuations from a database of high resolution CGYRO runs.Only the flux surface average of the model is used in the quasi-linear flux calculations.The linear eigenvalues and quasi-linear weights computed with CGYRO were used in the construction of SAT2.It achieves unprecedented accuracy (10% for energy fluxes) for a quasi-linear model when the CGYRO eigenmodes are used.
The magnetic field factor in the RMS width effectively changes the gyro-Bohm scaling.This can be incorporated into the normalization by using the outboard midplane magnetic field as the reference field, B 0 = B(0).This has a strong impact on the aspect ratio scaling of the normalization as well as flux surface shaping.The most natural reference length is the total pressure gradient length (|∇r| 0 /L p = −|∇r| 0 (∂p/∂r)/p) evaluated at the outboard midplane from the MHD force balance.Combining these gives a natural order parameter where ρ * CGYRO = T e /m i m i c/(eB unit a) and the ion thermal velocity has been used to compute the physical ion gyro-radius at the outboard midplane.
Using this natural unit to normalize the CGYRO fluxes in the SAT2 database yields fluxes that are of order 1 as expected from the ordering assumptions of gyrokinetic theory.As an example, the average energy flux ((Q e + Q i )/2) in the (regular) CGYRO system of normalization is plotted vs the case number for the database of CGYRO turbulence runs used to calibrate the SAT2 model in figure 4(left) and re-normalized using the natural units ρ * natural (right).The highest flux is for the case with R/r = 2.0 (case 43).The ratio of the fluxes from low (2) to high (12) values of R/r is reduced by a factor of 2 in the natural units due to the dependence of B(0)/B unit on R/r.Applying this natural gyro-Bohm normalization to a database of tens of thousands of DIII-D time and space points finds that the normalized power balance fluxes are also of order 1 consistent with weakly turbulent gyrokinetic transport.

Resolution of the L-mode shortfall problem
The SAT0 model [4] is 1D in k y , with a saturation at each value of k y that is independent of all the other k y modes, and depends on the linear growthrate for that mode.The SAT1 model [18] fits the 2D (k x , k y ) flux surface average saturated potential fluctuation spectrum.It includes coupling between different k y modes and to the zonal fluctuations.The SAT1 model can have larger electron transport due to ETG modes than SAT0 when the non-linear ETG threshold is exceeded.It can be seen that the shortfall for the DIII-D discharge 150 139 in figure 1 is greatly diminished for the SAT1 model (green).Even though ETG modes only produce electron energy flux, the collisional energy exchange impacts the ion power balance energy flux when the ratio of T i /T e is changed.The SAT2 model [19,20] is 3D (k x , k y , θ), fitting the poloidal angle (θ) variation of the saturated potential fluctuations intensity.This results in a flux surface geometry dependence that is close to the fully spectral gyrokinetic turbulence simulations.For the case in figure 1 the SAT2 model (blue) gives even more transport than SAT1 (green) near the separatrix.These same trends of higher thermal transport for SAT2 > SAT1 > SAT0 are seen in the simulations of the DIII-D current ramp discharge figure 2. In both of these discharges the temperatures are flatter near the separatrix for SAT2 than the DIII-D data indicating a need for continued refinement of SAT2.However, there is no sign of a shortfall in energy transport near the separatrix with SAT2.
A collection of DIII-D L-mode discharges was prepared by Kinsey et al to examine the parametric trends in the shortfall of the predicted TGLF-SAT0 energy fluxes compared to power balance [2].These discharges were carefully designed to be single parameter scans holding other dimensionless parameters fixed.The experiments collected cover scans in gyro-radius [39], collisionality [40], safety factor [41] and beta [42].This same set of DIII-D L-modes, plus a few more, was used to test the accuracy of the three saturation models.The fractional deviation (error ) 2 from the measured temperatures of the predicted ion and electron temperatures with TGLF from r/a = 0.98 to 0.5 (avoiding the sawtooth region (r/a < 0.5) is shown in figure 5.The data is organized into plasma current (top) and collisionality (bottom) scans.The errors in the electron (left) and ion (right) temperature predictions are shown.
For SAT0 (red squares) the shortfall increases as the plasma current is reduced or safety factor is increased.Uncertainty in the radiated power contribution to the power balance flux does not account for the shortfall.The shortfall is also worse at low collisionality.A clear trend of higher error at low current and low collisionality is observed.The trend of the TGLF-SAT0 predicted temperatures developing a steep gradient near the separatrix, as shown in figure 1, is repeated for all of these conditions, but is smaller for high current and high collisionality cases.
The SAT1 model (green circles), with higher ETG transport than SAT0, has lower errors in figure 5 than SAT0 for most cases.However, there is still a tendency of SAT1 to have larger errors at low current and/or low collisionality.The errors for SAT2 (blue triangles) in figure 5 show an improved error for SAT2 compared to the other two models.The trend of the error with plasma current and collisionality has been reduced for SAT2.Combining these results for DIII-D L-modes with the similarly good results for TGLF-SAT2 for AUG L-modes reported in [1] supports the conclusion that the SAT2 model resolves the L-mode edge energy transport shortfall.This has been accomplished making a more accurate fit to the 3D saturated electric potential fluctuation spectrum using the fully spectral CGYRO code.

Outlook for whole plasma integrated predictive modeling
Now that the L-mode edge shortfall problem has been resolved it is possible to model the current ramp and L-mode phase of a tokamak pulse as shown in figure 2. The inability to model the L-mode was an obstacle to predicting the L/H transition.This is the next problem that needs to be addressed.Once the Hmode transport barrier has formed, the pedestal stability and structure model EPED [43] can be used to predict the pedestal pressure and width.Coupling of the 1D pedestal (with FASTRAN) and 2D scrape-off layer (with SOLPS) has been demonstrated [44].This approach can now be extended to Lmode and Ohmic plasmas.A simple 2-point models for the boundary condition at the separatrix, allows fully predictive modeling of L-modes [1], and H-modes with a pedestal model [45], using TGLF-SAT2.Other than the L/H transition, the elements of a plasma pulse design simulator (PDS) capability [46], using validated theory based transport models, has now been demonstrated.
ITER would greatly benefit from a validated PDS to plan experimental discharges and ensure safe operation.Predictions of operational programming and contingencies will build confidence, and enable higher performance to be achieved in discharges with identified control and stability margins.There are still many open questions about the accuracy requirements of the PDS models and integrated modeling methods.The ITPA can help resolve the open physics gaps of whole plasma integrated modeling before ITER operations commence.This will require collaboration between the ITPA topical groups and the teams developing integrated modeling frameworks.Using the PDS tools to plan tokamak experiments now will quantify the accuracy of the predictions and establish the best practices for successful operations.This field testing can be done in parallel with the research to resolve outstanding physics/modeling gaps.The experience using a PDS, in the way it will be used on ITER, will guide the priorities for the model improvement research.Working together, we can deliver a robust, uncertainty quantified PDS tool for planning ITER operations.

Figure 3 .
Figure 3.Comparison of the electron (left top) and ion (left bottom) power balance energy flux (black) with gyrokinetic turbulence simulations from GYRO (blue dots) and GEM (green diamonds) [28] and the TGLF-SAT0 quasilinear model (red) and the measured (black) and simulated (light blue) electron density (right top) and electron temperature (right bottom) fluctuation amplitude for DIII-D L-mode 128 913 at 1500 ms.These figures are from an invited talk given by Chris Holland at the 2012 Transport Task Force Workshop.A published version of the power balance is in [8] and the measured fluctuations in [30].(Left) Reproduced from [8] © 2011 IAEA, Vienna.All rights reserved.

Figure 5 .
Figure 5. fractional deviation (sigma) of the predicted electron (left) and ion (right) temperatures for a set of DIII-D L-mode discharges scanning the plasma current (top) and collisionality (bottom).Note that the SAT2 cases were run with an early version of SAT2 and may differ from the published version.