4.1 Local, linear scans

To characterise the instabilities present in the given pedestal we have performed scans with linear local simulations in the binormal wavenumber $k_y$ from $0.05{\rho _i}^{-1}$ to $350{\rho _i}^{-1}$ (corresponding to e.g. $N_{\rm tor} =[4,29954]$ at $\rho _{\rm tor}=0.86$ or $N_{\rm tor} = [7,48365]$ at $\rho _{\rm tor}=0.99$), the radial wavenumber at the outboard midplane $k_{x,\text {centre}}$ from $-$40 to 40 (connected to the ballooning angle $\theta _0=k_{x,\text {centre}}/(\hat {s}k_y)$) and radial position $\rho _{\rm tor}$ from 0.85 to 1. Figure 4 shows the growth rate spectra at four positions (pedestal shoulder, intermediate, pedestal centre and pedestal foot) maximised for each $\rho _{\rm tor}$ and $k_y$ over the ballooning angle. At selected wavenumbers, ballooning angles and radial positions we have additionally performed scans in temperature and density gradients ($\pm 30\,\%$), collisionality (0–18$\nu ^*_e$) and plasma $\beta$ (0 %–1 %). In total about 7000 simulations were performed for the characterisation of linear instabilities.

Figure 4. Growth rate spectra for four radial positions. Here $c_s=\sqrt {T_e/m_i}$ and $\rho _i$ are local normalisations evaluated at each radial position, and $L_{\rm ref}=0.65$ m. Shaded regions indicate the wavenumber ranges used in nonlinear heat flux simulations.

We use the following criteria to distinguish between the different instabilities: parity of the parallel mode structure in ballooning representation (tearing or ballooning), size (ion or electron gyroradius), drift direction (ion or electron diamagnetic), sensitivity on gradients ($T_e$, $T_i$, $n$), dependence on collisionality and plasma $\beta$ as well as diffusivity and heat flux ratios following the fingerprint approach (Kotschenreuther et al. Reference Kotschenreuther, Liu, Hatch, Mahajan, Zheng, Diallo, Groebner, Hillesheim, Maggi, Giroud, Koechl, Parail, Saarelma, Solano and Chankin2019). An alternative way to identify tearing modes is the tearing or parity factor (see Patel et al. Reference Patel, Dickinson, Roach and Wilson2022), which has been used to verify the identifications. An electromagnetic mode on scales ${>}\rho _i$ with tearing parity is called MTM. An electromagnetic mode on scales ${>}\rho _i$ with ballooning parity and drift in ion diamagnetic direction is called KBM. An electrostatic mode on scales ${\approx }\rho _i$, which is stabilised by collisionality and is not destabilised by $\boldsymbol {\nabla } T_i$, is called TEM. An electrostatic mode on scales ${\approx }\rho _i$, which is destabilised by $\boldsymbol {\nabla } T_i$ and propagates in ion diamagnetic direction, is called ITG. An electrostatic mode on scales ${\lessapprox }\rho _i$ that is destabilised by $\boldsymbol {\nabla } T_e$ is called ETG. However, it should be emphasised that a clear categorisation of modes on ion scales in the region of very strong drive close to the separatrix is particularly challenging. Different drive mechanisms interact and excite instabilities with characteristics that do not fall neatly in the mode prototypes developed in the study of core turbulence.

At the pedestal shoulder ($\rho _{\rm tor}=0.86$; violet triangles in figure 4) we find MTMs on largest scales, TEMs at ion scales, a region of stable wavenumbers and then ETGs on electron scales. In agreement with Kotschenreuther et al. (Reference Kotschenreuther, Liu, Hatch, Mahajan, Zheng, Diallo, Groebner, Hillesheim, Maggi, Giroud, Koechl, Parail, Saarelma, Solano and Chankin2019), the linear TEMs produce significant particle transport ($D_e/(\chi _i+\chi _e)\approx 0.15$, with particle diffusivity $D$, heat diffusivity $\chi$) compared with MTMs. At the intermediate position ($\rho _{\rm tor}=0.94$; blue squares in figure 4), where the magnetic shear has a minimum, a growth rate gap without any clear mode on ion scales exists. On intermediate scales, more ETG modes occur, which were not present at the pedestal shoulder. In the steep gradient region ($\rho _{\rm tor}=0.97$; green crosses in figure 4) on ion scales we find ITG modes that are close to a KBM transition. At smaller scales, ETG-driven modes are present that extend to larger scales towards the pedestal centre and foot. At the pedestal foot ($\rho _{\rm tor}=0.99$; yellow circles in figure 4) on ion scales modes show an ETG/TEM character but tend to be destabilised with increasing collision frequency. The ETG modes at pedestal centre and foot tend to peak increasingly at finite ballooning angle, indicating the presence of toroidal ETG modes (Told et al. Reference Told, Jenko, Xanthopoulos, Horton and Wolfrum2008; Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba, Ball and Ivanov2020). While the fastest-growing ETG modes we find tend to be toroidal in nature, subdominant slab ETGs turn out to play an important role in the nonlinear simulations (see § 5). Overall growth rates increase from pedestal top to foot.

Scans over plasma $\beta$ at different radial positions show that the pedestal in these linear local simulations sits close to a KBM threshold – being closer at the pedestal foot than at the pedestal top (cf. figure 5). The closeness to the KBM threshold in the pedestal centre resembles the use of a KBM constraint in the EPED model (Snyder et al. Reference Snyder, Groebner, Leonard, Osborne and Wilson2009) for the prediction of the pedestal width. It has, however, been reported that the radial structure of KBMs may not be compatible with the narrow pedestal region (Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Merlo, Field, Giroud, Hillesheim, Maggi, Perez von Thun, Roach and Saarelma2019; Predebon et al. Reference Predebon, Hatch, Frassinetti, Horvath, Saarelma, Chapman-Oplopoiou, Görler and Maggi2023), suggesting that the details of the KBM threshold may be skewed by the local approximation.

Figure 5. Growth rate scan over plasma $\beta$ at $k_y\rho _i=0.075$ and $k_{x,\text {centre}}=0$ at four radial positions. Nominal $\beta$ value indicated in orange.

4.2 Pressure and magnetic shear effect: second stability region

The low growth rates on ion scales in the intermediate region ($\rho _{\rm tor}=0.94$) between pedestal top and steep gradient are caused by the interplay of low magnetic shear and already increased pressure gradient in this region. Their interplay locally stabilises ballooning modes (second stability region). This can be shown by a scan over the magnetic shear at this position (see figure 6b). At nominal parameters (black crosses) no clear ion-scale mode is present, but with increasing magnetic shear $\hat {s}$ (green triangles and yellow circles) the system leaves the second stability region and an ion-scale mode becomes unstable. At the pedestal shoulder, where a lot of transport is driven on ion scales, lowering the magnetic shear, unfortunately, does not decrease growth rates significantly, since the pressure gradient is too low to access the second stability region (figure 6a, green triangles). In contrast, a high magnetic shear lets the pedestal shoulder enter the first stability region, where TEMs are suppressed (yellow circles). Further scans with, for example, modified profiles are not within the scope of this paper, but these results illustrate the potential of reducing micro-turbulence instabilities through careful tailoring of safety factor and pressure profiles.

Figure 6. Linear, local growth rates as a function of the toroidal mode number $N_{\rm tor}$ for three magnetic shear ($\hat {s}$) values at two different radial positions:(a) pedestal shoulder $\rho _{\rm tor}=0.86$ and (b) low-shear region $\rho _{\rm tor}=0.94$. Black crosses are nominal growth rates.

In the next section, we analyse the nonlinear turbulent system that is driven by the presented instabilities and compare to what extent linear mode signatures prevail in the nonlinear state. To reduce computational cost, ion scales and electron scales are treated separately.

5.1 Global, ion-scale simulations

To calculate heat fluxes we have performed gradient-driven, global, nonlinear, electromagnetic simulations on ion scales at experimental $\beta$ values, which have been enabled by the code upgrade presented in § 2. These simulations cover $\rho _{\rm tor} = 0.85$–$0.995$ in radius (almost the full width shown in figure 1) and $k_y\rho _i = 0.05\unicode{x2013}1.6$ in binormal wavenumber (corresponding to a toroidal wavenumber $N_{\rm tor} \approx 4\unicode{x2013}124$), indicated by the left shaded region in figure 4. They are two-species simulations (deuterium and electrons) with correct mass ratio ($m_e/m_D=1/3670$), collisions (Landau collision operator) and perpendicular magnetic fluctuations $\bar {A}_{1\parallel }$, but without compressional magnetic perturbations $B_{1,\parallel }$. When indicated, the simulations include the effect of $E\times B$ rotation due to the radial electric field $E_r$ (cf. figure 2). For numerical reasons the background flow in GENE is restricted to the toroidal direction. As an approximation to the $E\times B$ shear effect we retain the magnitude of $v_{E\times B}$, but rotate it to be purely toroidal. In total, the nonlinear global simulations required about $2.5 \times 10^6$ CPUh. More details on, for example, grids are specified in Appendix A.

The obtained heat fluxes are shown in figure 7 as a function of time averaged over real space (flux surface and radius, including the radial buffer zones; see Appendix A). In both ion and electron channels, the electrostatic heat flux dominates. But while the ions show vanishing electromagnetic heat flux, the electron heat flux is about 1/4 electromagnetic. When including $E\times B$ shear in the simulations the heat flux is strongly damped in all channels: the electrostatic heat flux by a factor of three and the electromagnetic heat flux even more strongly.

Figure 7. Heat fluxes (electrostatic $Q_{\rm es}$ and electromagnetic $Q_{\rm em}$) for electrons (blue, left) and ions (red, right) as a function of time. The green continuations are performed including an external background velocity shear corresponding to experimentally measured $E\times B$ shear.

5.2 Connecting linear and nonlinear results: frequencies and cross-phases

To identify the dominant turbulent transport mechanisms, results from linear and nonlinear simulations have to be connected, to test if the linearly fastest-growing modes remain important in the nonlinear saturated state. This can be achieved by comparing frequencies and cross-phases.

Figure 8 shows the mode frequencies of the linear and nonlinear simulations. Figure 8(a) shows the frequency spectrum of the local, linear simulations for three radial positions (pedestal shoulder at $\rho _{\rm tor}=0.86$, centre at $\rho _{\rm tor}=0.97$ and foot at $\rho _{\rm tor}=0.99$), in the co-moving frame of reference. In the pedestal shoulder spectrum (violet triangles) the transition from MTM to TEM is visible and at the pedestal centre (yellow circles) the transition from ITG to TEM/ETG. Figure 8(b) shows the frequency spectrum of the global, nonlinear simulation analysed at two positions (pedestal shoulder and centre) overlaid by the local, linear spectra. Usually frequencies from local simulations are specified in a frame co-moving with the $E\times B$ rotation. Since $E\times B$ rotation depends on radial position, frequencies of global simulations are specified in the laboratory frame and for the comparison of both we transform the local frequencies also to the laboratory frame. The comparison shows that at the pedestal shoulder and even at the centre the linear frequencies persist in the nonlinear simulations. This indicates that the linearly fastest-growing modes remain important in the saturated turbulent state. However, one important difference appears. At the pedestal shoulder ($\rho _{\rm tor}=0.86$) and $k_y\rho _i\approx 0.3$ the nonlinear simulation (figure 8b, upper red line) does not show the mode transition that linear, local simulations present (figure 8b, violet triangles). This suggests that MTMs are suppressed or at least restricted to the very largest scales in global, nonlinear simulations compared with the local, linear ones.

Figure 8. Frequency spectrum on ion scales for linear, local simulations in the co-moving frame (a) and global, nonlinear results compared with the linear, local scans in the laboratory system (b). The blue colour map (b) shows the distribution of frequencies at two positions from the nonlinear simulations. The red lines are the mean of these distributions. The dashed line in (a) indicates the transition between ion and electron diamagnetic drift direction.

The cross-phases of electric potential and electron density fluctuations (see figure 9) support this picture. Linear mode structures survive in the nonlinear simulations, but on largest scales at the pedestal shoulder differences are visible, corroborating the suppression of MTMs in global, nonlinear simulations observed in the frequency comparison. Overall, the remarkable agreement seen in the frequency and cross-phase comparison between local/linear and global/nonlinear simulations at pedestal shoulder and centre encourages the extension of quasi-linear models to the pedestal region.

Figure 9. Cross-phases of electric potential $\phi$ and electron density fluctuations $n$ from nonlinear simulations (blue background) and linear simulations (orange circles).

5.3 Local, electron-scale simulations

Local, nonlinear electron-scale simulations have been performed at multiple radial positions between pedestal top and foot. They complement the global, nonlinear simulations, which only cover ion scales up to $k_y\rho _i = 1.6$ at the reference flux surface. The dedicated ETG simulations cover $k_y\rho _i = 2.5\unicode{x2013}157.5$ (i.e. $k_y\rho _e = 0.04\unicode{x2013}2.6$ as indicated by the right shaded region in figure 4). The simulations use adiabatic ions. Resolutions and further simulation settings are given in Appendix A. Figure 10 shows the ETG heat flux obtained in the steep gradient region ($\rho _{\rm tor}=0.97$) with nominal profiles and tuned profiles that are designed to increase ETG drive but remain within experimental uncertainties at all locations (see figure 13 in Appendix A). At $\rho _{\rm tor}=0.97$ the tuned profiles have approximately 30 % increased 1/$L_{\rm Te}$ and 30 % decreased 1/$L_{\rm ne}$, resulting in a change of $\eta _e=L_{\rm ne}/L_{\rm Te}$ from approximately 1.5 to 3. To make a connection to recent ETG transport scalings for JET and DIII-D (Guttenfelder et al. Reference Guttenfelder, Groebner, Canik, Grierson, Belli and Candy2021; Chapman-Oplopoiou et al. Reference Chapman-Oplopoiou, Hatch, Field, Frassinetti, Hillesheim, Horvath, Maggi, Parisi, Roach, Saarelma and Walker2022; Field et al. Reference Field, Chapman-Oplopoiou, Connor, Frassinetti, Hatch, Roach and Saarelma2023), heat fluxes are shown in modified electron gyro-Bohm units (see figure caption) and megawatts. At nominal parameters the ETG heat flux is very small, whereas with tuned profiles ETGs drive substantial transport of about 0.9 MW. This indicates a strong stiffness of ETG transport in the pedestal centre. By artificially neglecting magnetic drifts in the simulation we can selectively turn off the drive of the toroidal ETG branch to asses the relative importance of slab and toroidal ETG for the heat transport. We find that the pure slab ETG simulation produces slightly less than half of the original transport. This shows that nonlinearly, both slab and toroidal branch contribute substantially to the ETG heat transport in this scenario.

Figure 10. Heat flux due to ETGs at $\rho _{\rm tor}=0.97$ as a function of time. Here $L_{\rm ref}=0.65$ m, $v_{\rm Te}=\sqrt {2T_e/m_e}$. Left axis in units of modified electron gyro-Bohm units $Q_{{\rm gB},e,{\rm mod}}=Q_{{\rm gB},e}\times 1/{L_{\rm Te}}^2$ following Chapman-Oplopoiou et al. (Reference Chapman-Oplopoiou, Hatch, Field, Frassinetti, Hillesheim, Horvath, Maggi, Parisi, Roach, Saarelma and Walker2022); right axis in units of megawatts. Parameter $Q_{{\rm gB},e,{\rm mod}}$ is evaluated at $\rho _{\rm tor}=0.97$ for the tuned profiles.

Figure 11 establishes a connection between the linear and nonlinear ETG simulations by comparing different cross-phases. As expected ETGs do not contribute to particle transport and $\phi \times n = \pi$ linearly and nonlinearly. In the potential temperature cross-phases $\phi \times T$ linear and nonlinear simulations show some resemblance but no clear agreement.

Figure 11. Cross-phases of electric potential $\phi$ and electron density fluctuations $n$ (a), parallel temperature component (b) and perpendicular temperature component (c) from nonlinear simulations (blue background) and linear simulations (orange circles).

5.4 Heat flux profile

To study how turbulent transport changes across the pedestal we consider the heat flux averaged over time and flux surface as a function of the radius. Figure 12 shows important aspects of the turbulent heat flux profile in the studied AUG pedestal. Particular care has to be taken with regard to the normalisation of the heat flux. Due to the strong temperature changes across the pedestal, the commonly used gyro-Bohm heat flux based on mixing-length estimates $Q_{\rm gB}=c_s p_i {\rho ^*}^2=(T_i/m_i)^{1/2}n_iT_i{\rho ^*}^2$ changes by two orders of magnitude across the pedestal (figure 12a). A modified gyro-Bohm heat flux $Q_{\rm gb,mod}=Q_{\rm gb}\times \max {(a/L_{\rm Ti},a/L_{\rm Te})}^2$, in which the minor radius is replaced by the gradient length as the macroscopic length scale, exhibits strong variations across the pedestal as well. We therefore focus the analysis on the heat flux in SI units.

Figure 12. Turbulent heat flux profile in an AUG pedestal from pedestal top to foot. (a) Gyro-Bohm heat flux profile. (b) Components of the ion-scale heat flux profile without $E\times B$ shear. (c) Total heat flux ($Q_{\rm es}$+$Q_{\rm em}$) due to ion-scale turbulence with and without $E\times B$ shear. (d) Total heat flux due to ion-scale turbulence from global simulations (red and blue solid lines) as well as ETG heat fluxes from local simulations at nominal values (light blue stars) and increased ETG (dark blue stars) compared with power balance calculations (broad lines). Region of increased measurement uncertainty in grey.

Figure 13. Tuned electron temperature and density profiles (dashed) used to test the sensitivity of ETG heat flux in comparison with nominal profiles (solid).

Figure 12(b) shows the profiles of the different heat flux components (electrostatic, electromagnetic, electron, ion). The dominant component over most of the pedestal is electrostatic electron heat flux (blue solid line), with the exception of $\rho _{\rm tor}\approx 0.97$, where the electrostatic ion heat flux has a local peak. This peak corresponds to the ITG/KBM mode identified in the linear analysis, which occurs at the peak of the ITG. Interestingly, the turbulent ion-scale heat flux is strongly reduced in all channels from pedestal top to pedestal centre and foot – even without $E\times B$ shear. The onset of this reduction coincides with the region of linear stabilisation, discussed in § 4.2.

Figure 12(c) illustrates the influence of $E\times B$ shear on the heat flux profiles. It reduces the ion-scale heat flux strongly and widens the region of almost vanishing turbulent heat flux.

Figure 12(d) shows that the turbulent heat flux on electron scales behaves oppositely across the pedestal compared with the ion-scale turbulent heat flux. It vanishes at the pedestal top and strongly increases down the pedestal (stars).

A comparison with power balance and neoclassic calculations shows qualitative agreement in heat flux structure and trends (figure 12d). The electron channel is dominant throughout the pedestal. Our gyrokinetic simulations reveal that the electron heat flux transitions in scale. At the pedestal top/shoulder it is driven by ion-scale TEM/MTM turbulence while at the pedestal centre and foot it is driven by small-scale ETG turbulence. While power balance considerations suggest a roughly constant electron heat flux across the pedestal, our gyrokinetic simulations find a reduction towards the pedestal centre. The discrepancy might, for example, be due to $T_i$ and $T_e$ profile uncertainties influencing simulations as well as the implicit uncertainties of the interpreted power balance results.

The region around $\rho _{\rm tor}=0.92$ where the total gyrokinetic electron heat flux reaches a minimum likely indicates a limitation of our separate scale ansatz. At this location heat flux is likely driven by scales that are resolved neither in our global simulations nor in the electron-scale simulations (cf. figure 4) and would require multi-scale simulations to be properly resolved. The ion channel contributes substantially to the total heat flux at the pedestal top/shoulder and reduces to neoclassic values towards the pedestal centre. At the pedestal foot our gyrokinetic simulations do not show the increase of ion heat flux suggested by power balance. The differences in the grey-shaded region are possibly due to increased measurement uncertainties in the ion profile that affect both power balance (see Viezzer et al. Reference Viezzer, Cavedon, Cano-Megias, Fable, Wolfrum, Cruz-Zabala, David, Dux, Fischer, Harrer, Laggner, McDermott, Plank, Pütterich and Willensdorfer2020) and gyrokinetic simulations.

For the power balance comparison, the turbulent ion heat flux component was estimated by subtracting neoclassic heat flux calculated with NEOART from the ASTRA result (both from Viezzer et al. Reference Viezzer, Cavedon, Cano-Megias, Fable, Wolfrum, Cruz-Zabala, David, Dux, Fischer, Harrer, Laggner, McDermott, Plank, Pütterich and Willensdorfer2020). The total ion heat flux due to power balance is constant (Viezzer et al. Reference Viezzer, Cavedon, Cano-Megias, Fable, Wolfrum, Cruz-Zabala, David, Dux, Fischer, Harrer, Laggner, McDermott, Plank, Pütterich and Willensdorfer2020) and the minimum in the turbulent heat flux is compensated by an increased neoclassic heat flux, which has a roughly constant diffusivity across the pedestal so that its heat flux follows the increasing ITG.

The following picture emerges for the investigated scenario. Turbulent heat flux at the pedestal top is dominated by electrostatic TEMs with electromagnetic MTM contributions. The interplay of low magnetic shear and increasing pressure gradient stabilises modes locally before the steep gradient region ($\rho _{\rm tor}\approx 0.94$) and reduces turbulent heat flux. The $E\times B$ shear suppresses turbulent heat flux strongly in all channels and widens the region of vanishing turbulent ion-scale heat flux. Turbulent electron heat flux changes scale across the pedestal and turbulent ion heat flux strongly reduces towards the pedestal centre.