2. Overview of argon-seeded EDA H-mode

ASDEX Upgrade discharge #36330 is an excellent candidate for gyrokinetic study, having well-diagnosed temperature and rotation profiles, and having previously been studied experimentally in Kallenbach et al. (Reference Kallenbach, Bernert, David, Dunne, Dux, Fable, Fischer, Gil, Görler and Janky2021). The discharge operates in a lower-single-null magnetic configuration with a toroidal magnetic field $B_T = -2.5\ {\rm T}$ on axis and a plasma current of 0.8 MA. The net power through the separatrix is feedback-controlled in real time via argon seeding, maintaining the absence of ELMs in EDA H-mode during auxiliary heating scans. The H-mode confinement phase is only lost upon the discontinuation of external heating. The stationary EDA H-mode period from 6.1–6.2 s with a constant and even mix of ECRH and NBI power totalling about 5 MW is chosen for the simulations. Furthermore, the QCM is experimentally present from around 20–40 kHz during this time phase in discharge 36330 (Kallenbach et al. Reference Kallenbach, Bernert, David, Dunne, Dux, Fable, Fischer, Gil, Görler and Janky2021). The regime presented in this work has higher core temperature and pressure compared with lower-power unseeded EDA H-mode discharges in ASDEX Upgrade. However, the pedestal profiles remain low despite higher input power, offering an explanation for the absence of ELMs.

Figure 1 displays experimental radial profile measurements available for this discharge together with cubic spline fits used as simulation inputs shown by a solid line in each plot. The radial coordinate along the top axis of figure 1, $\rho _{\rm tor}= \sqrt {(\varPhi _{\rm tor}-\varPhi _{\rm a})/(\varPhi _{{\rm tor},{\rm sep}} - \varPhi _{a})}$ is used throughout this work, where $\varPhi _{\rm tor}$ is the toroidal magnetic flux, and the subscripts ‘${a}$’ and ‘${\rm {sep}}$’ refer to the magnetic axis and separatrix, respectively. The profiles are mapped to the outer midplane using the IDE equilibrium (Fischer et al. Reference Fischer, Bock, Dunne, Fuchs, Giannone, Lackner, McCarthy, Poli, Preuss and Rampp2016) which considers the bootstrap current. The profiles are radially aligned, with shifts up to a few millimetres to account for certain equilibrium and measurement uncertainties, improving consistency between diagnostics.

Figure 1. Experimental radial profile measurements of ASDEX Upgrade discharge #36330 from 6.115–6.190 s, as a function of normalised poloidal (bottom) or toroidal (top) flux radius: (a) electron density from Thomson scattering, (b) electron and ion ($N^{7^{+}}$ impurity) temperature from Thomson scattering and CXRS, respectively, (c) impurity toroidal and poloidal rotation velocity from CXRS and (d) radial electric field from CXRS. Cubic spline fits used as GENE simulation inputs are shown as solid lines and local simulation locations are marked with dashed vertical lines. Spline knot locations in $\rho _{\rm {pol}}$ are: $a_{{\rm ne}}$, $[0.22,0.50,0.70,0.92,0.98,0.99,1.00,1.01,1.08]$; $b_{\rm Te}$, $[0.22,0.35,0.55,0.98,0.99,1.00,1.01,1.05]$; $b_{{\rm Ti}}$, $[0.00,0.50,0.70,0.98,0.99,1.01]$; $c_{{\rm vtor}}$, $[0.51,0.60,0.90,0.98,1.00]$; $c_{{\rm vpol}}$, $[0.93,0.97,0.99,1.00]$; and $d_{{\rm Er}}$, $[0.93,0.97,0.99,1.00]$.

Electron density and temperature are measured via Thomson scattering in the edge and core as shown in figures 1(a) and 1(b) and radially shifted in order for the separatrix temperature to be consistent with the two-point model of the divertor scrape-off layer (Stangeby Reference Stangeby2000). The density profile is cross-validated by line-integrated interferometry measurements, also used to recalibrate the core Thomson scattering measurements by $-10\,\%$.

The regime outlined in this work includes heating from a continuous, full-power and full-voltage (60 kV) NBI source, allowing for CXRS-based $N^{7^{+}}$ temperature and rotation measurements in the core and pedestal, from $\rho _{\rm {tor}}= 0.35 \unicode{x2013}1$, as shown in figure 1. Poloidal and toroidal impurity rotation are also measured by CXRS, but only the edge system is capable of measuring the poloidal rotation, as shown in figure 1(c). Ion temperatures are taken from core and edge CXRS systems as shown in figure 1(b) and radially aligned for the position of maximum gradient (Viezzer et al. Reference Viezzer, Pütterich, Conway, Dux, Happel, Fuchs, McDermott, Ryter, Sieglin and Suttrop2013) to match that of the electron temperature. The radial electric field, shown in figure 1(d), is computed from the radial force balance equation, $E_r = \boldsymbol {\nabla p_a} /( e Z_a n_a) - \boldsymbol{v_a \times B}$, where $\boldsymbol{p_a}$, $Z_a$, $n_a$ and $\boldsymbol{v_a}$ are the impurity pressure, charge number, density and velocity, respectively, $e$ is the elementary charge and $\boldsymbol{B}$ is the magnetic field. The impurity density profile of the fully stripped $N^{7+}$ population from CXRS measurements (McDermott et al. Reference McDermott, Dux, Pütterich, Geiger, Kappatou, Lebschy, Bruhn, Cavedon, Frank and den Harder2018) is included in the first term. The resulting ${\boldsymbol{E\times B}}$ shearing profile is used as input into the GENE simulations as discussed in detail in § 3, with the contribution from the poloidal velocity assumed to be negligible in the core.

Not included in figure 1 is a radially uniform effective ion charge estimate, $Z_{\rm {eff}}=2.05$, based on measurements of the background bremsstrahlung spectrum (Fischer et al. Reference Fischer, Hanson, Dose and von der Linden2000; Rathgeber et al. Reference Rathgeber, Fischer, Fietz, Hobirk, Kallenbach, Meister, Pütterich, Ryter, Tardini and Wolfrum2010) and used throughout this work. The $Z_{\rm {eff}}$ parameter primarily influences collisionality in two species simulations while in three species simulations this parameter is discarded in favor of more realistic impurity profiles, the inclusion of which also affects density gradients and profiles in addition to collisionality.

Finally, an argon impurity density profile is reconstructed via modelling in a similar manner to argon analyses in Dux et al. (Reference Dux, Cavedon, Kallenbach, McDermott and Vogel2020) and Kallenbach et al. (Reference Kallenbach, Bernert, David, Dunne, Dux, Fable, Fischer, Gil, Görler and Janky2021), using profiles from figure 1 as partial inputs to the impurity transport code STRAHL (Dux Reference Dux2006). The argon density profile, labelled Ar profile 1 in figure 2(a), is derived in the core from the soft x-ray (SXR) diagnostic after tomographic inversion of measured line-of-sight integrals (Odstrčil et al. Reference Odstrčil, Pütterich, Odstrčil, Gude, Igochine and Stroth2016). The contribution of tungsten and deuterium ions to the SXR emission is included in the calculation, being rather small when compared to the argon contribution in this discharge. In the pedestal where SXR emission is almost non-existent due to the lower temperature, the argon density is determined by fitting the modelled radiation, including charge exchange with neutral deuterium atoms, to spectral lines in the vacuum ultraviolet range measured by the SPRED spectrometer, and to the total radiated power measured by bolometry (David et al. Reference David, Bernert, Pütterich, Fuchs, Glöggler and Eich2021) and shown in figure 2(b). Previous analysis has clearly shown that argon seeding increases local radiation only slightly in the core, but significantly in the pedestal, contributing both to ELM avoidance and a colder divertor compared to an identical unseeded discharge.

Figure 2. (a) Two argon impurity profile fits used as GENE inputs plotted against $\rho _{\rm {tor}}$. Here $n_{\rm {Ar}}$ is the sum of the densities of all charged states of Ar. Profile 1 is derived from STRAHL modelling while profile 2 is based on a crude $Z_{\rm {eff}}$ estimate. (b) Experimental radiated, net, and total power flux integrated up to each flux surface.

A second crude argon profile estimate in figure 2(a), labelled Ar profile 2, is generated from assuming a uniform $Z_{\rm {eff}}=2.05$ and assuming argon as the only impurity while enforcing quasineutrality together with existing density profiles. This crude estimate is used as an alternate upper argon density estimate to further explore the effects of argon impurities. These argon profiles are used as impurity study inputs in §§ 4 and 5. GENE does not reproduce atomic radiation effects, however the net heat flux after excluding the experimentally measured radiated power, $P_{\rm {tot}} - P_{\rm {rad}}$, shown in figure 2(b), is used to compare simulated nonlinear heat fluxes with experimentally measured values in § 5. The total heat flux $P_{\rm {tot}}$ is computed with an interpretative TRANSP run (Pankin et al. Reference Pankin, McCune, Andre, Bateman and Kritz2004) that includes the local deposited heating power using the real geometry of the ECRH and NBI sources as well as profiles from figure 1 and argon profile 1.

3. Set-up of gyrokinetic modelling

Utilising the fits to the experimental data presented in figure 1, linear and nonlinear gyrokinetic computations are carried out locally at four radial positions and globally for one pedestal scenario using the GENE code. GENE employs a five-dimensional phase-space grid consisting of radial, binormal, field line following, parallel velocity (with respect to the magnetic field direction) and magnetic moment coordinates labelled as $x, y, z,v_{\parallel }$ and $\mu$, respectively. All simulations are local flux tube simulations self-consistently retaining $A_{||}$ and $B_{||}$ fluctuations, with the exception of the global gyrokinetic pedestal simulation, which is electrostatic due to computational limitations. Furthermore, all simulations use a realistic deuteron to electron mass ratio. All simulations use periodic boundary conditions, with the exception of the global simulations, which use Dirichlet boundary conditions with Krook-type (Görler et al. Reference Görler, Lapillonne, Brunner, Dannert, Jenko, Merz and Told2011) buffer zones in the outer $5\,\%$ of the simulation domain. All simulations are two species (deuteron and electron) with the exception of the linear impurity study in § 4 and the nonlinear case at $\rho _{\rm {tor}}=0.70$ in § 5.

Throughout this work the ‘Landau–Boltzmann’ collision operator is utilised (Doerk Reference Doerk2012) with the collisionality, $\nu _{{ei}}$, calculated (Neiser et al. Reference Neiser, Jenko, Carter, Schmitz, Told, Merlo, Bañón Navarro, Crandall, McKee and Yan2019) from the other input parameters, and with finite Larmor radius (FLR) effects not included. Preliminary tests superfluous to this work reveal that using the ‘Sugama’ collision operator (Sugama, Watanabe & Nunami Reference Sugama, Watanabe and Nunami2009) without FLR effects did not result in significant changes to linear growth-rates.

Linear growth rates are calculated for each wavenumber investigated, $k_y$, which is normalised by the respective gyroradius of the simulation $\rho _s = c_s/\varOmega _i$ with $c_s$ being the ion sound speed and $\varOmega _i$ being the ion cyclotron frequency. The gyroradius to tokamak minor radius is given as $\rho ^{*} = \rho _s/a$. Linear growth rates are normalised by $c_s/a$, where $c_s=\sqrt {T_e/m_{i}}$, $T_e$ is the electron temperature and $m_{i}$ is the ion mass. The total simulated heat flux for a given radial location is composed of contributions from two-species ion-scale, single-species electron-scale (adiabatic ions) and neoclassical estimates. Neoclassical calculations in GENE are based on a linearised solver described previously in Doerk (Reference Doerk2012), Oberparleiter (Reference Oberparleiter2015) and Stimmel et al. (Reference Stimmel, Bañón Navarro, Happel, Told, Görler, Wolfrum, Martin Collar, Fischer, Schneider and Jenko2019) and are here employing the Landau collision operator.

Important parameters used throughout this work which are taken from the experimental fits to the data are presented in table 1. Two core and two pedestal top positions are chosen for local simulations to allow for characterisation of core and pedestal top. Pedestal top and global simulation domain locations are limited radially by numerical prerequisites to ensure code stability (such as buffer zones) and experimental uncertainty. Ion and electron temperatures are presented in kiloelectronvolts, density $n_{\rm {ion}}$ is assumed to be equal to $n_e$ to enforce simulation quasineutrality in two species simulations, and all normalised logarithmic gradients are defined as $\omega _{{k}} = - ({1}/{k})({{\rm d}k}/{\rho _{\rm {tor}}})$ where $k$ represents the ion temperature $T_i$, electron temperature $T_e$ or the electron density $n_e$ in table 1.

Table 1. Nominal GENE input parameters taken from experimentally measured values from ASDEX Upgrade discharge #36330 from $t \sim 6.1 \unicode{x2013}6.2\ {\rm s}$ for four radial positions.

The thermal plasma pressure to magnetic pressure ratio is given as $\beta _e = {8 {\rm \pi}n_e T_e}/{B^{2}_{0}}$. The magnetic shear is defined as $\hat {s} = (\rho _{\rm {tor}}/q)({\rm d}q/{\rm d}\rho _{\rm {tor}})$ where $q$ is the tokamak plasma safety factor. The effective ion charge $Z_{\rm {eff}}$ is set to be 2.05 for all radial positions. In simulations where an impurity species profile is used, this parameter is ignored. The GENE ‘Landau–Boltzmann’ collision value, $\nu _{ei}$, is also listed in table 1 for reference.

Radial electric field shearing is included in nonlinear simulations by shifting the radial Fourier mode grid in time radially with aliasing compensation (Hammett et al. Reference Hammett, Dorland, Louriero and Tatsuno2006; McMillan, Ball & Brunner Reference McMillan, Ball and Brunner2019). The normalised local external shearing rate in GENE, $\gamma _{{E\times B}}$, is given as

(3.1)\begin{equation} \gamma_{{E\times B}} ={-}\sigma_{I_p} \sigma_{\omega} \sigma_{\partial_x |\omega|} \left| \frac{\rho_{\rm{tor,0}}}{q_0} \frac{\partial \omega_{\rm{tor}}}{\partial \rho_{\rm{tor}}} \right| \frac{a}{c_s}. \end{equation}

Here $\sigma _{I_p}$, $\sigma _{\varOmega }$ and $\sigma _{\partial _x |\varOmega |}$ are the sign of the plasma current, plasma rotation and radial derivative of the plasma rotation profile, respectively, with $+1/{-}1$ corresponding to clockwise/anticlockwise as seen from above. The reference flux surface, $\rho _{\rm {tor,0}}$, is taken at the simulation domain centre, $q_0$ is the safety factor at the simulation domain centre and $\omega _{\rm {tor}}$ is the toroidal angular velocity in radians per second taken from one of the two profiles displayed in figure 3.

Figure 3. Two background velocity profile fits used as GENE inputs for ASDEX Upgrade 36330. The blue $\omega _{\rm {tor}}$ fit represents the toroidal component of the total measured velocity. The dashed orange $\omega _{\rm {pseudo}}$ fit represents the $\boldsymbol{E\times B}$ velocity based on the radial electric field profile multiplied by the prefactor $|B_{\rm {Tot}}|/B_{\rm {pol}}$ for use via the Hammett–McMillan shearing mechanism in GENE.

Using the toroidal angular velocity profile to estimate the radial electric field is an appropriate estimate for core simulations. However, in the pedestal this assumption is no longer valid as the magnitude and sign of the $\boldsymbol{E\times B}$ velocity is no longer primarily represented by the toroidal component alone in the radial electric field force balance equation. The current implementation of $\boldsymbol{E\times B}$ shearing in the GENE code only considers a toroidal component, and so a pseudo-velocity, $\omega _{\rm {pseudo}}$, shown in figure 3, is used to replace $\omega _{\rm {tor}}$ in (3.1) as an ad-hoc representation for the full magnitude of the $\boldsymbol{E\times B}$ shearing in some nonlinear pedestal scenarios.

This pseudo-velocity is the full magnitude of the $\boldsymbol{E\times B}$ velocity calculated from the radial electric field shown in figure 1(d) (rather than just the toroidal component) multiplied by the prefactor $|B_{\rm {Tot}}|/B_{\rm {pol}}$ to seamlessly integrate with the toroidal $\boldsymbol{E\times B}$ shearing scheme already implemented in GENE. The resulting GENE $\gamma _{{E\times B}}$ shearing rates used in this work from both toroidal angular velocity profiles are listed in table 1. As is shown in detail in § 6, the utilisation of a pseudo-velocity to represent the full effect of the radial electric field shearing elucidates the importance of $\boldsymbol{E\times B}$ shearing in the tokamak pedestal, however it also demonstrates the limit of local GENE simulations to accurately represent shearing, and ultimately serves as a motivation for a more complete numerical methodology from the full radial electric field vector.

4. Linear instability study

Linear characterisation of the ELM-free regime presented in this paper reveals a core scenario similar to standard H-Mode, and a more complex picture of gyrokinetically modelled pedestal top turbulence. All linear simulations presented here are tested for numerical convergence by increasing spatial resolutions, box sizes and by varying the edge$\_$opt parameter, used to redistribute GENE parallel grid points to better adhere to tokamak magnetic geometry (Told Reference Told2012), until a growth-rate difference of less than 5 % or a difference of less than 0.01 $k_y\rho _s$ is achieved. Linear simulations utilised a resolution range of $N_x \times N_z \times N_{v\parallel } \times N_\mu = 31\unicode{x2013}127 \times 48\unicode{x2013}128 \times 32\unicode{x2013}48 \times 16$ where $N_x,N_z,N_{v\parallel }$ and $N_\mu$ are the number of radial, parallel coordinate, parallel velocity and magnetic moment grid points, respectively. Generally simulations with core parameters at ion scale $k_y\rho _s$ tend to require lower numeric resolutions while parameters at $k_y\rho _s < 0.2$ or at $k_y\rho _s > 10$ require additional resolution. Linear velocity and magnetic moment box size parameters of $L_v \times L_\mu = 3 \times 9 \unicode{x2013} 11$ are appropriate for all linear simulations. Values of ${\tt {edge\_opt}}= 1\unicode{x2013}2$ suffice for all simulations with the exception of $k_y\rho _s > 100$ calculations at $\rho _{\rm tor}= 0.95$ utilising ${\tt {edge\_opt}} = 4$. Linear simulations are carried out with deuteron and electron species unless noted otherwise and results shown are evaluated at a ballooning angle of zero.

A linear, local, initial-value scan across a wide range of $k_y \rho _s$ for four radial positions with two species gives a first characterisation of the regime as shown in figure 4.

Figure 4. Linear GENE growth rates ($\gamma$) and frequencies ($\omega$) for four radial positions across a range of $k_y\rho _s$. Negative frequencies in the GENE simulation co-moving frame correspond to the electron diamagnetic direction. A subplot shows ion-scale frequencies on a smaller scale for visibility. Frequencies from negative or zero growth rate $\gamma$ are omitted.

An ion-scale instability peak is present in all four cases at $k_y\rho _s < 1$ while an electron scale peak exists for core scenarios at $k_y\rho _s > 3$ and for pedestal top scenarios at $k_y\rho _s > 60$. When considering the heuristic rule (Görler & Jenko Reference Görler and Jenko2008a; Görler & Jenko Reference Görler and Jenko2008b; Jenko Reference Jenko2004) often applied to multiscale simulations of $\gamma _e/ \gamma _i \gg \sqrt {(m_i/m_e)} \sim 60$, the two core scenarios as well as the pedestal top scenario do not exceed $\gamma _e/ \gamma _i \gg \sim 60$ and, therefore, multiscale simulations are assumed to not be important. The $\rho _{\rm {tor}}=0.95$ case does exceed this value ($167 > 60$). However, all four scenarios do not meet the criteria for multiscale importance using ${\rm {max}}(\gamma _{\rm {ETG}}/k_y) \geqslant {\rm {max}}(\gamma _{\rm {ITG}}/k_y)$ from Staebler et al. (Reference Staebler, Howard, Candy and Holland2017) and the ion and electron scale peaks are separated by more than two orders of magnitude; therefore, the assumption that these two peaks do not interact nonlinearly is made for this work.

Two intermediate range peaks from $k_y\rho _s=2\unicode{x2013}10$ with low growth rate relative to the respective ion-scale growth-rate peaks are also present for the two pedestal top scenarios. Gradient variation test scans (not shown) are performed for all four radial locations to help identify the dominant instability.

The two core scenarios are easily identifiable as ion temperature gradient (ITG) dominated at low $k_y\rho _s$ due to positive growth-rate frequencies and a sensitivity to varying $\omega _{{\rm Ti}}$ at low $k_y\rho _s$. The electron-scale core peaks are identified as electron temperature gradient (ETG) dominated due to negative growth-rate frequencies and a strong sensitivity to variations in $\omega _{{\rm Te}}$ at the upper $k_y\rho _s$ range.

Linear pedestal calculations present a more complex picture of the argon-seeded EDA H-mode. Negative-frequency instabilities sensitive to variation primarily in $\omega _{{\rm Te}}$ but also $\omega _{{\rm Ti}}$ and $\omega _{n}$ are found at ion-scale $k_y\rho _s$ for $\rho _{\rm tor}= 0.90$, with a positive frequency mode beginning to dominate at $k_y\rho _s = 1$, hinting at TEM transitioning to ITG as the dominant instabilities as $k_y\rho _s$ increases (Kotschenreuther et al. Reference Kotschenreuther, Liu, Hatch, Mahajan, Zheng, Diallo, Groebner, Hillesheim, Maggi and Giroud2019). Although linear flux values are not meaningful on their own, by comparing the ratios of finite linear flux channels in converged linear simulations, it is possible to discern which flux channel is dominant to help identify an instability. These modes at $\rho _{\rm tor}= 0.90$ are dominated by electron electrostatic heat flux channel when considering the ratio of electrostatic to electromagnetic heat flux channel by species.

The ion-scale peak at $\rho _{\rm tor}= 0.95$ is characterised by negative frequencies and electron electrostatic dominated heat flux with the exception of the lowest positive growth rate at $k_y\rho _s = 0.1$ matching a microtearing mode (MTM) ballooning mode structure like in Hatch et al. (Reference Hatch, Told, Jenko, Doerk, Dunne, Wolfrum, Viezzer and Pueschel2015) for example, and electron electromagnetic-dominated heat flux.

Upon decreasing the density gradient, the growth rate increases by roughly a factor of two in the ion-scale and intermediate ranges. Finally, the electron scale peak at $\rho _{\rm tor}= 0.95$ is ETG-like, and is subject to an increased growth rate upon decreasing the density gradient as expected (Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Merlo, Field, Giroud, Hillesheim, Maggi, von Thun and Roach2019; Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020). While it has been previously demonstrated that ETG growth-rates in particular may peak at non-zero ballooning angles due to strong equilibrium shaping effects (Told et al. Reference Told, Jenko, Xanthopoulos, Horton and Wolfrum2008), the focus in this section is limited to the growth rate at the outboard midplane for instability identification and comparison with nonlinear results; nonlinearly it is found throughout this work that heat flux spectra peak at or near the outboard midplane.

Additional scans over $\beta _e$ and collisionality are included for the two pedestal top scenarios and are displayed in figures 5 and 6, respectively.

Figure 5. GENE linear $\beta _e$ scan at (a) $\rho _{\rm {tor}}= 0.90$ and (b) $\rho _{\rm {tor}}= 0.95$. The black vertical line denotes the nominal $\beta _e$ value for each scenario.

Figure 6. GENE linear collisionality scan for two pedestal top positions at $\rho _{\rm {tor}} = 0.90$ and $0.95$. Three values of collisionality are tested for each $k_y\rho _s$.

Unlike some previous recent GENE pedestal studies (Bonanomi et al. Reference Bonanomi, Angioni, Crandall, Di Siena, Maggi and Schneider2019; Stimmel et al. Reference Stimmel, Bañón Navarro, Happel, Told, Görler, Wolfrum, Martin Collar, Fischer, Schneider and Jenko2019), scanning along increasing $\beta _e$ with fixed geometry does not reveal instabilities which are bordering a highly unstable kinetic ballooning mode (KBM)-like range at $\rho _{\rm tor}= 0.90$ as shown in figure 5. While a critical threshold is observed at higher $\beta _e$ at $\rho _{\rm tor}= 0.95$, the proximity of the nominal $\beta _e$ is at least a factor of two away from increased linear growth rates which are indicative of a KBM in both pedestal top scenarios in figure 5.

Collisionality scans at the ion-scale pedestal top scenario peaks with $50\,\%$ and $150\,\%$ the nominal collisionality value suggest again trapped modes which are sensitive to collisionality as shown in figure 6(a). Linear frequencies are plotted in figure 6(b) for reference, illustrating the effects of collisionality on linear frequencies. Together with the frequency of the ion-scale peaks in figure 4, one can deduce that at $\rho _{\rm {tor}} = 0.90$ the dominant instability transitions from TEM to ITG, while at $\rho _{\rm {tor}} = 0.95$ the ion-scale peak is TEM with a very weak MTM detected at the lowest positive growth rate $k_y\rho _s < 1$.

A scan across $k_y\rho _s$ performed with two different argon impurity profiles is shown in figure 7, together with a linear scan using only electrons with adiabatic ions. Introducing a third argon species into the linear core simulations is found to have only a minimal effect on dominant growth rates at ion-scale wavenumbers and decreases electron-scale growth-rate peaks significantly relative to two species linear calculations, thus rendering significant cross-scale coupling more unlikely. While the decrease in electron scale growth rate from the addition of argon agrees with expectations, the nonlinear electron-scale heat flux contribution is minimal as shown in § 5. The relatively small changes in ion-scale behaviour from the addition of argon weights the assumption of using two species for nonlinear simulations. Furthermore, it should be noted that using the artificially high Ar profile 2 is shown to have a relatively greater effect nonlinearly than linearly at $\rho _{\rm {tor}}=0.70$ as shown in § 5.

Figure 7. Linear GENE argon impurity study at four radial locations. (a)–(d) The converged linear growth ratesfor no impurity, impurity profile 1, impurity profile 2 and adiabatic ion linear scans for $\rho _{\rm {tor}}= 0.60,0.70,0.90$ and $0.95$, respectively.

In addition, linear scans at $\rho _{\rm {tor}}=0.70$ over $k_y\rho _s$ with the inclusion of argon profiles 1 and 2 are reproduced with a newly implemented dilution model in GENE. The new dilution model allows for the modification of density and gradient profiles without the direct introduction of a third species; the dilution model is used to linearly reproduce the density dilution and gradient modification introduced by the addition of argon profiles 1 and 2 as shown in figure 8.

Figure 8. Linear GENE growth rate comparison between three-species and dilution-model $k_y\rho _s$ scans. The default two-species linear scan is shown in black, scans with argon profiles 1 and 2 are shown in red and magenta, respectively, and the dilution model scans are shown in blue and green, respectively.

Figure 8 illustrates that at ion scales the linear growth rate peak and shape are closely reproduced by the dilution model for both argon profiles for ion scale at $\rho _{\rm {tor}}=0.70$, demonstrating the significance of gradient modification and density dilution from impurities. However, as shown in § 5, considering that the nonlinear heat flux shift is greater than the linear growth rate shift upon the addition of argon, one must not rule out the nonlinear role the addition of argon plays on turbulent heat transport.

In both cases linear simulations qualitatively predict the increased heat flux trend when comparing between two-species and three-species cases.

5. Nonlinear GENE simulations

Using linear analysis as an initial guide for the characterisation of the argon-seeded EDA H-mode, nonlinear gyrokinetic results are presented here. Nonlinear heat fluxes for various nonlinear simulations as well as neoclassical estimates are given in table 2 together with experimentally measured heat flux $P_{\rm {tot}} - P_{\rm {rad}}$ from figure 2(b).

Table 2. Simulated and experimental heat flux values in megawatts. From left to right the columns are radial position, optional modification to the ion-scale simulation, specification of either standard toroidal angular velocity or pseudo-angular velocity profile input used for determination of $\gamma _{E\times B}$, ion-scale heat flux subdivided into four channels (six with impurities), electron-scale heat flux sum where asterisks denote three-species simulations, neoclassical heat flux estimate, sum of total heat flux as an upper estimate with statistical error included and the net experimental heat flux value defined as $P_{\rm {tot}} - P_{\rm {rad}}$.

5.5. Global pedestal simulations

Global effects are confirmed to be of importance in the pedestal regime as evidenced by two global, gradient-driven, electrostatic, ion-scale simulations Dirichlet boundary conditions that reduce local overpredictions of the experimental heat flux value at $\rho _{\rm tor}=0.90$, with the input $\gamma _{E\times B}$ profile proving crucial. Using past gyrokinetic findings such as those of Candy et al. (Reference Candy, Holland, Waltz, Fahey and Belli2009) and Görler et al. (Reference Görler, Lapillonne, Brunner, Dannert, Jenko, Aghdam, Marcus, McMillan, Merz and Sauter2011), one might expect that $\rho ^{*} = 0.0018 < 0.002$ is satisfied at $\rho _{\rm {tor}}=0.90$ and that flux tube simulations would agree with global simulations. Considering that more recent findings show this scaling is not necessarily valid for medium-sized tokamaks, however (Navarro et al. Reference Navarro, Told, Jenko, Görler and Happel2016), and considering that the finite width of the steep gradient region is of greater importance than the overall system size (McMillan et al. Reference McMillan, Lapillonne, Brunner, Villard, Jolliet, Bottino, Görler and Jenko2010), it is not unexpected that global simulations correct a local overprediction at $\rho _{\rm {tor}}=0.90$.

Local simulations with electromagnetic effects offer additional insight to how the inclusion of realistic $\beta _e$ may affect the heat flux. Global investigation is motivated by the local ion-scale nonlinear simulation at $\rho _{\rm tor}=0.90$ which yields a heat flux more than an order of magnitude higher than the experimental value. Such an overprediction in the pedestal has previously been associated with proximity to a nonlinearly shifted critical $\beta _e$ KBM threshold; such cases typically demonstrate a clear increase in linear growth rates with increasing $\beta _e$ and a strong heat flux sensitivity to $\beta _e$ nonlinearly (Bonanomi et al. Reference Bonanomi, Angioni, Crandall, Di Siena, Maggi and Schneider2019; Stimmel et al. Reference Stimmel, Bañón Navarro, Happel, Told, Görler, Wolfrum, Martin Collar, Fischer, Schneider and Jenko2019). However, in this case even setting $\beta _e$ to zero together with a high estimate of $\gamma _{E\times B}$ does not recover the experimental heat flux, suggesting the strong relevance of finite-size effects.

To investigate this, two global, electrostatic, nonlinear simulations with two different toroidal velocity inputs used to calculate $\boldsymbol{E \times B}$ shearing are presented in table 2. Global electromagnetic GENE simulations, having recently been implemented, would be of high value but are left for future work due to the high computational cost.

The importance of the $\gamma _{E\times B}$ profile is explored globally and comparisons are drawn between local and global simulations. Using only the toroidal velocity component to calculate $\boldsymbol{E \times B}$ shear gives a total heat flux of about 3.7 MW, which is a slight overestimate of pedestal heat transport after factoring in radiative losses. On the other hand, the pseudo-velocity effectively damps all turbulence in the pedestal, not even reaching the experimental 2–3 MW range for the pedestal as shown in figure 9. Global effects are necessary to correct the local heat flux overprediction at $\rho _{\rm {tor}}=0.90$; compared with the 6.91 MW from the local electrostatic simulation at $\rho _{\rm {tor}}=0.90$, the global simulation using $\omega _{\rm {psu}}$ as GENE input yields less than 0.1 MW at a local slice around $\rho _{\rm {tor}}=0.90$, as shown in figure 9. Furthermore, the choice of input rotational velocity is shown to have a dramatic effect as shown in figure 9. While the heat flux is not perfectly recovered by global simulations, even these electrostatic cases demonstrate the importance of global effects in the pedestal. Based on the comparisons between nonlinear electrostatic and electromagnetic local simulations shown here alone, it is likely that more physical realism via inclusion of electromagnetic effects may help corroborate the underprediction of heat fluxes when using $\omega _{\rm {psu}}$ in a global pedestal simulation.

Figure 9. Global electrostatic GENE heat flux profiles for two simulations with different background $\gamma _{E\times B}$ profiles, and one electrostatic local simulation for comparison. Global Dirichlet boundary buffer zones are denoted with grey patches. The experimental net power flux is plotted as a dashed magenta line for reference.

Of key importance to the argon-seeded EDA H-mode on ASDEX Upgrade is understanding the QCM. Despite extensive efforts, notable challenges were encountered while attempting to elucidate more about the QCM. Previous work from LaBombard et al. (Reference LaBombard, Golfinopoulos, Terry, Brunner, Davis, Greenwald and Hughes2014) shows the EDA H-mode spans the separatrix and into the steep gradient region of the pedestal; by considering this alone, the QCM would require a simulation which crosses the separatrix to fully capture the necessary physics. However, it is important to note that there is no evidence that the QCM does not exist further inside the pedestal top in ASDEX Upgrade ELM Free H-mode plasmas due to Doppler reflectometry limitations at higher densities inside the pedestal. Furthermore, evidence from Theiler et al. (Reference Theiler, Terry, Edlund, Cziegler, Churchill, Hughes, LaBombard and Golfinopoulos2017) also suggests that the QCM resides inside the $\boldsymbol{E_r}$ well minimum rather than outside.

Assuming that the QCM could be found at the pedestal top, the biggest challenge in such an analysis is that the local pedestal top heat flux discrepancy between experiment and simulation may be too large to make any meaningful frequency comparison. Finally, even if there was no such discrepancy, analysis of the nonlinear local simulation frequencies using the more realistic $\omega _{{\rm psu}}$ at $\rho _{{\rm tor}}=0.90$ reveals that different velocity profiles simultaneously alter resulting shearing as well as the applied Doppler shift to the simulations, resulting in a wide range of outcomes from even small changes in the velocity profile. Using the more realistic, higher-valued $\omega _{{\rm psu}}$ velocity profile value in conjunction with higher gradient inputs at $\rho _{{\rm tor}}=0.95$ requires significantly more numerical resolution in the radial direction than using the $\omega _{{\rm tor}}$ profile for $\boldsymbol{E\times B}$ input, thus creating an impractical computational challenge. This stems from the current local implementation of $\boldsymbol{E\times B}$ shearing in GENE which requires very high radial resolution for intense shearing rates to prevent artificial discontinuities from forming.

Conjecture on nonlinear global electrostatic frequencies is omitted from the electrostatic global analysis so as to not misrepresent the electromagnetic QCM frequencies (LaBombard et al. Reference LaBombard, Golfinopoulos, Terry, Brunner, Davis, Greenwald and Hughes2014).

5.6. Quasilinearity analysis

When characterising such a multidimensional parameter space it useful to determine which linear features persist nonlinearly. Considering the four scenarios, there is fair quasilinear agreement at peak $k_y\rho _s$ for the two core scenarios while the linear and nonlinear pedestal cross-phases do not align. Figure 10 displays potential perpendicular temperature fluctuation and potential density fluctuation cross-phases using the same definition as in Told et al. (Reference Told, Jenko, Görler, Casson and Fable2013) at $\rho _{\rm {tor}}=0.95$. Linear cross-phase peaks denoted with crosses do not match the nonlinear cross-phase peaks denoted with small circles at peak heat flux $k_y\rho _s$ shown in orange, nor do secondary linear cross-phase peaks (not shown here) for either $\phi \times n$ or $\phi \times T_{\perp }$.

Figure 10. Local nonlinear (NL.) and linear (lin.) simulation maximum cross-phase (a) ion and (b) electron amplitudes for $\phi \times n$ and $\phi \times T_{\perp }$ at $\rho _{\rm {tor}}=0.95$ for AUG discharge #36330 from 6–6.2 s. Nonlinear data are shown as circles while linear data are displayed as crosses. Nonlinear time-averaged electrostatic (ES) and electromagnetic (EM) heat flux spectra are overlaid as solid and dashed orange lines, respectively, for reference.

Linear frequency computations at $\rho _{\rm {tor}}=0.90$ qualitatively predict TEM/ITG nonlinearly, yet nonlinear frequencies are generally sensitive to the input toroidal velocity profile. Both nonlinear scenarios at this radial position are identified as ITG/TEM from the electron to ion diffusivity ratios (Kotschenreuther et al. Reference Kotschenreuther, Liu, Hatch, Mahajan, Zheng, Diallo, Groebner, Hillesheim, Maggi and Giroud2019) $\chi _e/ \chi _e \sim 1$ and $D_e/(\chi _i + \chi _e) \sim 0.16$ where $\chi _{i,e}$ is the ion or electron heat diffusivity and $D_e$ is the electron particle diffusivity. Using the stronger and more realistic $\omega _{\rm {psu}}$ profile as shearing input results in an increasingly negative-frequency (in the plasma frame) TEM-dominated frequency spectrum compared with $\omega _{\rm {tor}}$ at $\rho _{\rm {tor}}=0.90$.

Other recent research suggests that quasilinear agreement is attainable in lower-power scenarios (Bonanomi et al. Reference Bonanomi, Angioni, Crandall, Di Siena, Maggi and Schneider2019; Neiser et al. Reference Neiser, Jenko, Carter, Schmitz, Told, Merlo, Bañón Navarro, Crandall, McKee and Yan2019; Stimmel et al. Reference Stimmel, Bañón Navarro, Happel, Told, Görler, Wolfrum, Martin Collar, Fischer, Schneider and Jenko2019) such as L- or I-mode, and even to some degree in an H-mode inter-ELM pedestal scenario (Hatch et al. Reference Hatch, Told, Jenko, Doerk, Dunne, Wolfrum, Viezzer and Pueschel2015). However, these findings do not seem to translate to the pedestal top scenario here. Further considerations on quasilinear validity have previously been outlined in Casati (Reference Casati2012), and nonlinear effects for edge turbulence are detailed in Scott (Reference Scott2002).