Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas

We present gyrokinetic simulations with the GENE code addressing the near-edge region of an L-mode plasma in the DIII-D tokamak. At radial position $\rho=0.80$, simulations with the ion temperature gradient increased by $40\%$ above the nominal value give electron and ion heat fluxes that are in simultaneous agreement with the experiment. This gradient increase is consistent with the combined statistical and systematic uncertainty $\sigma$ of the Charge Exchange Recombination Spectroscopy (CER) measurements at the $1.6 \sigma$ level. Multi-scale simulations are carried out with realistic mass ratio and geometry for the first time in the near-edge. These multi-scale simulations suggest that the highly unstable ion temperature gradient (ITG) modes of the flux-matched ion-scale simulations suppress electron-scale transport, such that ion-scale simulations are sufficient at this location. At radial position $\rho=0.90$, nonlinear simulations show a hybrid state of ITG and trapped electron modes~(TEMs), which was not expected from linear simulations. The nonlinear simulations reproduce the total experimental heat flux with the inclusion of $\mathbf{E} \times \mathbf{B}$ shear effects and an increase in the electron temperature gradient by $\sim 23\%$. This gradient increase is compatible with the combined statistical and systematic uncertainty of the Thomson scattering data at the $1.3 \sigma$ level. These results are consistent with previous findings that gyrokinetic simulations are able to reproduce the experimental heat fluxes by varying input parameters close to their experimental uncertainties, pushing the validation frontier closer to the edge region.


I. INTRODUCTION
In order to improve the energy confinement time of magnetic fusion experiments, a thorough understanding of turbulent transport is necessary. The main carriers of cross-field transport in tokamaks are gyroradius-scale microinstabilities with scale lengths much smaller than the machine size. These microinstabilities are physically driven by electron and ion temperature gradients and density gradients. Since microturbulence is suspected to reduce energy confinement time with an increase in heating power, its mitigation is of intrinsic interest to the fusion community.
Paradoxically, a state of improved confinement arises together with steepened gradients in the edge region when the heating power is increased above a certain threshold power, P th . This transition from low confinement mode (Lmode) to high confinement mode (H-mode) was first discovered at the ASDEX tokamak in 1982 and has since been reproduced in all major tokamaks 1,2 . H-mode is the favored operational regime for nuclear fusion reactors and ITER. However, finding a self-consistent description of the L-H transition is a major unsolved problem [2][3][4] . An important first step towards understanding the L-H transition is to correctly describe Lmode plasmas in the near-edge and edge regions. This is also important for ITER, which will be in L-mode operation during plasma current ramp-up and ramp-down phases; correctly predicting the L-mode profiles is important for vertical stabilization of the plasma during these phases 5,6 . This motivates a study of microinstabilities in the L-mode near-edge just before an L-H transition.
Gyrokinetic theory provides an accurate description of mia) Electronic mail: tomneiser@physics.ucla.edu croturbulence in magnetically confined plasmas. Here, the assumptions of high background magnetic field, low frequencies relative to the ion cyclotron frequency and small fluctuation amplitudes typically apply 7 . Gyrokinetic codes such as GKV 8,9 , GEM 10,11 , GYRO 12 and GENE 13 have been in good agreement with the experiment and with each other in the core of both L-mode [14][15][16] and H-mode [17][18][19] plasmas. Similarly, the near-edge region of H-mode plasmas has been successfully modeled 17,18,20 . However, these codes have been in occasional disagreement in the near-edge of L-mode plasmas. For example, simulations have shown an underprediction of heat transport of ∼ 7 for GYRO 14 and ∼ 2 for GENE 16 for nominal parameters of DIII-D discharge #128913. This has raised fears of an apparent systematic shortfall of heat flux predictions 14,15 . However, these fears have recently been reduced 16,[21][22][23][24] . For example, increasing the ion temperature gradient within the experimental error bars has produced flux-matched simulations with GENE 16 . Recent GYRO simulations of a different L-mode discharge, namely DIII-D #101391, have revisited this shortfall problem and found good agreement with experiment 22 . Similarly, the CGYRO code 25 matches the experimental heat flux of this discharge 22 , as does the GENE code. Studies with GYRO on the Alcator C-Mod tokamak have also been in good agreement with experiment near the edge region of L-modes, with the use of multi-scale simulations in some cases 23,24 . These gyrokinetic validation exercises reduce fears that the shortfall is a universal feature of near-edge L-mode plasmas. However, these fears have not been completely removed yet. Moreover, the role of microturbulence and its interactions on multiple scales is still poorly understood in the near-edge. Therefore, there exists a continued need for code validation and microturbulence characterization in the near-edge of L-mode plasmas.
To address these issues, we present a gyrokinetic validation study of a DIII-D near-edge L-mode plasma just before an L-H transition with the gyrokinetic turbulence code GENE. Our primary finding is that gyrokinetic simulations are able to match the heat-flux in the near-edge b of the L-mode plasma at ρ = 0.8 and ρ = 0.9 within the uncertainty of the experiment at the 1.6σ and 1.3σ levels, respectively. In the course of this validation study, we make three secondary findings. First, current heuristic rules for the relevance of multi-scale effects appear to be on the cautious side; multi-scale simulations at ρ = 0.80 suggest that single-scale simulations can be sufficient in a scenario when multi-scale effects are expected, which could increase the realm of applicability of single-scale simulations. Second, the effect of edge E × B shear is found to be already important in the near-edge at ρ = 0.90, which was unexpected. Third, nonlinear simulations at ρ = 0.90 find a hybrid ion temperature gradient (ITG)/ trappend electron mode (TEM) scenario, which was not obvious from linear simulations. This successful validation exercise helps push the gyrokinetic validation frontier closer to the L-mode edge region.
This paper is structured as follows. We discuss the experimental data in section II and the gyrokinetic simulation method in section III. In section IV, we present linear and nonlinear simulation results at radial position ρ = 0.80. In section V, we present linear and nonlinear simulation results closer to the edge at ρ = 0.90. We reflect on the limitations and implications of our results in section VI, summarize the present work in section VII and outline future work in section VIII.

II. SUMMARY OF EXPERIMENTAL DATA
The subject of this study is DIII-D discharge #153624. This discharge exhibits three L-H transitions, where the first was induced with Electron Cyclotron Resonance Heating (ECRH) and the following two with Neutral Beam Injection (NBI) heating. Figure 1 shows highly resolved time-traces taken at ρ = 0.96 during the second L-H transition (since we focus our study on t 0 = 2400 ms, only data relevant to this transition is shown). Note the remarkable change in density fluctuation amplitudeñ/n and divertor D α recycling light emission as the plasma changes its operational state from Lmode to H-mode via an extended phase of Limit Cycle Oscillations (LCOs) 26,27 . The E × B velocity v E×B and density fluctuation amplitude were obtained with the Doppler Backscattering (DBS) diagnostic 28 , which probes a wavenumber range of 0.3 k y ρ s 0.6.
The plasma is in the lower single-null shape, where the ion ∇B-drift direction is towards the single active X-point and the power threshold for the L-H transition is relatively low (as compared to the upper single-null shape) 2 . Therefore, the neutral beam is deliberately operated at a relatively low beam b We loosely define the near-edge as 0.80 ≤ ρ ≤ 0.96, where ρ = (Φ/Φ edge ) 1/2 is the toroidal flux radius and Φ edge is the toroidal flux at the separatrix. power of 1.1 MW, which is in the vicinity of this power threshold. In steady-state operation, this heating power plus Ohmic heating is comparable to the total heat flux relevant for fluxmatching gyrokinetic simulations. Note that the heat flux is typically carried by radiation losses and by radial conductive and convective transport due to microturbulence in the electron and ion channels, which are sometimes difficult to separate empirically at moderate to high collisionality. This work will focus on a constant time t 0 = 2400 ms, because this is when the turbulence is still in a quasi-steady equilibrium state before the L-H transition begins (see Fig. 1). At t 0 the plasma has temperature and density profiles as shown below in Figure 2. Both the electron temperature and density are measured by the Thomson scattering diagnostic. This method produces highly resolved profile data for the electrons, but cannot be used effectively for the ions. This is due to their much lower Thomson scattering cross section σ t , which scales as σ t ∝ 1/m 2 j , where m j is the mass of the scattered charge. Therefore, a Charge Exchange Recombination Spectroscopy (CER) diagnostic is used for the impurity ions, which studies the emission lines from neutral beam injection. Specifically, the impurity ion CER diagnostic detects predominantly charge exchange recombination radiation from fully ionized impurity carbon ions (C 6+ ). The measured carbon temperature profile is assumed to be equal to the deuterium ion temperature pro- file in this discharge.
In general, the comparatively lower detection number statistics of CER data versus Thomson data causes the statistical uncertainty of the ion temperature to be larger than the uncertainty in the electron temperature (see Fig. 3). For instance, the statistical uncertainty for the ion temperature gradients at ρ = 0.80 is estimated to be σ ITG, stat ∼ 15% (see Appendix). On the other hand, the statistical uncertainty in the electron temperature gradient (ETG) at ρ = 0.90 is estimated to be σ ETG, stat ∼ 8%. Moreover, the systematic uncertainty for both electron and ion temperature gradients is estimated to be σ sys ∼ 10%. Recall that only slightly more than two thirds (68%) of normally distributed measurements fall within 1σ stat , while most (95%) fall within 2σ stat and nearly all (99.7%) fall within 3σ stat of the mean. Therefore, changes in temperature gradient of up to ≤ 2σ stat + 1σ sys can reasonably be attributed to combined uncertainties in the measurement and model assumptions.
The physical parameters at four radial positions in the nearedge and edge region are summarized in Table I. The variables in the table are defined as follows. The logarithmic gradients are defined as where ρ = (Φ/Φ edge ) 1/2 is the toroidal flux radius normalized by Φ edge , which is the toroidal flux at the separatrix. The shear parameter is given bŷ where q is the safety factor. The electron beta is defined as the ratio of thermal to magnetic pressure, where B ref is the magnetic field on axis and the other variables take their usual meaning. The effective atomic number of the plasma, Z eff = Σ i Z 2 i n i /n e , is greater than that of a pure deuterium plasma (Z = 1) mostly due to carbon impurities that enter the plasma from the divertor and the wall. The electron-ion collision frequency is defined as ν ei = Z eff n e e 4 ln Λ where ln Λ is the Coulomb logarithm and other variables take their usual meaning. Convenient reference lengths are given by For example, we will use ρ s to rescale the wavenumber of turbulent modes as k y ρ s . For a better comparison with other work, a conversion from k y ρ s to the toroidal mode number n is useful. The conversion factor is given by the relation 21 which is evaluated numerically in the bottom row of Table I. In this work, we also take into account the effect of the E × B shearing rate, which is a function of both the radial electric field gradient and the magnetic field geometry. For realistic geometry, the E × B shearing rate is given by the so-called Hahm-Burrell formalism 29 , where E r is the radial electric field, B θ is the poloidal component of the magnetic field B, and R is the major radius. Note that (E r /B θ R) is constant on the flux surface, but the factor cRB θ /B is not constant (it is larger on the outboard midplane). In the simulation community, a flux-surface-averaged formalism for the E × B shearing rate in shaped geometry is often used for simplicity 30 , where the factor preceding the derivative, r/q, is now a flux function. The present work will also employ this so-called Waltz-Miller formalism. The experimental shearing rate values are inferred using data from the magnetic field geometry and the Doppler Backscattering diagnostic (see Fig. 4).

III. GYROKINETIC SIMULATION METHOD
Throughout this work, we employ the gyrokinetic turbulence code GENE 13 . The gyrokinetic approximation reduces the six-dimensional phase space to five dimensions by averaging over a charged particle's gyro-motion perpendicular to the magnetic field, and removes several phenomena on small space-time scales. The resulting Vlasov equation can be coupled self-consistently to Maxwell's equations. GENE is a Eulerian code that solves the relevant equations on a field-aligned coordinate system (x, v , µ), which minimizes the necessary number of grid points 31,32 . Here, v represents the velocity along the field lines and µ is the magnetic moment resulting from the gyro-averaged motion of a charge.
GENE can model the plasma in a local (toroidal fluxtube 33 ) or a global (radial annulus) simulation domain. The local approximation is preferred where applicable, because periodic boundary conditions in both the radial (x) and binormal (y) directions invite numerically efficient spectral methods. The local approximation holds where the turbulent correlation lengths are smaller than the gradient scale lengths, such that the plasma parameters do not vary much across a typical turbulent structure. This condition appears to be satisfied in the near-edge region of the L-mode plasma considered here, so we will use the local approach throughout this work.
In our nonlinear simulations we consider two particle species (electrons and Deuterons), and thereby neglect the effect of carbon as motivated by linear simulations presented below in section IV and V. Electromagnetic effects are included in our simulations by solving for the parallel component of Ampère's law. Moreover, GyroLES techniques were used to model dissipation at unresolved wavenumbers, thus avoiding the unphysical build-up of energy at the highest resolved wavenumbers [34][35][36] . The flux-tube geometry is calculated with the TRACER-EFIT interface 37 . Due to the high collisionality in the near-edge region, a collision operator developed by Sugama et al. was used 38 .
Closer to the edge region, at ρ = 0.90, it becomes evident that the effects of sheared flow on turbulent heat transport need to be included, while the shearing rate at ρ = 0.80 is negligible (see Fig. 4). The radial flow shear can have an important effect, because it shears turbulent eddies in the poloidal direction, which increases their poloidal correlation length and reduces their radial correlation length [39][40][41] . This can lead to experimentally relevant improvements in particle and energy confinement. Specifically in GENE, the constant Waltz-Miller shearing rate throughout the flux-tube is implemented using a method developed by Hammett et al. 42 .
Here, a transformation into the co-moving coordinate system of the equilibrium flow and a discrete time evolution of the sheared radial wavenumber greatly reduce computational intensity while maintaining acceptable numerical accuracy 42 .
In order to investigate the interaction between strong ETG streamers and ion-scale modes, it is instructive to carry out multi-scale simulations resolving both ion and electron scales. These simulations are very computationally intensive and cannot currently be carried out resolving the full wavenumber domain of linearly unstable modes. Therefore, a reasonable reduction in the resolved wavenumber domain is sought with nonlinear single-scale simulations. For example, nonlinear electron-scale simulations with varying k y,max ρ s can help determine the reliability of GyroLES techniques; they help identify a reasonable maximum extent of the wavenumber domain that still captures the main nonlinear turbulent transport. A similar nonlinear scan in k y,min ρ s is carried out at the ion scales to determine a feasible multi-scale simulation domain that still captures the majority of the ion-scale physics. With this method, we are able to carry out multi-scale simulations with realistic electron-Deuteron mass ratio and realistic geometry for the first time in the near-edge, at ρ = 0.80. Throughout this work, the simulation domain in velocity space extends in the parallel direction up to v ,max = 3v th,j , where v th,j = 2T 0,j /m j is the thermal velocity. In the perpendicular direction, it extends to µ j = 9T 0,j /B ref . Further details of the simulation method, such as the number of grid points or the radial size of the simulation box, are described together with the results of the relevant nonlinear simulation. Throughout this work, the GyroLES methods, the Sugama collision operator and the inclusion of E × B shear effects mentioned above are of particular relevance.

IV. RESULTS AT FIRST RADIAL POSITION (ρ = 0.80)
All microinstabilities ultimately saturate due to nonlinear interactions between modes of differing wavenumbers. Nonetheless, linear simulations often give useful insight into the nature of the nonlinear instabilities. Generally, the distribution of linear growth rates in wavenumber space highlights the scales of the nonlinear turbulent drive. This can inform the size of the nonlinear simulation box. Moreover, sensitivity of these linear growth rates to changes in physical parameters such as temperature gradients can help identify the nature of modes encountered in nonlinear simulations. Therefore, we will first present linear simulation results, followed by nonlinear simulation results.

A. Linear Simulations
We use both the initial value solver and the eigenvalue solver in GENE to find linear growth rates of modes at the electron and ion scales. The resolution for all linear simulations in this work, which analyze each k y mode individually, is n x , n z , n v , n µ = (31,32,32,24). At ρ = 0.80, linearly unstable ion temperature gradient (ITG) modes can be identified by the positive frequencies associated with the diamagnetic drift direction of the ions. Similarly, trapped electron modes (TEM) and electron temperature gradient (ETG) modes are identified by the negative frequencies associated with the electron diamagnetic drift direction (see Fig. 5).
The frequencies of these modes are much smaller than the ion-cyclotron frequency such that the gyrokinetic approximation can be used to describe these modes. The unstable ITG and TEM/ETG modes are separated by a stable region in wavenumber space (0.68 < k y ρ s < 0.90). These separate domains allow us to clearly define ion scales (0.05 ≤ k y ρ s ≤ 0.80) and somewhat overlapping electron scales (0.70 ≤ k y ρ s ≤ 180) for separate nonlinear analysis. It is not clear how the heat fluxes found with nonlinear ion-scale and electron-scale simulations contribute to the collective heat flux. Generally, multi-scale simulations that simultaneously resolve both scale ranges are necessary to answer this question. These simulations are very expensive and may not be possible in all scenarios. A heuristic rule has emerged from pioneering work 43-46 using a reduced mass ratio ( m i /m e = 20) andŝ − α geometry, with α = 0. Namely, if the ratio of maximum growth rates at the electron and ion scales is much larger than the square root of the mass ratio, then multi-scale effects could be present. This is because the contributions by the electron-scale turbulence to the overall heat transport could be important. Otherwise, turbulent structures at the ion scales disrupt the efficient heat transport at the electron scales. In our case the mass ratio is m i /m e ≈ 60 and the ratio between the maximum growth rates is γ max ETG /γ max ITG = 338. Therefore our linear simulations with the nominal experimental parameters indicate that multiscale effects could be present. For an increase in the ion temperature gradient by ∼ 40%, we get γ max ETG /γ max ITG+40% = 161, so multi-scale effects could also be present at this point in parameter space, according to this rule-of-thumb. More recently, a model for saturation of multi-scale turbulence by zonal flow mixing has been proposed 47 . In this model, an important parameter is the RMS velocity of zonal flows (V ZF ), which saturates at V ZF = Max (γ/k y ), where γ is the linear growth rate of a turbulent mode with wavenumber k y . According to this model, multi-scale effects could be present when This criterion has been validated on recent multi-scale simulations using realistic mass ratio and geometry 23,24 . For our linear simulations at ρ = 0.80, this criterion is satisfied for the cases where the ion temperature gradient (ITG) is at its nominal value or increased by 10%; it is also satisfied when the ITG is increased by 40%, but by a small margin (see Fig. 6). There- For cases where the ITG is nominal (blue) or increased by 10% (green), the peak of said ratio is smaller at the ion scales than at the electron scales. For the case where the ITG is increased by 40% (orange), the difference between the peaks is much less pronounced.
fore, single-scale simulations are likely not sufficient for the cases considered here and multi-scale simulations will need to be carried out. This will be presented in the following subsection IV B. It is not a priori clear whether impurities significantly affect turbulence in the deuterium plasma. The DIII-D tokamak has a carbon divertor and wall that add carbon as the main impurity to the deuterium plasma. Using the impurity ion Charge Exchange Recombination (CER) diagnostic 48 , the plasma is observed to have an effective atomic number of Z eff = 1.80 (see Table I). In order to quantify the effect of this impurity, linear simulations are carried out with three particle species, namely deuterium and carbon ions, and electrons. We find that the carbon impurity has a negligible effect on the linear growth rates of ITG modes (see Fig. 7). Due to this observation, to first order in accuracy, carbon impurities can be neglected in our nonlinear simulations at ρ = 0.80. In summary, linear simulations at ρ = 0.80 identify coexisting ITG and ETG modes that could engage in multi-scale interactions. Linear growth rates and frequencies for a pure deuterium plasma (green) and for a plasma with an added carbon impurity species corresponding to Z eff = 1.8 (orange). For both cases, the ITG is increased by 10%, which is above the nonlinear critical gradient as shown in Figure 8. The carbon impurity has a small effect on the linear growth rates compared to the pure deuterium plasma.

B. Nonlinear Ion-scale Simulations
Fully nonlinear gyrokinetic simulations can be used to diagnose experiments in the hope to improve them in the future. When carrying out these simulations, several steps have to be taken to accurately extract the radial heat flux of the system. After each simulation, we ensure that the perpendicular box size (L x , L y ) accommodates several correlation lengths of the turbulence. This reduces the effect of boundary conditions on the turbulent structures and avoids their end-to-end connection across the boundaries. We check the grid resolution for convergence, (n x , n y , n z , n v , n µ ), by repeating a certain run with higher resolution in certain dimensions and checking for consistency with previous runs. This is particularly important closer to the edge, where high shear (ŝ > 2) demands high radial resolution 21 . Computational expense on the order of ten million central processing unit (CPU) hours (MCPUh) of multi-scale simulations presently restricts convergence tests to the single-scale simulation domain. To ensure an accurate reading of the simulated heat flux, it is averaged over a time greatly exceeding the turbulent correlation time.
Nonlinear ion-scale simulations are performed with (n x , n y , n z , n v , n µ ) = (256, 64, 24, 32, 24) grid points and perpendicular box size (L x , L y ) = (140ρ s , 126ρ s ). For the nominal experimental parameters as input we see a nonlinear quench of radial ion heat fluxes and the formation of a strong poloidal zonal flow (see lower inset, Fig. 8). Continuing the simulation for several hundred time units (∼ 400 L ref /c s ) to ensure nonlinear saturation, we find a time average ion heat flux that indicates nonlinear "stability" of ITG modes. In this case, the primary ITG instability leads to a secondary instability that generates a poloidal zonal flow and quenches the radial ion heat transport to Q i ∼ 10 kW (see lower inset of   8). Increasing the ion temperature gradient by ∼ 3.5% to ω Ti 2.9 marks the onset of nonlinear instability of ITG modes. This indicates a strong Dimits shift 49 , defined as the difference between the nonlinear and linear critical temperature gradients for onset of turbulent transport. Further increasing the ion temperature gradient leads to approximately linear increases in the electron and ion heat fluxes. These heat fluxes are carried by radially elongated turbulent structures (streamers) at the ion scales 33,50 (see upper inset of Fig. 8). This behavior is typical of plasmas in the core and has been used to infer the ion temperature gradient with gyrokinetic simulations 51 .
With the ion temperature gradient increased by 40%, our simulations recover the experimentally inferred heat fluxes for, remarkably, both the ion and electron heat channels. This ITG+40% scenario is consistent with the combined uncertainty of the ion temperature data at the 1.6σ level, where we write the combined uncertainty as σ = σ stat + σ sys (see Appendix). Specifically, the simulations give Q i = 610 kW and Q e = 190 kW, while the experimental values obtained with the ONETWO transport code are Q i = (600 ± 120) kW and Q e = (164 ± 33) kW (see Fig. 8).

C. Nonlinear Electron-scale Simulations
The nonlinear electron-scale simulations are performed with (n x , n y , n z , n v , n µ ) = (64, 512, 16,32,9) grid points and perpendicular box size (L x , L y ) = (9ρ s , 9ρ s ). Note that the box size is smaller relative to the ion-scale simulations because the electron-scale domain is defined by linear simulations as 0.70 ≤ k y ρ s ≤ 180. While finite Larmor radius (FLR) effects for k y ρ s 60 can be expected to significantly damp ETG modes, note that these high-k modes are mapped to smaller physical wavenumbers by geometric effects (encapsulated in the metric tensor). Thus, high-k modes can contribute meaningfully to the turbulent drive and it is advisable to extend the nonlinear electron-scale simulation domain over the entire wavenumber range of linearly unstable TEM/ETG modes. Note that for our electron-scale simulations we include fully kinetic ions (rather than only adiabatic ions).
For the nominal experimental input parameters we find a time average electron heat flux of Q e = 130 kW. Increasing the electron temperature gradient, which is the main driver of the electron heat flux, by its experimental error of ∼ 20%, we get a flux of Q e = 200 kW. This heat flux is within the experimentally inferred range of electron heat flux values obtained with the ONETWO transport code, namely Q e = 164 ± 33 kW (see Fig. 9). This suggests that electronscale heat transport could contribute to the overall heat transport, which will need to be studied with multi-scale simulations. Physically, we find that the electron heat flux is carried by radially elongated structures called streamers 52 . These structures are well-defined in a contour plot of electrostatic potential fluctuations,Φ(x, y) (see inset in Fig. 9).

D. Nonlinear Multi-scale Simulations
Linear and nonlinear simulations have indicated that multiscale interactions could be present. We have therefore carried out the first nonlinear gyrokinetic multi-scale simulations using a realistic mass ratio and experimental input parameters in the near-edge, at ρ = 0.80. These used on the order of 23 k processors and 10 MCPUh on NERSC supercomputers.
Resolving the full ion and electron scales is computationally prohibitive (for the full extent of the electron and ion scales, see section IV A). We therefore conducted a series of single-scale simulations with a sequentially reduced box size. This was done to identify an affordable domain that still resolves the main physical behavior of the plasma. For example, at the ion scales, we found that we could increase k y,min ρ s = 0.05 → 0.15 while maintaining the nonlinear heat flux to an accuracy of ∼ 10%. Similarly, at the electron scales, we were able to reduce k y,max ρ s = 180 → 40 with GyroLES techniques while maintaining a similar level of accuracy (∼ 15%) in the heat flux carried at the electron scales. This is significantly aided by the fact that most of the heat advection at the electron scales is carried by modes with k y ρ s ≈ 8. The flux-spectrum at the electron scales is plotted as the red dotted line in Figure 10. In this type of plot, adapted from Görler and Jenko [43][44][45][46] , the area under the curve roughly corresponds to the total heat flux carried by a range of wavenumbers. It is evident that, while the maximum growth rate in the linear flux spectrum is located at k y ρ s ≈ 50 (see Fig. 5), the nonlinear heat flux is carried predominantly by streamers associated with ETG modes with wavenumbers in the vicinity of k y ρ s ≈ 8 (see Fig. 10). This facilitates the above reduction in the resolved electron scales in the prepa-ration for multi-scale simulations. At the ion scales, the heat flux is carried predominantly by modes with k y ρ s > 0. 15. We thus resolve both the electron and ion scales in a carefully selected domain of poloidal wavenumbers of 0.15 ≤ k y ρ s ≤ 40 (see Fig. 10). Note that these nonlinear scans in simulation domain, while themselves computationally intensive, reduced the resource intensity of multi-scale simulations by a factor of 10, bringing them into the realm of the possible. The nonlinear multi-scale simulations are performed with resolution (n x , n y , n z , n v , n µ ) = (512, 512, 16,32,18) and perpendicular box size (L x , L y ) = (75ρ s , 42ρ s ). These simulations give the following qualitative results. First, we find that ETG-scale streamers co-exist with a zonal flow at ion scales when ITG modes are stable at nominal ω Ti . Second, we find that ETG-scale streamers are strongly sheared apart by ITG modes in the ITG+40% scenario (see Fig. 10). In this scenario, a brief coexistence of ETG-streamers and an ion-scale zonal flow can also be observed intermittently before both are disrupted by ITG-streamers (see Supplementary  Material). Therefore, electron-scale transport does not contribute significantly to the total transport when ITG modes are highly unstable, such as when ITG+40%. Thus, the heat-fluxmatching single-scale simulation in Fig. 8 is representative of the multi-scale heat flux.
In summary, simulations at ρ = 0.80 reproduce both the experimental ion and electron heat fluxes within the uncertainty of the CER data at the 1.6σ level. Multi-scale simulations suggest that turbulent structures on the ion scales strongly disrupt the streamers found on the electron scales. We now direct our attention further outwards to the nearedge region at ρ = 0.90. We first carry out linear simulations to identify the linear mode spectrum, which is dominated by TEM/ETG turbulence. Subsequent nonlinear simulations show high sensitivity of the total heat flux to changes in the electron temperature gradient. This could be due to a hybrid ITG/TEM scenario that was not predicted by linear simulations. With the inclusion of E × B shear and an increase in the electron temperature gradient by 23%, which is consistent with the experimental temperature data at the 1.3σ level, we are able to match the heat flux of the experiment. These results validate our gyrokinetic method and help push the gyrokinetic validation frontier closer to the edge region.

A. Linear Simulations
Linear simulations indicate that the main turbulent drive is carried by TEM/ETG modes at ρ = 0.90 (see Fig.'s 11 and  12). Note that this can be explained by the increased density gradient driving TEM turbulence closer to the plasma edge (see Fig. 2). Figure 12 shows the turbulent modes in lowwavenumber domain, which is most interesting for nonlinear simulations because the turbulent advection is most efficient at these large scales. To identify the subdominant modes that can play a role in nonlinear simulations, we employ the eigenvalue solver in GENE. Curiously, the subdominant mode is an ITG-type mode that is stable over all considered wavenumbers, which can be seen by its negative growth rates (see Fig. 12). Moreover, we find that an increase in electron temperature gradient by 23% further destabilizes the TEM/ETG modes. These linear simulations indicate that the TEM/ETG modes are likely to dominate in nonlinear simulations, with negligible effect of ITG modes. At ρ = 0.90, the carbon impurity has a small effect on the linear growth rates compared to the pure deuterium plasma (see Figure 13). Therefore, to reduce computational complex- 13. Linear growth rates and frequencies for a pure deuterium plasma (red) and for a plasma with Z eff = 1.8 due to a carbon impurity species (orange). The effect of carbon impurities on growth rates is small and will be neglected in nonlinear simulations.
ity by approximately 30%, we carry out nonlinear simulations at ρ = 0.90 with two rather than three particle species.

B. Nonlinear Ion-scale Simulations
Our linear simulations have shown that the growth rates at ρ = 0.90 are sensitive to changes in ω Te (see Fig. 12). This is due to the predominance of TEM/ETG-type modes with ITG modes subdominant and stable at nominal gradients. We thus study the sensitivity of nonlinear simulations to changes in ω Te .
The nonlinear ion-scale simulations are performed with resolution (n x , n y , n z , n v , n µ ) = (512, 64, 32, 32, 18) and perpendicular box size (L x , L y ) = (188ρ s , 126ρ s ). Note that convergence tests found that a higher resolution was required for the ion scales at ρ = 0.90 than at ρ = 0.80. Physically, this is due to the need to resolve higher magnetic shear here (see Table I). Moreover, the radial box size was increased because the simulated plasma was more susceptible to simulation boundary effects at ρ = 0.90 than at ρ = 0.80.
We define the ion-scale domain as 0.05 ≤ k y ρ s ≤ 1.60 and employ GyroLES techniques to avoid the unphysical build-up of free energy at k y ρ s 1.60. The individual ion and electron heat channels are difficult to distinguish experimentally with current techniques due to the high collisionality at ρ = 0.90, so that the observable here is the total heat flux.
Nonlinear simulations with an increase in ω Te by up to ∼ 30% were performed. Without the inclusion of experimental E × B shear, we found a saturated time-average flux of up to Q tot = 5 MW (see Fig. 14). Interestingly, we see high sensitivity to increases in ω Te between the +22% and +23% mark. The inset shows a contour plot of electrostatic potential fluctuations with high levels of turbulence at large scales for the +23% scenario without E × B shear. The observations FIG. 14. Total heat flux as a function of the electron temperature gradient (ETG) in nonlinear simulations at ρ = 0.90. Small changes in ETG from +22% to +23% show large changes in the heat flux, with high levels of turbulence for the +23% case (see inset). The sensitivity to the electron temperature gradient could be caused by a hybrid TEM/ITG scenario investigated in the following figures. The sensitivity is tempered by the experimental E × B shearing rate.
persisted with an increase in the radial box size and an independent numerical test of the validity of GyroLES techniques. This indicates that there may be a physical origin for this high sensitivity of the heat flux to changes in ETG.
The following analysis shows that TEM turbulence may nonlinearly destabilize the linearly stable ITG modes, leading to a hybrid ITG/TEM scenario in our nonlinear simulations. Recall that the ITG modes are stable for all wavenumbers (see Fig. 12). Moreover, the ETG/TEM branch with ETG increased by 23% (ETG+23%) is linearly unstable for k y ρ s > 0.18. Based on these observations, we expected nonlinear simulations to show that the majority of the heat flux is carried by electrons in the vicinity of k y ρ s ∼ 0.30. However, the heat-flux spectrum of nonlinear simulations with ETG+23% shows that the majority of the heat flux is carried by ions in the vicinity of k y ρ s ∼ 0.15 (see Fig. 15). This indicates high ITG mode activity, which was not predicted by linear simulations. To test the suspected ITG-dependence of the heat flux, we reduced the ion temperature gradient by 100% and uncovered the purely TEM/ETG turbulence with electrondominated heat flux carried predominantly at k y ρ s ∼ 0.30 as originally expected.
To further identify the nature of nonlinear turbulence at ρ = 0.90, we have plotted the frequency spectrum of electrostatic potential fluctuations (see Fig. 16). We find that the modes that carry most of the heat flux (i.e. k y ρ s ∼ 0.15) have positive frequency and therefore are associated with ITG turbulence. Moreover, the nonlinear frequencies at k y ρ s = 0.15 and k y ρ s = 0.30 are in good agreement with linear frequencies. This indicates that both TEM and ITG turbulence are active in nonlinear simulations. This seems to suggest that there exists a hybrid ITG/TEM scenario with nominal ITG and ETG+23% (even though linear simulations did not pre-   dict this). We now study the cross-phases between fluctuations in the electrostatic potentialΦ and the perpendicular electron tem-peratureT ⊥,e (see Fig. 17). For the wavenumbers carrying most of the heat flux, k y ρ s ∼ 0.15, the nonlinear cross-phases agree with the linear cross-phases of the ITG modes. Moreover, at k y ρ s ∼ 0.15 the nonlinear cross-phases are out of phase by approximately π/2, which is characteristic of sig- FIG. 17. Cross-phases between fluctuations in the electrostatic potential and perpendicular electron temperature for the scenario with nominal ITG and ETG+23%. At wavenumbers responsible for heat transport, kyρs ∼ 0.15 (horizontal line), the nonlinear cross-phases agree with those of the linear ITG modes (blue circles). For higher wavenumbers, the nonlinear cross-phases continue to agree with linear cross-phases associated with ITG modes rather than the TEM modes (red triangles).
nificant electrostatic heat transport. For higher wavenumbers, the nonlinear cross-phases agree with the linear cross-phases of ITG modes rather than TEM/ETG modes. A similar hybrid ITG/TEM scenario was previously found in Ref. 21, with the difference that both the dominant TEM and subdominant ITG modes were linearly unstable. Our results indicate a hybrid ITG/TEM scenario, with linearly subdominant and stable ITG modes carrying most of the heat flux. This scenario could also be relevant to spherical tokamaks, where ITG modes are more often linearly stable than in conventional tokamaks due to the lower aspect ratio 53 .
One possible explanation for this scenario could be that the turbulent fluctuations of TEM turbulence nonlinearly excite the linearly stable ITG modes. This is reasonable by elimination, as there is insufficient linear drive and no alternative nonlinear drive for the ITG modes other than the TEM turbulence. Therefore, it appears that ETG+23% is the critical gradient not only for TEM turbulence, but for hybrid ITG/TEM turbulence. This could explain the high heat-flux stiffness in our simulations at ρ = 0.90 (see Fig. 14).
We now introduce electric field shear due to the L-mode E r -well present in the edge region of most tokamaks. Generally, E r -wells are much more pronounced in H-mode configurations, but are also present in L-mode plasmas 2 (see Fig. 4). We find that a value of ω E×B = 0.5 c s /L ref , which is in the middle of the experimentally inferred range, is able to reduce the total heat flux approximately to the experimentally inferred values (see Fig. 14). Therefore, the E × B shearing rate may be an important simulation parameter to accurately model the heat flux in the edge region of L-mode plasmas. We conclude that we are able to reproduce the total experimental heat flux at ρ = 0.90 with an increase of ω Te by ∼ 23%. This increase is within the combined uncertainty of the Thomson scattering data at the 1.3σ level (see Fig. 3).
Lastly, multi-scale effects may be present at ρ = 0.90, but their investigation is computationally too expensive to fit within the scope of this work. Based on our experience with highly unstable ITG modes at ρ = 0.80, we predict that the effect of electron-scale streamers would likely be strongly reduced by the highly unstable ITG/TEM turbulence found at the ion scales at ρ = 0.90.

VI. DISCUSSION
To address an apparent shortfall problem 14,15 , recent validation exercises studying L-mode plasmas in Alcator C-Mod, ASDEX Upgrade, and DIII-D tokamaks have reduced fears that the shortfall is a universal feature of near-edge L-mode plasmas 16,21,23,24 (see section I). The results presented here are consistent with previous work that has been able to reproduce the experimental heat flux by changing input parameters close to their experimental uncertainty 16,21 .
Results from multi-scale simulations with realistic Deuteron-electron mass ratio and geometry at ρ = 0.80 have been presented. An early heuristic rule, found with pioneering multi-scale simulations [43][44][45][46] , has suggested that ETG-modes can contribute experimentally relevant heat flux if γ max ETG /γ max ITG m i /m e . This rule of thumb was found using a reduced mass ratio and simplified geometry in the core, and is not expected to apply universally. Nevertheless, this rule appears to hold in recent multi-scale simulations using more realistic parameters with the GKV code 8,9 and GYRO 54,55 . In the present work, an example of the limit of applicability of this rule may have been found. Recall that linear simulations give γ max ETG γ max ITG+40% = 161 (see section IV A). However, multi-scale simulations have qualitatively found very little ETG contribution to the overall heat flux with an increase in ITG by 40% (see section IV B). Physically, largescale turbulent structures of ITG-modes are able to shear ETG streamers apart. Thus, when ITG modes are highly unstable, they strongly reduce the flux carried by high-k modes. More recently, the condition Max (γ ETG /k y ) ≥ Max (γ ITG /k y ) has been identified as a simple test to predict possible multi-scale effects 47 . This test has been validated with recent multi-scale simulations and seems to be less conservative than using only the ratio of maximum growth rates. For our linear simulations at ρ = 0.80, this linear condition is only marginally satisfied for the ITG+40% scenario; namely, we get Max (γ ETG /k y ) /Max (γ ITG /k y ) = 1.16 > 1 (see Fig. 6). Since this linear test is intended as a heuristic rule and not a hard-and-fast rule, our results of negligible multi-scale effects are not significantly at odds with this linear test. We therefore conclude that the recently proposed test appears to be more useful than relying only the ratio of maximum linear growth rates to instruct whether multi-scale simulations might be necessary.
There are several limitations associated with our study. First, it became necessary to include the E × B shear for nonlinear simulations at ρ = 0.90. As mentioned in sec-tion III, the realistic shear in real geometry is given by the Hahm-Burrell formalism, which is not a flux function. For the purposes of our study, we assume that the flux-surfaceaverage shear given by the Waltz-Miller formalism is representative of the total shear effect. This is a common assumption in the simulation community 30 . The more accurate Hahm-Burrell formalism has not been implemented in GENE yet. Nevertheless, our simplifying assumption does not affect the generality of our finding that E × B shear from the L-mode E r -well is already important for simulations at ρ = 0.90. Second, we have carried out multi-scale simulations in the near-edge for the first time. These simulations would not be possible without simplifying assumptions. For instance, we have relied on GyroLES techniques to replace unresolved dissipation with a model. Moreover, we have relied on a large number of single-scale simulations to empirically determine the validity of our multi-scale simulation parameters. However, multi-scale convergence tests may be helpful to test the sensitivity of multi-scale effects to simulation parameter changes. Since these are computationally expensive they are beyond the scope of the present work. Third, note that particularly the ion temperature data used for this work has a large uncertainty associated with it, likely due to low carbon density (see Fig. 3). We have estimated the relevant statistical uncertainties for the flux-matching temperature gradients to be σ ITG, stat ∼ 15% and σ ETG, stat ∼ 8% (see Appendix). We have estimated reasonable systematic uncertainties as σ sys ∼ 10%. As a result of these error estimates, we find that our flux-matching gradients of ITG+40% and ETG+23% both fall within ≤ 2σ stat + 1σ sys . It is worth repeating here that only slightly more than two thirds (68%) of normally distributed measurements fall within 1σ stat , while most (95%) fall within 2σ stat and nearly all (99.7%) fall within 3σ stat . Therefore, our deviations from the experimental measurements of ≤ 2σ stat + 1σ sys are acceptable from a purely statistical perspective and sufficient for gyrokinetic validation. Nonetheless, there may be limitations in our validation method that exclude some relevant physics. For example, we are only considering a local flux-tube domain for our simulations and we are assuming a purely Maxwellian background distribution. While we have ruled out multi-scale effects, using a global model or including non-Maxwellian fastions could contribute sufficient missing physics to affect our conclusions. Including these phenomena is beyond the scope of the present work.

VII. SUMMARY
We have presented results from a study of a DIII-D L-mode plasma in the near-edge. At ρ = 0.80, the radial ion flux is quenched by strong poloidal zonal flows for nominal input parameters. In the ITG+40% scenario, nonlinear single-scale simulations give remarkable agreement with both the ion and electron heat fluxes of the experiment. This change in gradient is compatible with the combined statistical and systematic uncertainty in the ion temperature gradient at the 1.6σ level (see Fig. 3). At the electron-scales, radially elongated streamers are found to carry significant electron heat flux that is comparable to the experiment. This motivates multi-scale simulations, which were carried out for the first time in the near-edge with realistic mass ratio and geometry. Results suggest that the highly unstable ITG modes in the flux-matched ion-scale simulations strongly suppress turbulent transport at the electron-scales. Therefore, single-scale simulations are sufficient to match the experimentally inferred heat flux by changing the ion temperature gradient within the uncertainty of the experiment at ρ = 0.80. At ρ = 0.90, nonlinear simulations uncover a hybrid ITG/TEM scenario, which was not predicted by linear simulations. Moreover, nonlinear simulations are able to match the total experimental heat flux in the ETG+23% scenario when E × B flow shear (as evaluated from Doppler Backscattering measurements) is taken into account. This is consistent with the combined uncertainty of the electron temperature measurements at the 1.3σ level. Therefore, our primary conclusion is that gyrokinetic simulations are able to match the heat-flux in the near-edge of the L-mode plasma by varying input parameters close to their experimental uncertainties at ρ = 0.8 and ρ = 0.9.

VIII. FUTURE WORK
The present work invites several avenues for future research. For instance, future work could study the nature of multi-scale effects for the ITG+3.5% scenario: When ITG modes are marginally unstable, previous multi-scale simulations in the core have found that (i) ETG streamers can contribute experimentally significant transport at small scales 23,24 and (ii) ETG streamers can dampen poloidal zonal flows and enhance ion-scale transport 55 . Similar effects might be found with our simulations in the near-edge when ITG modes are brought close to marginal stability. Building upon the results presented here, future work could quantify the effect of multiscale interactions in the near-edge of L-mode plasmas.
The main ion temperature profile is often assumed to be equal to the impurity ion temperature measured with the CER diagnostic (see section II). This could introduce systematic errors. However, recent diagnostic development at DIII-D is now able to extract the main ion temperature directly and is currently quantifying this known source of uncertainty 56 . A main ion CER (MICER) diagnostic has been developed to study the deuterium ion (D + ) charge exchange signal. In order to test the conclusions of this work, future work could study a similar L-mode discharge using data from the more precise and more accurate MICER diagnostic currently under development at DIII-D 56 .
Closer to the edge region, global simulations are likely required, as the shearing rate changes substantially within a narrow radial region just inside the last closed flux surface. Resistive ballooning modes can also be expected to potentially contribute to thermal edge transport [57][58][59] . Future work could extend the present gyrokinetic simulations closer to the edge of plasmas just before the L-H transition 26,27 .